Function mcpc3ds_combine 

pub fn mcpc3ds_combine<P: AsRef<[f64]>, M: AsRef<[f64]>>(
    phases: &[Vec<P>],
    mags: &[Vec<M>],
    tes: &[f64],
    sigma: [f64; 3],
    echoes: [usize; 2],
    unwrap_method: UnwrapMethod,
    grid: &Grid,
) -> CoilCombinationResult

MCPC-3D-S multi-coil phase combination (Eckstein et al., MRM 2018), as in MriResearchTools.jl’s mcpc3ds for 5D (multi-echo, uncombined) input.

Each receive coil carries its own TE-independent phase offset, so uncombined channels cannot simply be summed. The algorithm estimates the offsets from the coil-summed Hermitian inner product (HIP) of two echoes and combines the channels coherently:

1. HIP = Σ_c |S_{1,c}| |S_{2,c}| exp(i(φ_{2,c} − φ_{1,c})) — the coil offsets cancel in the inter-echo phase difference, so the HIP phase is pure field evolution over ΔTE.
2. The HIP phase is unwrapped once (weighted by √|HIP|, on a robust mask of that weight).
3. Per coil, po_c = φ_{1,c} − (TE₁/ΔTE)·unwrapped_HIP, smoothed in the complex domain with the same masked box-filter Gaussian approximation as MriResearchTools (nan_box_smooth_3d_phase; sigma in voxels, the reference default is [10, 10, 5]).
4. Per echo, S_e = Σ_c |S_{e,c}|² exp(i(φ_{e,c} − po_c)); output phase arg(S_e) and magnitude sqrt(|S_e|).

With a single coil this reduces to phase_offset_removal (with the robust HIP mask).

Arguments

• phases - Wrapped phase, coil-major: phases[coil][echo] (each nx*ny*nz)
• mags - Magnitude, same layout as phases
• tes - Echo times (any consistent unit; only the ratio TE₁/ΔTE is used)
• sigma - Gaussian smoothing sigma in voxels [x, y, z] for the phase offsets
• echoes - Which two echoes form the HIP, default [0, 1]
• unwrap_method - How to unwrap the HIP phase (ROMEO or Laplacian)
• grid - Volume grid

Panics

If there are no coils, the coils do not all have the same number of echoes, magnitude and phase layouts differ, or echoes is out of range.