fn artefact_solve(
n: [usize; 3],
x1: &[f64],
m: &[f64],
inside: &[usize],
gw: Vec<Vec<f64>>,
dk: &[f64],
crop_box1: &Box3,
tol: f64,
max_iter: usize,
tracer: &impl Tracer,
) -> (Vec<f64>, LsqrOutcome)Expand description
Step 6 of ilsqr: the streaking artefacts, by LSQR (tol) over
the k-space cone |D| < 0.1 against the weighted gradient (weights gw, full grid,
zero off the mask) of the initial solution x1. Returns the real part of F⁻¹ of the
solution (the artefacts divided by √N) and how LSQR ended.
A y = G ∇ (√N F⁻¹(M_ic y)), Aᴴ r = M_ic F(∇ᴴ(G r)) / √N. Every vector LSQR forms in the
domain is exactly zero off the cone, and every one in the range is exactly zero off the
mask (G = 0 there), so both are stored compactly: the unknown on the cone positions, the
residual on the mask positions of each axis. With norms split as for the dense vectors
(compact_parts) the iterates have the same values as a dense solve. The image-space data
lives in crop_box1 (the mask and its +e_a neighbours), so the transforms are pruned to it.