Inati's method for csm is walsh / inefficient phase correction

[fzimmermann89]

We implement an algorithm called "inati" and refernce both the 2013 and 2014 abstract.
The 2013 abstract is mostly walsh method + using the a local phase refernce. (According to the abstract Walsh should use a reference coil for the phase, which we also dont do in walsh)

The two algorithms we implement also differ slightly in how they acutally perform the task - via filtering or sliding windows. But this should not matter.

Bart, tensorflow mri, ismrmrd-pythontools and most papers i find on scholar (as well as the LLMs) instead refer to his 2014 abstract as "inati's method".

This confusion stems from gadetron, which implements two algorithms, inati_iter (~2014 abstract) and inati (~2013).

I would check which variant (inati or walsh) uses less memory and is faster.
Both should provide similar results if one comments out the "mrpro inati" phase correction.
Then we should keep the faster/memory efficient variant and add a switch for the 2013 inati phase correction.

And add a proper inati method that uses the inati 2014 approach to get smother non-singular phases of images and sensitivity maps.

[fzimmermann89]

as a starting point:

"""Inati method for coil sensitivity map calculation."""

import torch
from einops import einsum

from mrpro.data.SpatialDimension import SpatialDimension
from mrpro.utils.filters import uniform_filter


def inati(
   coil_img: torch.Tensor,
   smoothing_width: SpatialDimension[int] | int,
   n_iterations: int = 5,
) -> torch.Tensor:
   """Calculate a coil sensitivity map (csm) using the Inati method [INA2013]_ [INA2014]_.

   This function computes CSMs using the iterative method described in [INA2014]_.
   It enforces smoothness in both magnitude and phase and aligns the phase to a
   global reference to avoid phase singularities.

   The input should be a tensor with dimensions (..., coils, z, y, x). The output
   will have the same dimensions.

   .. [INA2013] Inati S, Hansen M, Kellman P (2013) A solution to the phase problem in adaptvie coil combination.
      in Proceedings of the 21st Annual Meeting of ISMRM, Salt Lake City, USA, 2672.

   .. [INA2014] Inati S, Hansen M, Kellman P (2014) A Fast Optimal Method for Coil Sensitivity Estimation and
      Adaptive Coil Combination for Complex Images. in Proceedings of Joint Annual Meeting ISMRM-ESMRMB, Milan,
      Italy, 7115.

   Parameters
   ----------
   coil_img
       images for each coil element
   smoothing_width
       Size of the smoothing kernel
   n_iterations
       Number of iterations to run

   Returns
   -------
   csm
       coil sensitivity map
   """
   eps = 1e-12

   if isinstance(smoothing_width, int):
       smoothing_width = SpatialDimension(
           z=smoothing_width if coil_img.shape[-3] > 1 else 1,
           y=smoothing_width,
           x=smoothing_width,
       )

   if any(ks % 2 != 1 for ks in [smoothing_width.z, smoothing_width.y, smoothing_width.x]):
       raise ValueError('kernel_size must be odd')

   d_sum = torch.sum(coil_img, dim=(-3, -2, -1))
   d_sum /= d_sum.norm(dim=-1, keepdim=True) + eps

   combined_img = einsum(d_sum.conj(), coil_img, '... c, ... c z y x -> ... z y x')

   for _ in range(n_iterations):
       csm = coil_img * combined_img.conj().unsqueeze(-4)
       csm = uniform_filter(csm, width=smoothing_width.zyx, dim=(-3, -2, -1))
       csm /= csm.norm(dim=-4, keepdim=True) + eps

       combined_img = einsum(coil_img, csm.conj(), '... c z y x, ... c z y x -> ... z y x')

       d_sum = (csm * combined_img.unsqueeze(-4)).sum(dim=(-3, -2, -1))
       d_sum /= d_sum.norm(dim=-1, keepdim=True) + eps

       phase = einsum(d_sum.conj(), csm, '... c, ... c z y x -> ... z y x').angle()
       combined_img *= torch.exp(1j * phase)
       csm *= torch.exp(-1j * phase).unsqueeze(-4)

   return csm