Function resharp 

pub fn resharp(
    field: &[f64],
    mask: &[u8],
    grid: &Grid,
    params: &ResharpParams,
    progress: impl FnMut(usize, usize),
) -> (Vec<f64>, Vec<u8>)

RESHARP background field removal, as resharp.m (sunhongfu/QSM c75b138).

Like resharp.m, it zero-pads the field and mask by the kernel radius (l. 52-53), erodes the mask to the voxels whose kernel-weighted mean of the mask exceeds 0.999999 (l. 57-60), solves (H'H + λ)x = H'H y (H = M·ifft(DKER·fft(·)), DKER = δ − kernel) with MATLAB’s cgs at tolerance 1e-6 ([crate::solvers::cgs_solve_with_progress]: conjugate gradients squared with cgs’s stopping rules and returned iterate, l. 92), masks the result and crops the padding. resharp.m builds DKER from the kernel stored at index (r+1) and undoes the shift with circshift; this centres the kernel at index 0, the same operator.

Parity (QSMxT/QSM-CI resharp pair, run in astewartau/QSM.rs#404): the port previously ran conjugate gradients, which a numpy resharp.m showed to be exactly the 3.06 % (KKI) / 0.65 % (sim) gap: CG stops unconverged at 200 iterations where cgs converges in 199.

Arguments

• field - Unwrapped total field (nx * ny * nz)
• mask - Binary mask (nx * ny * nz), 1 = inside ROI
• grid - Volume dimensions and voxel sizes
• params - RESHARP algorithm parameters
• progress - Progress callback (current_iteration, max_iterations)

Returns

(local_field, eroded_mask)