diff --git a/limri/norm/__init__.py b/limri/norm/__init__.py index bd0dccf..93078c0 100644 --- a/limri/norm/__init__.py +++ b/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 diff --git a/limri/norm/from_mean.py b/limri/norm/from_mean.py new file mode 100644 index 0000000..451bc74 --- /dev/null +++ b/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 diff --git a/limri/regtools.py b/limri/regtools.py index 4e94bad..765b553 100644 --- a/limri/regtools.py +++ b/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): diff --git a/limri/workflows/__init__.py b/limri/workflows/__init__.py index f8b67ea..4d18c7f 100644 --- a/limri/workflows/__init__.py +++ b/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) diff --git a/limri/workflows/normalization.py b/limri/workflows/normalization.py index bd0483f..0b1e90e 100644 --- a/limri/workflows/normalization.py +++ b/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) diff --git a/limri/workflows/registration.py b/limri/workflows/registration.py index d37ecd4..3b95712 100644 --- a/limri/workflows/registration.py +++ b/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)