pub fn cosmos_weighted<F: AsRef<[f64]>, W: AsRef<[f64]>>(
local_fields: &[F],
weights: &[W],
bdirs: &[(f64, f64, f64)],
mask: &[u8],
grid: &Grid,
params: &CosmosParams,
) -> Vec<f64>Expand description
COSMOS dipole inversion with per-orientation spatial weighting.
This is the original formulation: each orientation’s residual is weighted by its signal
magnitude, so voxels where that orientation had poor SNR (or was corrupted, or fell
outside its own coverage) contribute less. Weighting is a spatial-domain operation and
the dipole convolution is a k-space one, so the two no longer diagonalize together and
the closed form of cosmos is unavailable. Instead the normal equations
( sum_t D_t W_t^2 D_t + lambda I ) chi = sum_t D_t W_t^2 f_tare solved with conjugate gradient. The operator is symmetric positive semi-definite
(D_t is a real, even k-space multiplier, hence self-adjoint), so CG applies directly.
Each iteration costs two FFTs per orientation, making this substantially more expensive
than cosmos — prefer the closed form unless the weighting earns its cost.
weights are multiplied by mask internally, so they need not be pre-masked. Typically
they are the per-orientation magnitude images, normalized to a comparable scale across
orientations (their relative size between orientations is what matters; a global factor
only rescales lambda).
§Panics
As cosmos, and additionally if the weight count or any weight map size is wrong.