Struct Problem 

struct Problem {

    dims: (usize, usize, usize),
    n: usize,
    dims0: (usize, usize, usize),
    pad: [usize; 3],
    vs: [f64; 3],
    m: Vec<u8>,
    inside: Support,
    drp: f64,
    drn: f64,
    d: Kernel,
    w_i: Vec<f64>,
    b1: Vec<f64>,
    xm: BoxVol,
    rm: BoxVol,
    wg_thr: [[(f64, f64); 3]; 2],
    cone: Support,
}

The arrays of one run, on the padded grid. What is fixed across the NNLS iterations is kept in the least memory that gives the same numbers (see Performance in the module docs).

Fields

dims: (usize, usize, usize)

n: usize

dims0: (usize, usize, usize)

input grid and padding

pad: [usize; 3]

vs: [f64; 3]

m: Vec<u8>

mask (0/1)

inside: Support

the mask’s voxels

drp: f64

D_r,± in ppm

drn: f64

d: Kernel

dipole kernel (dipole_kernel(..., 'kspace'))

w_i: Vec<f64>

step 1’s data weight W_I at the mask’s voxels: w_mag = W_I.*m there, 0 elsewhere, and w_phs = 2π·w_mag

b1: Vec<f64>

step 1’s right-hand side, [b+; b−]

xm: BoxVol

step 2’s gradient weights are recomputed where used from the masked fast QSM and smoothed R2′ (zero outside the mask, so kept on its bounding box) and their percentiles (lo, hi) per source and axis

rm: BoxVol

wg_thr: [[(f64, f64); 3]; 2]

cone: Support

step 2’s cone |D| < D2_thr

Implementations

impl Problem

fn new( inputs: &ChiSepIlsqrInputs<'_>, grid: &Grid, bdir: (f64, f64, f64), params: &ChiSepIlsqrParams, ) -> (Self, Workspace, Vec<f64>)

The problem, the FFT workspace, and R2′ (ppm, clamped) on the padded grid for the start.

fn pad_input<T: Copy + Default>(&self, a: &[T]) -> Vec<T>

padarray(a, pad, 0, 'both') of an input-grid volume

fn on_mask(&self, x: &[f64], clamp: bool) -> Vec<f64>

[x+; x−] of a padded [x+; x−] at the mask’s voxels, with a.*(a > 0) if clamp

fn padded_from_mask(&self, v: &[f64], x: &mut [f64])

The padded [x+.*m; x−.*m] of [x+; x−] at the mask’s voxels, into x

fn crop_from_mask(&self, v: &[f64]) -> (Vec<f64>, Vec<f64>)

(χ+, χ−) on the input grid from [x+; x−] at the mask’s voxels (χ− = −x−)

fn zero_ring(&self, a: &mut [f64])

Zero the padding ring (the toolbox crops each step’s output and pads the next step’s input)

fn load_w( &self, c: &mut [Complex64], f: impl Fn(usize, f64) -> f64 + Sync + Send, )

c ← f(l, w) as a real volume, with w = W_I at the mask’s voxels and 0 elsewhere (the value of w_mag there)

fn kmul(&self, c: &mut [Complex64])

c ← D.*c

fn dc_raw(&self, ws: &mut Workspace)

c ← ifftn(D.*fftn(c))

fn dc_c(&self, ws: &mut Workspace)

c ← Dc(c), c holding a real volume: real(ifftn(D.*fftn(c))), as a real volume

fn dc_out(&self, ws: &mut Workspace, out: &mut [f64])

out ← Dc(c), c holding a real volume

fn af_c(&self, ws: &mut Workspace)

c ← Af(c) = Dc(w_phs.*Dc(c)), c holding a real volume

fn a1(&self, ws: &mut Workspace, x: &[f64], out: &mut [f64])

Step 1’s normal operator on [x+; x−]

fn step1( &self, ws: &mut Workspace, start: Start<'_, f64>, tol: f64, maxit: usize, ) -> Vec<f64>

Step 1: pcg from start; returns the solution with the padding ring zeroed.

fn cone_ifft(&self, ws: &mut Workspace, y: &[Complex64], out: &mut [f64])

out ← real(ifftn(M_ic.*y)), y on the cone’s voxels

fn g_pm(&self, s: usize, c: usize, p: [usize; 3]) -> f64

The fgrad component c (voxel size vs) of the masked fast QSM and smoothed R2′ at p, combined into |g+| (s = 0) or |g−| (s = 1)

fn wg2(&self, s: usize, c: usize, p: [usize; 3], l: usize) -> f64

W_G² of source s, axis c at voxel p (l): conj(W).*W with W = W_G.*M

fn cone_fft_reg( &self, ws: &mut Workspace, s: usize, x: &[f64], out: &mut [Complex64], )

out ← M_ic.*fftn(gT(W_G².*fgrad(x, 1), 1)) on the cone’s voxels, for source s

fn step2(&self, ws: &mut Workspace, x: &mut [f64], params: &ChiSepIlsqrParams)

Step 2 on step 1’s [x+; x−], in place: each source minus its streaking component, masked.

fn fast_qsm(&self, ws: &mut Workspace, f: Vec<f64>, vs: [f64; 3]) -> Vec<f64>

The fast QSM of step 2 (see the module docs), on the padded grid.

fn tukey_lpf(&self, ws: &mut Workspace, r: &[f64]) -> Vec<f64>

R2′ (ppm) smoothed by a Tukey window in k-space: real(ifftn(fftshift(ifftshift(fftn(r)).*w)))