Fit nlinv maps on a low-resolution grid, with optional phase-pole correction

[mcencini]

What failed

• Shared phase poles. One set of nlinv maps (apps.nlinv_maps, nlinv -m 1) can share a phase vortex. The maps are smooth, so every coil's sensitivity is driven to zero there, and any combination over them goes dark: the adjoint, or a pics solve.
• Full-resolution Cartesian fit. The Cartesian fit zero-filled the central size lines into the full matrix, so nlinv and the pole detection ran at full resolution: 14 s for a 220x220 slice with 32 coils. BART's pole circle is a twentieth of the grid's longest axis, rounded up, so on that grid it refused any slab of up to 4 partitions on a 64-point plane.
• Full-resolution non-Cartesian fit. The non-Cartesian fit zeroed the samples beyond radius but kept their trajectory. nlinv -t therefore sized its grid from the whole trajectory, ran at full resolution, and took the zeroed samples for measured zeros, which pulls the maps away from the coils. Three-dimensional trajectories were not supported at all: the reshape assumed a plane.

Change

• Phase-pole option. apps.nlinv_maps(..., phase_pole=None): None keeps BART's behaviour; an int d is nlinv --phase-pole d, detection after Gauss-Newton step d, with 0 meaning every step. apps.nlinv_pics forwards it.
• Cartesian. The central min(size, n) lines of every axis are cropped (centre at n // 2), and nlinv runs on that grid.
• Non-Cartesian.
  • Only the samples within radius are kept, and nlinv -t runs on an explicit grid of 2 ceil(radius) points per axis (-x, which BART honours).
  • Two- and three-dimensional trajectories are both supported.
  • New shape argument: the image grid of the maps, (y, x) or (z, y, x). By default it is the grid BART estimates from the whole trajectory, as before.
• Common to both paths.
  • The maps are zero-filled in k-space to the image grid and normalised to unit root sum of squares. Output shape and normalisation are unchanged.
  • BART's Sobolev weight (1 + a |k/n|^2)^(-b/2) is relative to the grid, so a is scaled by (m/n)^2 along the shortest cropped axis (m points kept of the axis's n).
  • With BART's own a, small-grid Cartesian maps came out much smoother than before: 0.38 from parallel to the full-grid maps at the worst voxel of a 64-point phantom.
• Thin grids. BART still refuses the detection on a grid with an axis no longer than a twentieth of its longest, rounded up. On the calibration grid at the default size=24, that means a slab of fewer than 3 partitions. This is documented, under Raises too.
• Tests:
  • The option reaches BART (its pole-detection log line appears with it and never without), for a slice, a volume and a trajectory.
  • The thin-volume edge, 2 partitions refused and 3 accepted, on the 24 central lines of a 64-point plane.
  • Cartesian calibration-grid maps match the full-grid fit: within 2e-2 of parallel at every object voxel, 1e-3 at the median.
  • Radial-plane small-grid maps match the whole-grid fit: within 2e-2 of parallel everywhere, 1e-3 at the median.
  • On a 32-cube kooshball, small-grid maps agree with the whole-grid fit to 0.99 at the median, and are closer to the true coils at the worst voxel (0.77 vs 0.74 from parallel).
  • The nlinv_pics trajectory test now runs on one thread. Its 32 spokes on a 16-point grid fit loosely enough that the order of a threaded sum moved the maps by 2e-2. It failed once at 4.4e-3 against a 1e-3 tolerance, and is bit-exact on one thread.
  • None of the synthetic coil sets I tried gave maps with a shared pole, so the correction itself was checked only on real data (below).

Measurements (CPU)

Old means the full-grid fit; new means the small-grid fit. |<old,new>| is the per-voxel |sum_c conj(old_c) new_c| (1 = parallel). "vs. true coils" is the worst voxel of the same quantity against the sensitivities that made the data.

Cartesian

Case	Fit time, old / new	|<old,new>| min / median	adj/rss min, old / new
64x64 phantom, 8 coils	0.28 s / 0.06 s	0.986 / 0.9991	0.907 / 0.946
16x64x64 slab, 8 coils	1.44 s / 0.21 s	0.984 / 0.995	0.681 / 0.697
7 T phantom slice, 220x220, 32 coils	13.7 s / 0.14 s	0.10 / 0.9995	0.105 / 0.769
same slice, phase_pole=0	29.0 s / 0.36 s	0.944 / 0.9995	0.938 / 0.982

Non-Cartesian (synthetic, radius=12)

Case	Fit time, old / new	|<old,new>| min / median	vs. true coils, old / new
radial plane 128, 8 coils, 201 spokes	2.16 s / 0.12 s	0.985 / 0.9996	0.988 / 0.990
radial plane 256, 8 coils, 403 spokes	7.64 s / 0.23 s	0.991 / 0.9997	0.994 / 0.993
same, phase_pole=0	25.2 s / 0.18 s	0.993 / 0.9996	0.995 / 0.997
kooshball 48³, 8 coils, 2000 spokes	13.2 s / 1.19 s	0.852 / 0.968	0.639 / 0.816
kooshball 96³, 8 coils, 9000 spokes	96.9 s / 2.20 s	0.607 / 0.897	0.531 / 0.929
same, phase_pole=0	171.8 s / 3.57 s	0.605 / 0.896	0.529 / 0.929

In three dimensions the two fits disagree because the old one is the worse of the two: it fitted the zeroed samples as data.

Phase-pole detection on the Cartesian path

Volume	phase_pole=None	phase_pole=0
2x64x64	0.09 s	refused
3x64x64	0.19 s	0.30 s
16x64x64	0.19 s	0.85 s
32x128x128	0.31 s	1.30 s

The default stays off, which keeps BART's behaviour, avoids the cost where the maps share no pole, and avoids the refusal on thin grids. Callers whose data show the spots ask for it; pulserver/pulserver#282 does, and needs this in a release.

Tests

These ran against the v0.0.15 library (no C changed):

• tests/test_apps.py, tests/test_docstrings.py (with numpydoc), tests/test_docs.py, tests/test_namespace.py, tests/test_options.py, tests/test_tools.py: 293 passed, 3 skipped (2 need CUDA, 1 in test_tools). The nlinv tests passed 3 runs in a row.
• scripts/lint.sh (ruff, format, codespell): clean.

I did not run the docs build.

🤖 Generated with Claude Code

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.

📢 Thoughts on this report? Let us know!