Expand description
Scan geometry from the NIfTI affine: B0 direction, obliquity, and resampling to an axial grid.
The dipole kernel lives in the voxel grid, so an oblique acquisition must either supply the
true B0 direction or be resampled to a cardinal-aligned grid. Wrapped phase has to be
resampled in the complex domain — see geometry::resample_complex_to_axial.
Cropping reconstruction to the region that carries signal, and putting the answer back.
FFT-based stages cost O(N log N) in the whole grid, not in the brain. See crop for why
the box is also rounded up to FFT-friendly sizes, and for the wrap-around caveat.
Cropping a volume to the region that actually carries signal, and putting it back.
Every FFT-based stage — background removal, dipole inversion, Laplacian unwrapping — costs
O(N log N) in the whole grid, not in the brain. A brain typically fills a fifth of an
acquired volume and a tenth of one that has been resampled to a cardinal grid (resampling
wraps a box around a tilted slab, so half the result is empty corner). Reconstructing inside a
box around the mask and padding the answer back afterwards is therefore close to free, and on
a UK Biobank SWI it takes the resampled grid from 4.5× the cost of the acquired grid down to
parity.
Two things make the difference, and they are worth separating because their risk differs:
- Fewer voxels ([
crop_box_for_mask]). Mask bounding box plus a margin, rather than the full field of view. This moves the FFT’s periodic boundary closer to the object. - Sizes an FFT likes ([
fft_pad_box]). An extent with a large prime factor pushesrustfftonto Bluestein’s algorithm. An axially-resampled UK Biobank grid comes out 272×339×77 — that is 2⁴·17, 3·113 and 7·11, awkward on every axis, which is what resampling to a bounding box tends to produce. Padding to 280×343×80 costs 8% more voxels and takes the transform from 131 ms to 71 ms, a 1.85× speedup. Padding moves the boundary further from the object, so unlike cropping it carries no wrap-around risk.
Both are expressed as a [CropBox], which may sit inside the grid (cropping), extend beyond
it (padding), or do both on different axes.
§This changes the numbers, not just the speed
This caveat applies to cropping, not to padding. FFT-based reconstruction is periodic, so
moving the boundary closer to the object brings wrap-around with it, and the dipole kernel has
infinite support. Measured on a UK Biobank acquisition, a crop that actually removed voxels
changed χ by ~0.6% of its dynamic range at the median and ~4% at the 99th percentile.
[margin_voxels] takes the margin in millimetres so anisotropic voxels get a
geometrically equal margin on every side, but a caller should still validate a cropped
reconstruction against an uncropped one rather than assume the two agree.
[fft_pad_box] has no such caveat: it discards nothing and only moves the boundary outward.
Structs§
- CropBox
- A box within a larger grid: where reconstruction actually happens.
Functions§
- crop_
box_ for_ mask - The box to reconstruct in: the mask’s bounding box, grown by
margin_mmon every side, each axis then rounded up to an FFT-friendly size and clamped to the grid. - crop_
volume - Copy the box out of a full-grid volume.
- crop_
volume_ with - Copy the box out of a full-grid volume, filling anything outside the grid with
fill. - fft_
pad_ box - A box covering the whole grid, each axis grown outward to an FFT-friendly size.
- margin_
voxels - Margin in voxels per axis for a margin given in millimetres, at least one voxel where the margin is positive. Anisotropic voxels get a geometrically equal margin rather than an equal voxel count — 8 voxels is 6.4 mm in-plane but 24 mm through-plane at 0.8 × 0.8 × 3 mm.
- next_
fft_ friendly_ size - Smallest size
>= nwhose prime factors are all at most 7 — the radicesrustffthas dedicated butterflies for. Sizes with a large prime factor fall back to Bluestein’s algorithm and cost several times more despite holding the same data. - uncrop_
volume - Put a cropped volume back into a full-grid volume, filling everything outside the box with
fill(zero for field maps and χ, which are undefined outside the mask anyway).