Merge pull request #3 from LoicDorval/ld_changes

Ld changes

limri/norm/__init__.py

@@ -13,3 +13,4 @@
 
 from .hist import hist_matching
 from .minmax import minmax_matching
+from .from_mean import from_mean_matching

limri/norm/from_mean.py

@@ -0,0 +1,37 @@
+# -*- 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.
+##########################################################################
+
+"""
+From mean send by Newcastle Team matching.
+"""
+
+# Imports
+import numpy as np
+
+
+def from_mean_matching(source, ref_val, concentration=2.):
+    """ Adjust the pixel values of a grayscale image.
+
+    Parameters
+    ----------
+    source: np.ndarray
+        the image to transform.
+    ref_val: int
+        reference value of phantom intensity for the corresponding site.
+    concentration: float, default 2.
+        the compartment concentration in milli mols / litre.
+
+    Returns
+    -------
+    matched: np.ndarray
+        the transformed source image.
+    """
+    # Normalize data
+    matched = (source * concentration) / ref_val
+    return matched

limri/regtools.py

@@ -237,7 +237,7 @@ def normalize2field():
     raise NotImplementedError
 
 
-def antsregister(template_file, lianat_file, hanat_file, outdir,
+def antsregister(template_file, li_file, lianat_file, hanat_file, outdir,
                  mask_file=None):
     """ Compute the deformation field with Ants from a T1w image to a template.
     """
@@ -248,6 +248,14 @@ def antsregister(template_file, lianat_file, hanat_file, outdir,
                           "function.")
 
     print_subtitle("Load data...")
+    li = ants.image_read(li_file)
+    print_result(f"li spacing: {li.spacing}")
+    print_result(f"li origin: {li.origin}")
+    print_result(f"li direction: {li.direction}")
+    filename = os.path.join(outdir, "li.png")
+    li.plot_ortho(
+        flat=True, xyz_lines=False, orient_labels=False,
+        title="li", filename=filename)
     lianat = ants.image_read(lianat_file)
     print_result(f"lianat spacing: {lianat.spacing}")
     print_result(f"lianat origin: {lianat.origin}")
@@ -279,13 +287,27 @@ def antsregister(template_file, lianat_file, hanat_file, outdir,
     template = ants.iMath_normalize(template)
 
     print_subtitle("Rigid: lianat -> hanat...")
-    li2h = ants.registration(
+    lianat2h = ants.registration(
         fixed=hanat, moving=lianat, type_of_transform="Rigid",
-        outprefix=os.path.join(outdir, "li2h"))
-    print_result(f"rigid transforms: {li2h['fwdtransforms']}")
+        outprefix=os.path.join(outdir, "lianat2h"))
+    print_result(f"rigid transforms: {lianat2h['fwdtransforms']}")
+    lianat2hanat = ants.apply_transforms(
+        fixed=hanat, moving=lianat, transformlist=lianat2h["fwdtransforms"],
+        interpolator="bSpline")
+    print(f"{lianat2hanat=}")
+    filename = os.path.join(outdir, "lianat2hanat.nii.gz")
+    lianat2hanat.to_filename(filename)
+    print_result(f"lianat2h T1: {filename}")
+    filename = os.path.join(outdir, "lianat2hanat.png")
+    lianat2hanat.plot_ortho(
+        hanat, flat=True, xyz_lines=False, orient_labels=False,
+        title="lianat2hanat", filename=filename, overlay_alpha=0.5)
+
+    print_subtitle("Rigid: li -> hanat...")
     li2hanat = ants.apply_transforms(
-        fixed=hanat, moving=lianat, transformlist=li2h["fwdtransforms"],
+        fixed=hanat, moving=li, transformlist=lianat2h["fwdtransforms"],
         interpolator="bSpline")
+    print(f"{li2hanat=}")
     filename = os.path.join(outdir, "li2hanat.nii.gz")
     li2hanat.to_filename(filename)
     print_result(f"li2h T1: {filename}")
@@ -316,35 +338,35 @@ def antsregister(template_file, lianat_file, hanat_file, outdir,
     jac = ants.create_jacobian_determinant_image(
         domain_image=hanat, tx=h2mni["fwdtransforms"][0])
     jac -= 1
-    h2mnianat = ants.apply_transforms(
+    hanat2mni = ants.apply_transforms(
         fixed=template, moving=hanat, transformlist=h2mni["fwdtransforms"],
         interpolator="bSpline")
     h2mnijac = ants.apply_transforms(
         fixed=template, moving=jac, transformlist=h2mni["fwdtransforms"],
         interpolator="bSpline")
-    li2mnianat = ants.apply_transforms(
+    lianat2mni = ants.apply_transforms(
         fixed=template, moving=lianat, interpolator="bSpline",
-        transformlist=h2mni["fwdtransforms"] + li2h["fwdtransforms"])
+        transformlist=h2mni["fwdtransforms"] + lianat2h["fwdtransforms"])
     filename = os.path.join(outdir, "hjac.nii.gz")
     jac.to_filename(filename)
     print_result(f"h jacobian: {filename}")
     filename = os.path.join(outdir, "h2mnijac.nii.gz")
     h2mnijac.to_filename(filename)
     print_result(f"h2mni jacobian: {filename}")
-    filename = os.path.join(outdir, "li2mnianat.nii.gz")
-    li2mnianat.to_filename(filename)
+    filename = os.path.join(outdir, "lianat2mni.nii.gz")
+    lianat2mni.to_filename(filename)
     print_result(f"li2mni T1: {filename}")
-    filename = os.path.join(outdir, "h2mnianat.nii.gz")
-    h2mnianat.to_filename(filename)
+    filename = os.path.join(outdir, "hanat2mni.nii.gz")
+    hanat2mni.to_filename(filename)
     print_result(f"h2mni T1: {filename}")
-    filename = os.path.join(outdir, "li2mnianat.png")
-    li2mnianat.plot_ortho(
+    filename = os.path.join(outdir, "lianat2mni.png")
+    lianat2mni.plot_ortho(
         template, flat=True, xyz_lines=False, orient_labels=False,
-        title="li2mnianat", filename=filename, overlay_alpha=0.5)
-    filename = os.path.join(outdir, "h2mnianat.png")
-    h2mnianat.plot_ortho(
+        title="lianat2mni", filename=filename, overlay_alpha=0.5)
+    filename = os.path.join(outdir, "hanat2mni.png")
+    hanat2mni.plot_ortho(
         template, flat=True, xyz_lines=False, orient_labels=False,
-        title="h2mnianat", filename=filename, overlay_alpha=0.5)
+        title="hanat2mni", filename=filename, overlay_alpha=0.5)
 
 
 def apply_transforms(fixed_file, moving_file, transformlist, filename):

limri/workflows/__init__.py

@@ -48,6 +48,6 @@ def li2mni_all(li_file, lianat_file, hanat_file, outdir, thr_factor=2,
     transformlist = [
         os.path.join(outdir, "h2mni1Warp.nii.gz"),
         os.path.join(outdir, "h2mni0GenericAffine.mat"),
-        os.path.join(outdir, "li2h0GenericAffine.mat"),
+        os.path.join(outdir, "lianat2h0GenericAffine.mat"),
         os.path.join(outdir, "li2lianat0GenericAffine.mat")]
     applytrf(ref_file, li_file, transformlist, shiftedli2mni_file)

limri/workflows/normalization.py

@@ -14,17 +14,19 @@
 # Imports
 import os
 import nibabel
-from limri.norm import hist_matching, minmax_matching
-from limri.color_utils import print_title, print_subtitle, print_result
+from limri.norm import hist_matching, minmax_matching, from_mean_matching
+from limri.color_utils import print_title, print_result
 
 # Global parameters
 NORM_MAP = {
     "hist": hist_matching,
-    "minmax": minmax_matching
+    "minmax": minmax_matching,
+    "mean": from_mean_matching
 }
 
 
-def li2mninorm(li2mni_file, li2mniref_file, mask_file, outdir, norm="hist"):
+def li2mninorm(li2mni_file, mask_file, outdir, ref_value=None,
+               li2mniref_file=None, norm="hist"):
     """ Normalize intensities using histogram matching.
 
     Parameters
@@ -33,26 +35,39 @@ def li2mninorm(li2mni_file, li2mniref_file, mask_file, outdir, norm="hist"):
         path to the Li image.
     li2mniref_file: str
         path to the reference Li image.
+    ref_value: int
+        reference value of phantom intensity for the corresponding site.
     mask_file: str
         the brain mask image.
     outdir: str
         path to the destination folder.
     norm: str, default 'hist'
-        the normalization method.
+        the normalization method, can be: hist, minmax, mean.
     """
     print_title("Load data...")
     li2mni = nibabel.load(li2mni_file)
     li2mni_arr = li2mni.get_fdata()
-    li2mniref = nibabel.load(li2mniref_file)
-    li2mniref_arr = li2mniref.get_fdata()
+    if li2mniref_file:
+        li2mniref = nibabel.load(li2mni_file)
+        li2mniref_arr = li2mniref.get_fdata()
     mask = nibabel.load(mask_file)
     mask_arr = mask.get_fdata()
+    if norm == "hist" or norm == "minmax" and li2mniref_file is None:
+        raise ValueError("We need the path to the reference Li image with"
+                         "li2mniref_file argument for this type of "
+                         "normalization method.")
 
     print_title("Normalization...")
     norm_fn = NORM_MAP.get(norm)
     if norm_fn is None:
         raise ValueError("Normalization method not defined.")
-    norm_arr = norm_fn(li2mni_arr, li2mniref_arr, mask_arr)
+    if norm != "mean":
+        norm_arr = norm_fn(li2mni_arr, li2mniref_arr, mask_arr)
+    elif norm == "mean":
+        if ref_value is None:
+            raise ValueError("We need the ref_value for this type of "
+                             "normalization method")
+        norm_arr = norm_fn(li2mni_arr, ref_value)
     norm = nibabel.Nifti1Image(norm_arr, li2mni.affine)
     norm_file = os.path.join(outdir, "li2mninorm.nii.gz")
     nibabel.save(norm, norm_file)

limri/workflows/registration.py

@@ -53,6 +53,13 @@ def li2mni(li_file, lianat_file, hanat_file, outdir, li2lianat=None):
     else:
         print_warning("hanat already reoriented")
     print_result(hanat_reo_file)
+    li_reo_file = os.path.join(outdir, "li.nii.gz")
+    if not os.path.isfile(li_reo_file):
+        gzfile(li_file, li_reo_file)
+        fslreorient2std(li_reo_file, li_reo_file, save_trf=True)
+    else:
+        print_warning("li already reoriented")
+    print_result(li_reo_file)
 
     print_title("Bias field correction...")
     lianat_bcorr_file = os.path.join(outdir, "lianat_restore.nii.gz")
@@ -76,7 +83,7 @@ def li2mni(li_file, lianat_file, hanat_file, outdir, li2lianat=None):
     print_result(hanat_bcorr_file)
 
     print_title("Coregistration & normalization...")
-    rigid_transforms = [os.path.join(outdir, "li2h0GenericAffine.mat")]
+    rigid_transforms = [os.path.join(outdir, "lianat2h0GenericAffine.mat")]
     deform_transforms = [
         os.path.join(outdir, "h2mni1Warp.nii.gz"),
         os.path.join(outdir, "h2mni0GenericAffine.mat")]
@@ -88,11 +95,12 @@ def li2mni(li_file, lianat_file, hanat_file, outdir, li2lianat=None):
     ref_file = os.path.join(os.path.dirname(limri.__file__), "resources",
                             "MNI152_T1_2mm.nii.gz")
     mask_file = os.path.join(os.path.dirname(limri.__file__), "resources",
-                            "MNI152_T1_2mm_brain.nii.gz")
+                             "MNI152_T1_2mm_brain.nii.gz")
     if not is_generated:
         antsregister(
             template_file=ref_file, lianat_file=lianat_bcorr_file,
-            hanat_file=hanat_bcorr_file, outdir=outdir, mask_file=mask_file)
+            li_file=li_file, hanat_file=hanat_bcorr_file, outdir=outdir,
+            mask_file=mask_file)
     else:
         print_warning("li2mni transformation already computed")
     print_result(deform_transforms + rigid_transforms)