Module cosmos 

COSMOS: calculation of susceptibility through multiple orientation sampling

Single-orientation dipole inversion is ill-posed because the dipole kernel vanishes on a pair of cones at the magic angle, so those spatial frequencies carry no information and have to be regularized back in. COSMOS removes the ill-posedness at acquisition time instead: the cones sit at a fixed angle to B0, so re-imaging the same object at a different orientation moves them somewhere else. Combining N orientations, every frequency is sampled well by at least one of them, and the inversion becomes a well-conditioned least-squares problem needing no regularization at all.

Minimizing the sum of squared residuals over orientations,

    chi = argmin_x  sum_t || F^-1 D_t F x  -  f_t ||^2

diagonalizes in k-space (every operator involved is a k-space multiplication), giving a closed form with no iteration:

    chi = Re . IFFT(  sum_t D_t . FFT(f_t)  /  ( sum_t |D_t|^2 + lambda )  )

cosmos is that expression. cosmos_weighted is the original formulation, which weights each orientation’s residual by its signal magnitude; spatial weights do not commute with the k-space operators, so that variant loses the closed form and solves the normal equations with conjugate gradient instead.

Orientations

Each orientation contributes its own B0 direction expressed in the voxel frame of the common (registered) grid — the object is resampled into one shared frame, so what varies between orientations is where B0 points relative to it. Derive each one with crate::geometry::b0_direction_from_affine from that acquisition’s affine, which factors the voxel sizes out before taking the direction and so stays correct for anisotropic voxels.

How many orientations, and what limits accuracy

k-space coverage saturates fast. Measured on a 164x205x205 1 mm grid, the number of k-points where every orientation’s kernel is negligible goes 16952 (one orientation, the magic-angle cone) to 16 (two orientations 60 degrees apart, where the two cones still intersect along curves) to 0 (three). Even the real in-vivo orientations of the reference dataset, which span only 0 to 25 degrees from B0 because that is as far as a head turns inside a coil, leave nothing uncovered once there are twelve of them — a 17-degree rotation already shifts the cone clear of where it started.

So beyond about three orientations, extra scans buy noise averaging rather than coverage. That still matters a great deal: with field noise at 5% of the field RMS, reconstruction error on the same phantom runs roughly 205 / 0.09 / 0.03 (relative NRMSE) for one, two and four orientations. The single-orientation case is catastrophic precisely because the cone is uncovered and there is no redundancy to average against — which is the whole argument for acquiring more than one.

What is not fixed by adding orientations:

• Anisotropy. COSMOS models susceptibility as a scalar. White matter is measurably anisotropic, and that mismatch does not shrink with more orientations; it is the motivation for susceptibility tensor imaging.
• Registration. Misalignment between orientations propagates straight into the result.
• Rotations about a single axis. These leave residual streaking that L2 regularization suppresses (Schweser & Zivadinov 2024) — one reason lambda is exposed.

COSMOS is widely treated as a susceptibility ground truth; those three are what that status actually rests on, not the orientation count.

When the orientations are too close

The closed form divides by sum_t |D_t|^2, so its noise gain at a frequency is 1 / sqrt(sum_t |D_t|^2). A single orientation’s cone of |D| < 0.1 (CONE_THRESHOLD) is where that gain exceeds 10, and with lambda = 0 nothing damps it. Orientations only help where they move that cone, so the question for a set is how much of it the set still leaves ill-conditioned: uncovered_cone, the share of a single orientation’s cone where sum_t |D_t|^2 is still below CONE_THRESHOLD^2. One orientation leaves all of it. Three pointing the same way leave 57%: repeats average the noise down by sqrt(3), which narrows the cone, but they cover none of it. The share falls further as the directions tilt apart.

The largest angle between directions is not the measure: two orientations always leave the curves where their cones intersect, and three that tilt about one axis leave more than three spread around B0. Measured with uncovered_cone:

Directions	Uncovered
Three runs of one head position, re-registered (B0 within 1.1°, 2.6°, 3.3° of each other)	55%
Two orientations 5° / 10° / 15° / 20° apart	64% / 35% / 22% / 16%
Three of the twelve in-vivo orientations of Bilgic et al. (2015): median, worst	0.4%, 28%
Three at 30° to each other, all twelve of Bilgic et al.	0%

The first is what three runs of a single-orientation protocol give: registration recovers the head’s small movements between them, and COSMOS from those directions is a single-orientation inversion with no regularization, its cone filled with amplified noise. MAX_UNCOVERED_CONE is the share above which a set is that rather than a multi-orientation acquisition: it refuses the first two rows below about 13°, and passes designed acquisitions and most real triples.

Units and gauge

Local (tissue) field in ppm in, susceptibility in ppm out, matching the rest of the crate. The dipole kernel is zero at DC, so the recovered chi has zero mean by construction; referencing to a tissue region is a separate pipeline stage (crate::pipeline::apply_reference).

References

Liu T, Spincemaille P, de Rochefort L, Kressler B, Wang Y. “Calculation of susceptibility through multiple orientation sampling (COSMOS): a method for conditioning the inverse problem from measured magnetic field map to susceptibility source image in MRI.” Magnetic Resonance in Medicine. 2009;61(1):196-204. https://doi.org/10.1002/mrm.21828

Bilgic B, Fan AP, Polimeni JR, et al. “Fast quantitative susceptibility mapping with L1-regularization and automatic parameter selection” / “Rapid multi-orientation quantitative susceptibility mapping.” NeuroImage. 2015;111:622-630. https://doi.org/10.1016/j.neuroimage.2015.02.036

Schweser F, Zivadinov R. “Chaos and COSMOS - considerations on QSM methods with multiple and single orientations and effects from local anisotropy.” Magnetic Resonance Imaging. 2024;109:135-146. https://doi.org/10.1016/j.mri.2024.04.001

Structs

CosmosParams

COSMOS parameters.

Constants

CONE_THRESHOLD

The dipole kernel magnitude below which a frequency is ill-conditioned: where one orientation’s noise gain, 1 / |D|, exceeds 10. The cone uncovered_cone measures.

MAX_UNCOVERED_CONE

The most of a single orientation’s cone a direction set may leave ill-conditioned and still be treated as multi-orientation (module docs, When the orientations are too close).

Functions

check_inputs 🔒

Validate the shared preconditions of both entry points and return (n_orientations, n_voxels).

cosmos

COSMOS dipole inversion, closed form.

cosmos_weighted

COSMOS dipole inversion with per-orientation spatial weighting.

uncovered_cone

The share of a single orientation’s ill-conditioned cone (|D| < CONE_THRESHOLD) that the directions bdirs still leave ill-conditioned: where sum_t |D_t|^2 < CONE_THRESHOLD^2. 1 for one direction; 0 when the set covers the cone. Repeats of one direction only average noise: three leave 57%. Directions need not be unit length; a zero one is skipped.