Module iharperella 

HARPERELLA and iHARPERELLA: integrated phase unwrapping and background field removal

Both methods remove the background field through the Laplacian. The background field is harmonic inside the brain, so the Laplacian there holds only tissue sources. They estimate the Laplacian outside the brain, add it to the brain’s, and invert the Laplacian over the whole field of view.

iHARPERELLA is a port of STI Suite 3.0’s iHARPERELLA (its protected P-code was recovered by black-box probing: logging the calls it makes to its helper functions and to MATLAB’s lsqr, and probing the lsqr operator it builds). HARPERELLA is not in STI Suite 3.0; it follows the paper, and uses STI’s building blocks for what the paper leaves open (see harperella).

Building blocks (STI Suite)

• Laplacian: spectral, L(x) = F⁻¹[k² F(x)] with k in cycles/mm (STI’s fftnc grid), which is −∇²x / 4π². Its inverse divides by k² and sets the DC term to zero. The scale cancels between the two, so the result is in the input’s units.
• Boundary rim: brain voxels with more than one voxel outside the brain among the 146 neighbours of STI’s ball3D(3) (the offsets with 1 ≤ |d|² ≤ 10 voxels²), i.e. SMVFiltering(mask, 3) < 0.99. Within the rim, the 20 % of Laplacian values largest in magnitude are set to zero (the paper’s “large Laplacian values near the boundary of the brain … were considered inaccurate and set to zero”). Outside the brain the Laplacian is zero.
• Solver: MATLAB’s lsqr (no preconditioner, zero initial guess) at tolerance tol for at most max_iter iterations; STI passes 0.001 and niter (default 40).
• Padding: pad zero voxels on each side (STI’s padsize, default 0).

Past about ten iterations LSQR on these operators amplifies rounding: two implementations whose residuals agree to 1e-14 for the first ten iterations differ by ~1e-3 by the twelfth. So results agree with STI Suite to rounding for a few iterations and to a few percent at 40.

Input

• HarperellaInput::FieldPpm: an unwrapped total field in ppm, e.g. a multi-echo field map. This is what STI’s iHARPERELLA does with its input: it takes the spectral Laplacian of it as it is (it does not unwrap).
• HarperellaInput::WrappedPhase: one wrapped echo in rad, as in the HARPERELLA paper: the Laplacian is cos φ L(sin φ) − sin φ L(cos φ) (the paper’s Eq. 1) with the spectral L, and the result is divided by rad_per_ppm(te_s, b0_t).

Either way the output is the local field in ppm, inside the mask.

Reference pipeline

• Input: STI’s help text and STI_Templates.m call it on the “raw phase”, but probed in MATLAB R2023b, STI 3.0’s iHARPERELLA does not unwrap: on a synthetic head whose phase wraps in 59 % of the brain, its output from the wrapped phase correlates −0.03 with the true local phase, from the unwrapped phase 0.98. Hence HarperellaInput::FieldPpm.
• Defaults: niter 40; padsize 0 (probed: bitwise equal to [0 0 0], while STI_Templates.m lists [12 12 12]). IharperellaParams::default matches the function.
• Chain: STI Suite 2.01’s manual has the GUI combine the Laplacians of all echoes (“Combine Laplacian”) and then run iHARPERELLA, as the alternative to “Unwrapping” + V-SHARP, before iLSQR; the mask is the GUI’s magnitude threshold of the 3rd or 4th echo. STI Suite 3.0’s GUI has no iHARPERELLA button (only the function).

References:

• HARPERELLA: Li, W., Avram, A. V., Wu, B., Xiao, X., Liu, C. (2014). “Integrated Laplacian-based phase unwrapping and background phase removal for quantitative susceptibility mapping.” NMR in Biomedicine, 27(2):219-227. doi:10.1002/nbm.3056
• iHARPERELLA: Li, W., Wu, B., Liu, C. (2015). “iHARPERELLA: an improved method for integrated 3D phase unwrapping and background phase removal.” Proc. ISMRM 23, p.3313. Code: STI Suite 3.0, Core_Functions_P/iHARPERELLA.p.

Structs

HarperellaParams

HARPERELLA parameters.

IharperellaParams

iHARPERELLA parameters (STI Suite’s iHARPERELLA options). It has no SMV kernel.

Work 🔒

The padded working volume: data, mask, Laplacian (rim cleaned), and the padded grid.

Enums

HarperellaInput

What HARPERELLA / iHARPERELLA start from. See the module docs.

Functions

apply_inv_kernel 🔒

F⁻¹[F(x) / k] with the k = 0 term set to zero, real part.

apply_kernel 🔒

F⁻¹[k F(x)], real part.

apply_smv 🔒

Apply SMV filter in k-space: S(x) = ifft(s_fft * fft(x))

compute_smv_fft 🔒

Compute FFT of SMV kernel (returns real part, since kernel is symmetric)

drop_large_rim_values 🔒

Set the 20 % of rim Laplacian values largest in magnitude to zero: keep |L| ≤ the sorted rim magnitude at 1-based index round(0.8 N).

erode 🔒

Erode mask by a sphere of radius mm: keep the voxels whose whole sphere lies in the mask.

harperella

HARPERELLA: the exterior Laplacian makes the SMV of the total Laplacian uniform.

iharperella

iHARPERELLA: the exterior Laplacian minimises the weighted resulting field in the brain.

k2_kernel 🔒

STI’s spectral Laplacian kernel k² = kx² + ky² + kz², k in cycles/mm, in FFT order.

pad_centred 🔒

Pad with pad zero voxels on each side (STI’s padarray(x, padsize)).

shift 🔒

Circular shift by s voxels.

sti_rim 🔒

STI’s boundary rim: mask voxels where SMVFiltering(mask, 3), the mean of the mask over the 146 neighbours of ball3D(3) (offsets with 1 ≤ |d|² ≤ 10 voxels², periodic), is below 0.99, i.e. more than one of those neighbours is outside the mask.

to_complex 🔒

unpad_centred 🔒

Inverse of pad_centred.