limri: update pipeline.

doc/requirements.txt

@@ -4,4 +4,5 @@ ipython
 numpy
 nibabel
 antspyx
+dipy
 sphinx==3.5.4

limri/denoising.py

@@ -0,0 +1,52 @@
+# -*- coding: utf-8 -*-
+##########################################################################
+# NSAp - Copyright (C) CEA, 2023
+# Distributed under the terms of the CeCILL-B license, as published by
+# the CEA-CNRS-INRIA. Refer to the LICENSE file or to
+# http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+# for details.
+##########################################################################
+
+"""
+Denoising tools.
+"""
+
+# Import
+from dipy.denoise.nlmeans import nlmeans
+from dipy.denoise.noise_estimate import estimate_sigma
+
+
+def nlm_denoising(arr, n_coils=0):
+    """ Non-local means for denoising 3D and 4D images, using blockwise
+    averaging approach.
+
+    Parameters
+    ----------
+    arr: 3D or 4D ndarray
+        The array to be denoised
+    n_coils: int, default 0
+        Number of coils of the receiver array. Use N = 1 in case of a SENSE
+        reconstruction (Philips scanners) or the number of coils for a GRAPPA
+        reconstruction (Siemens and GE). Use 0 to disable the correction
+        factor, as for example if the noise is Gaussian distributed.
+
+    Returns
+    -------
+    denoised_arr: ndarray
+        the denoised ``arr`` which has the same shape as ``arr``.
+
+    References
+    ----------
+    .. [Coupe08] P. Coupe, P. Yger, S. Prima, P. Hellier, C. Kervrann, C.
+                 Barillot, An Optimized Blockwise Non Local Means Denoising
+                 Filter for 3D Magnetic Resonance Images, IEEE Transactions on
+                 Medical Imaging, 27(4):425-441, 2008
+    .. [Coupe11] Pierrick Coupe, Jose Manjon, Montserrat Robles, Louis Collins.
+                Adaptive Multiresolution Non-Local Means Filter for 3D MR Image
+                Denoising IET Image Processing, Institution of Engineering and
+                Technology, 2011
+    """
+    sigma = estimate_sigma(arr, N=n_coils)
+    denoised_arr = nlmeans(arr, sigma=sigma, patch_radius=1, block_radius=2,
+                           rician=True)
+    return denoised_arr

limri/info.py

@@ -67,7 +67,8 @@
 REQUIRES = [
     "numpy",
     "nibabel",
-    "antspyx"
+    "antspyx",
+    "dipy"
 ]
 SCRIPTS = [
     "limri/scripts/limri"

limri/norm/__init__.py

@@ -0,0 +1,15 @@
+# -*- coding: utf-8 -*-
+##########################################################################
+# NSAp - Copyright (C) CEA, 2023
+# Distributed under the terms of the CeCILL-B license, as published by
+# the CEA-CNRS-INRIA. Refer to the LICENSE file or to
+# http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+# for details.
+##########################################################################
+
+"""
+Normalization/calibration tools.
+"""
+
+from .hist import hist_matching
+from .minmax import minmax_matching

limri/norm/hist.py

@@ -0,0 +1,129 @@
+# -*- coding: utf-8 -*-
+##########################################################################
+# NSAp - Copyright (C) CEA, 2023
+# Distributed under the terms of the CeCILL-B license, as published by
+# the CEA-CNRS-INRIA. Refer to the LICENSE file or to
+# http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+# for details.
+##########################################################################
+
+"""
+Histogram matching.
+"""
+
+# Imports
+import numpy as np
+
+
+def hist_matching(source, template, mask, plot=False):
+    """ Adjust the pixel values of a grayscale image such that its histogram
+    matches that of a target image.
+
+    Parameters
+    ----------
+    source: np.ndarray
+        the image to transform: the histogram is computed over the flattened
+        array.
+    template: np.ndarray
+        the template image: same dimensions as the source image.
+    mask: np.ndarray
+        the mask image: same dimensions as the source image.
+    plot: bool, default False
+        plot the matched histograms.
+
+    Returns
+    -------
+    matched: np.ndarray
+        the transformed source image.
+    """
+    # Compute the source and template histograms
+    mask_indices = np.where(mask == 1)
+    h1, c1 = np.histogram(source[mask_indices], bins=65536)
+    h1t, c1t = np.histogram(template[mask_indices], bins=65536)
+    vec_size = len(h1)
+    vec_size_f = float(vec_size)
+
+    # Compute the normalized cumulated density functions
+    cdf1 = h1.cumsum().astype(np.float64) / np.sum(h1)
+    cdf2 = h1t.cumsum().astype(np.float64) / np.sum(h1t)
+
+    # Compute the mapping
+    M = np.zeros((vec_size, ))
+    for idx in range(vec_size):
+        ind = np.argmin(np.abs(cdf1[idx] - cdf2))
+        M[idx] = ind
+    M /= vec_size_f
+    B = np.asarray(range(vec_size)) / vec_size_f
+
+    # Interpolate the data by fitting a piecewise linear interpolant
+    data1x = (source - c1.min()) / c1.max()
+    data1x[np.where(data1x < 0)] = 0
+    data1x[np.where(data1x > 1)] = 1
+    indices = np.where((data1x > 0) & (data1x <= 1))
+    xnew = data1x[indices]
+    ynew = np.interp(xnew, B, M)
+
+    # Plot if verbose
+    if plot:
+        import matplotlib.pyplot as plt
+        fig, ax = plt.subplots()
+        ax.plot(B, M, 'k-')
+        ax.plot(xnew[::10000], ynew[::10000], 'ro')
+        plt.show()
+
+    # Reshape the result
+    matched = np.zeros(source.shape, dtype=source.dtype)
+    matched[indices] = ynew * (c1t.max() - c1t.min()) + c1t.min()
+
+    return matched
+
+
+def _hist_matching(source, template, mask, plot=False):
+    """ Adjust the pixel values of a grayscale image such that its histogram
+    matches that of a target image.
+
+    Parameters
+    ----------
+    source: np.ndarray
+        the image to transform: the histogram is computed over the flattened
+        array.
+    template: np.ndarray
+        the template image: same dimensions as the source image.
+    mask: np.ndarray
+        the mask image: same dimensions as the source image.
+    plot: bool, default False
+        plot the matched histograms.
+
+    Returns
+    -------
+    matched: np.ndarray
+        the transformed source image.
+    """
+    # Image information
+    shape = source.shape
+    dtype = source.dtype
+    mask_indices = np.where(mask == 1)
+    source = source[mask_indices]
+    template = template[mask_indices]
+
+    # Get the set of unique pixel values and their corresponding indices and
+    # counts
+    s_values, bin_idx, s_counts = np.unique(source, return_inverse=True,
+                                            return_counts=True)
+    t_values, t_counts = np.unique(template, return_counts=True)
+
+    # Take the cumsum of the counts and normalize by the number of pixels to
+    # get the empirical cumulative distribution functions for the source and
+    # template images (maps pixel value --> quantile)
+    s_quantiles = np.cumsum(s_counts).astype(np.float64)
+    s_quantiles /= s_quantiles[-1]
+    t_quantiles = np.cumsum(t_counts).astype(np.float64)
+    t_quantiles /= t_quantiles[-1]
+
+    # Interpolate linearly to find the pixel values in the template image
+    # that correspond most closely to the quantiles in the source image
+    interp_t_values = np.interp(s_quantiles, t_quantiles, t_values)
+    matched = np.zeros(shape, dtype=dtype)
+    matched[mask_indices] = interp_t_values[bin_idx]
+
+    return matched

limri/norm/minmax.py

@@ -0,0 +1,53 @@
+# -*- coding: utf-8 -*-
+##########################################################################
+# NSAp - Copyright (C) CEA, 2023
+# Distributed under the terms of the CeCILL-B license, as published by
+# the CEA-CNRS-INRIA. Refer to the LICENSE file or to
+# http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+# for details.
+##########################################################################
+
+"""
+Histogram matching.
+"""
+
+# Imports
+import numpy as np
+from sklearn import mixture
+
+
+def minmax_matching(source, template, mask, concentration=2.):
+    """ Adjust the pixel values of a grayscale image such that its dynamic
+    matches that of a target image.
+
+    Parameters
+    ----------
+    source: np.ndarray
+        the image to transform.
+    template: np.ndarray
+        the template image: same dimensions as the source image. In this case
+        the template is a phantom with 1 compartment.
+    mask: np.ndarray
+        the mask image: same dimensions as the source image.
+    concentration: float, default 2.
+        the compartment concentration in milli mols / litre.
+
+    Returns
+    -------
+    matched: np.ndarray
+        the transformed source image.
+    """
+    # Segment the compartment
+    clf = mixture.GaussianMixture(n_components=2, covariance_type="full")
+    clf.fit(template.reshape(-1, 1))
+    m1, m2 = clf.means_
+    thr = (m1 + m2) / 2.
+    template[template < thr] = 0.
+
+    # Get the reference value
+    ref_val = np.mean(template[template > 0])
+
+    # Normalize data
+    matched = (source * concentration) / ref_val
+
+    return matched

limri/regtools.py

@@ -237,7 +237,8 @@ def normalize2field():
     raise NotImplementedError
 
 
-def antsregister(template_file, lianat_file, hanat_file, outdir):
+def antsregister(template_file, lianat_file, hanat_file, outdir,
+                 mask_file=None):
     """ Compute the deformation field with Ants from a T1w image to a template.
     """
     try:
@@ -294,9 +295,23 @@ def antsregister(template_file, lianat_file, hanat_file, outdir):
         title="li2hanat", filename=filename, overlay_alpha=0.5)
 
     print_subtitle("Rigid + Affine + deformation field: hanat -> template...")
-    h2mni = ants.registration(
-        fixed=template, moving=hanat, type_of_transform="SyNRA",
-        outprefix=os.path.join(outdir, "h2mni"))
+    if mask_file is None:
+        h2mni = ants.registration(
+            fixed=template, moving=hanat, type_of_transform="SyNRA",
+            outprefix=os.path.join(outdir, "h2mni"))
+    else:
+        h2mni = ants.registration(
+            fixed=template, moving=hanat, type_of_transform="Affine",
+            outprefix=os.path.join(outdir, "_h2mni"))
+        mask = ants.image_read(mask_file)
+        print_result(f"mask spacing: {mask.spacing}")
+        print_result(f"mask origin: {mask.origin}")
+        print_result(f"mask direction: {mask.direction}")
+        h2mni = ants.registration(
+            fixed=template, moving=hanat, type_of_transform="SyNOnly",
+            mask=mask, initial_transform=h2mni["fwdtransforms"][0],
+            outprefix=os.path.join(outdir, "h2mni"))
+
     print_result(f"deform transforms: {h2mni['fwdtransforms']}")
     jac = ants.create_jacobian_determinant_image(
         domain_image=hanat, tx=h2mni["fwdtransforms"][0])

limri/scripts/limri

@@ -17,5 +17,6 @@ fire.Fire({
     "li2mni-all": wf.li2mni_all,
     "li2mni": wf.li2mni,
     "applytrf": wf.applytrf,
-    "li2mnieyes": wf.li2mnieyes
+    "li2mnieyes": wf.li2mnieyes,
+    "li2mninorm": wf.li2mninorm
 })

limri/workflows/__init__.py

@@ -15,6 +15,7 @@
 import limri
 from .registration import li2mni, applytrf
 from .maskeyes import li2mnieyes
+from .normalization import li2mninorm
 
 
 def li2mni_all(li_file, lianat_file, hanat_file, outdir, thr_factor=2,

limri/workflows/maskeyes.py

@@ -22,6 +22,7 @@
 from collections import Counter
 import matplotlib.pyplot as plt
 import limri
+from limri.denoising import nlm_denoising
 from limri.regtools import save_translation
 from limri.color_utils import print_title, print_subtitle, print_result
 
@@ -50,6 +51,13 @@ def li2mnieyes(li2mni_file, outdir, thr_factor=2, bins=300):
     ref_im = nibabel.load(ref_file)
     ref_arr = ref_im.get_fdata()
 
+    print_title("Denoising...")
+    arr = nlm_denoising(arr, n_coils=0)
+    denoise_im = nibabel.Nifti1Image(arr, im.affine)
+    li2mnidenoised_file = os.path.join(outdir, "li2mnidenoised.nii.gz")
+    nibabel.save(denoise_im, li2mnidenoised_file)
+    print_result(li2mnidenoised_file)
+
     print_title("Last peak extraction: GMM...")
     data = arr[arr > 0]
     data.shape += (1, )