Function fit_voxel 

fn fit_voxel(
    te: &[f64],
    yc: &[(f64, f64)],
    logyc: &[(f64, f64)],
    b0: f64,
    upper: f64,
    params: &DecomposeParams,
    ws: &mut StageSolvers,
) -> Fit

The authors’ three-stage alternating fit for one voxel (DECOMPOSE_analysis.m):

chi_0 = [0.05 -0.05]; R0 = 25; C0 = [0.3 0.3 0.4];
for it = 1:N_inner
    C0    = lsqcurvefit(objective1 (C),  C0,    TE, complex_to_real(ydata),      [0 0 0], [inf inf inf]);
    R0    = lsqcurvefit(objective2 (R0), R0,    TE, complex_to_real(log(ydata)), 0,       inf);
    chi_0 = lsqcurvefit(objective3 (χ),  chi_0, TE, complex_to_real(log(ydata)), [0 -0.5], [0.5 0]);
end

with lsqcurvefit‘s default options, which [lsq_trust_region_reflective] reproduces. The χ bound upper is DecomposeParams::chi_bound, resolved (the authors’ 0.5 by default). lsqcurvefit’s residual is fun(x, TE) − ydata, packed by complex_to_real as [Re; −Im].