Function r2_epg 

pub fn r2_epg(
    magnitude: &[f64],
    mask: &[u8],
    echo_times: &[f64],
    grid: &Grid,
    params: &R2EpgParams,
    b1_map: Option<&[f64]>,
) -> (Vec<f64>, Vec<f64>)

R2 mapping from multi-echo spin-echo magnitude via EPG dictionary matching.

For each masked voxel the (normalized) measured echo train is matched against a dictionary of EPG-simulated trains over a (T2, B1) grid; the R2 = 1/T2 of the best match is returned. Normalization removes the proton-density / receive gain, so only the decay shape is fit. This corrects the stimulated-echo bias that a mono-exponential fit ([r2star_arlo]) suffers on imperfectly-refocused (B1 < 1) data.

When a per-voxel b1_map is supplied (e.g. from a separate B1 acquisition), the B1 dimension of the search is fixed: each voxel’s B1 is snapped to the nearest params.b1_grid entry and only that column of the dictionary is searched (a T2-only fit). This keeps EPG’s stimulated-echo bias correction while removing the free B1 parameter that overfits noise at low SNR.

Arguments

• magnitude - MESE magnitude, flattened [v0_e0, v0_e1, ..., v1_e0, ...] (row-major (n_voxels, n_echoes), same layout as [r2star_arlo])
• mask - Binary brain mask [nx*ny*nz] (1 = process, 0 = skip)
• echo_times - Spin-echo times in seconds [n_echoes] (equi-spaced, ≥3)
• grid - Volume grid (dimensions and voxel sizes)
• params - Dictionary grids and assumed T1
• b1_map - Optional known refocusing-efficiency map [nx*ny*nz]; None fits B1 per voxel over the full b1_grid

Returns

(r2_map, b1_map) - R2 in Hz and the B1 (refocusing efficiency) used per voxel — fitted, or the snapped input when b1_map was given — both [nx*ny*nz].

Panics

Panics if echo_times.len() < 3, echo times are not equi-spaced, or the magnitude/mask/b1_map lengths are inconsistent with grid.