Merge pull request #133 from KrisThielemans/NeuroLF_Esser_Dataset

NeuroLF Esser dataset
SyneRBI · Oct 15, 2024 · 86c6565 · 86c6565
2 parents b3f8802 + 8771dc3
commit 86c6565
Show file tree

Hide file tree

Showing 6 changed files with 186 additions and 69 deletions.
diff --git a/SIRF_data_preparation/NeuroLF_Esser_Dataset/VOI_prep.py b/SIRF_data_preparation/NeuroLF_Esser_Dataset/VOI_prep.py
@@ -0,0 +1,76 @@
+# %% Description
+# We use  the Siemens mMR ACR VOIs from Pawel, as they match this phantom.
+# However, the phantom is scanned in a different position, and reg_aladin does not
+# figure it out, so we have to help it by doing manual repositioning first.
+
+# %%
+import os
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import sirf.Reg as Reg
+import sirf.STIR as STIR
+from SIRF_data_preparation.data_utilities import the_data_path, the_orgdata_path
+from SIRF_data_preparation.registration_utilities import STIR_to_nii, nii_to_STIR, register_it, resample_STIR
+
+# %% set paths filenames
+scanID = 'NeuroLF_Esser_Dataset'
+intermediate_data_path = the_orgdata_path(scanID, 'processing')
+orgdata_path = the_orgdata_path(scanID, scanID)
+os.makedirs(intermediate_data_path, exist_ok=True)
+data_path = the_data_path(scanID)
+os.makedirs(data_path, exist_ok=True)
+out_path = the_data_path(scanID, 'PETRIC')
+os.makedirs(out_path, exist_ok=True)
+mumap_filename = the_orgdata_path('Siemens_mMR_ACR', 'ACR_data_design/synth_mumap/acr-complete-umap.nii.gz')
+mMR_data_path = the_data_path('Siemens_mMR_ACR')
+
+# %% write current OSEM image as Nifti
+OSEM_image = STIR.ImageData(os.path.join(data_path, 'OSEM_image.hv'))
+OSEM_image_filename_nii = os.path.join(intermediate_data_path, 'OSEM_image.nii')
+OSEM_image_nii = STIR_to_nii(OSEM_image, OSEM_image_filename_nii)
+# %%
+mMR_OSEM_image = STIR.ImageData(os.path.join(mMR_data_path, 'OSEM_image.hv'))
+mMR_OSEM_image_filename_nii = os.path.join(intermediate_data_path, 'mMR_OSEM_image.nii')
+mMR_OSEM_image_nii = STIR_to_nii(mMR_OSEM_image, mMR_OSEM_image_filename_nii)
+# %% do a rotation over 180 degrees around the x-axis to get the phantom roughly in the correct place
+mMR_OSEM_image_nii_rot = Reg.ImageData(mMR_OSEM_image_nii)
+mMR_OSEM_image_nii_rot.fill(np.flip(mMR_OSEM_image_nii.as_array(), axis=(1, 2)))
+mMR_OSEM_image_nii_rot.write(os.path.join(intermediate_data_path, 'mMR_OSEM_image_rot.nii'))
+# %% do a rotation over 45 degrees around the z-axis to align the features
+alpha = 45 / 180. * np.pi
+m = Reg.AffineTransformation(
+    np.array(((np.cos(alpha), np.sin(alpha), 0, 0), (-np.sin(alpha), np.cos(alpha), 0, 0), (0, 0, 1, 0), (0, 0, 0, 1))))
+init_resampler = Reg.NiftyResampler()
+init_resampler.add_transformation(m)
+init_resampler.set_interpolation_type_to_nearest_neighbour()
+init_resampler.set_reference_image(mMR_OSEM_image_nii)
+init_resampler.set_floating_image(mMR_OSEM_image_nii)
+# %%
+mMR_OSEM_image_nii_rotrot = init_resampler.forward(mMR_OSEM_image_nii_rot)
+mMR_OSEM_image_nii_rotrot.write(os.path.join(intermediate_data_path, 'mMR_OSEM_image_rotrot.nii'))
+# %%
+plt.subplot(131)
+plt.imshow(mMR_OSEM_image_nii_rot.as_array()[:, :, 20])
+plt.subplot(132)
+plt.imshow(mMR_OSEM_image_nii_rotrot.as_array()[:, :, 20])
+plt.subplot(133)
+plt.imshow(OSEM_image_nii.as_array()[:, :, 20])
+# %% register
+(reg_mMR_nii, resampler, reg) = register_it(OSEM_image_nii, mMR_OSEM_image_nii_rotrot)
+
+# %% write
+reg_mMR_filename = os.path.join(intermediate_data_path, 'reg_mMR')
+reg_mMR = nii_to_STIR(reg_mMR_nii, reg_mMR_filename)
+
+# %%
+VOI_names = ['VOI_whole_object', 'VOI_background', 'VOI_cold_cylinder', 'VOI_hot_cylinder', 'VOI_rods']
+for VOI in VOI_names:
+    mMR_VOI = STIR.ImageData(os.path.join(mMR_data_path, 'PETRIC', VOI + '.hv'))
+    mMR_VOI_nii = STIR_to_nii(mMR_VOI, os.path.join(intermediate_data_path, 'mMR_' + VOI + '.nii'))
+    mMR_VOI_nii_rot = Reg.ImageData(mMR_VOI_nii)
+    mMR_VOI_nii_rot.fill(np.flip(mMR_VOI_nii.as_array(), axis=(1, 2)))
+    mMR_VOI_nii_rotrot = init_resampler.forward(mMR_VOI_nii_rot)
+    VOI_im = resample_STIR(resampler, mMR_VOI_nii_rotrot, os.path.join(intermediate_data_path, VOI))
+    VOI_im.write(os.path.join(out_path, VOI + '.hv'))
diff --git a/SIRF_data_preparation/NeuroLF_Esser_Dataset/prepare.py b/SIRF_data_preparation/NeuroLF_Esser_Dataset/prepare.py
@@ -0,0 +1,25 @@
+import os
+
+import sirf.STIR as STIR
+from SIRF_data_preparation.data_utilities import the_data_path, the_orgdata_path
+
+# %% set paths filenames
+scanID = 'NeuroLF_Esser_Dataset'
+intermediate_data_path = the_orgdata_path(scanID, 'processing')
+orgdata_path = the_orgdata_path(scanID, scanID)
+os.makedirs(intermediate_data_path, exist_ok=True)
+data_path = the_data_path(scanID)
+os.makedirs(data_path, exist_ok=True)
+# %%
+acquired_data = STIR.AcquisitionData(os.path.join(orgdata_path, 'neuroLF_prompts.hs'))
+norm_factors = STIR.AcquisitionData(os.path.join(orgdata_path, 'neuroLF_normalisation_factors.hs'))
+randoms = STIR.AcquisitionData(os.path.join(orgdata_path, 'neuroLF_random.hs'))
+scatter = STIR.AcquisitionData(os.path.join(orgdata_path, 'neuroLF_scatter.hs'))
+acfs = STIR.AcquisitionData(os.path.join(orgdata_path, 'neuroLF_acfs.hs'))
+
+# %%
+mult_factors = (acfs * norm_factors).power(-1.)
+additive_term = (scatter+randoms) * norm_factors
+acquired_data.write(os.path.join(data_path, 'prompts.hs'))
+mult_factors.write(os.path.join(data_path, 'mult_factors.hs'))
+additive_term.write(os.path.join(data_path, 'additive_term.hs'))
diff --git a/SIRF_data_preparation/create_Hoffman_VOIs.py b/SIRF_data_preparation/create_Hoffman_VOIs.py
@@ -29,10 +29,10 @@
 import scipy.ndimage as ndimage
 from docopt import docopt
 
-import sirf.Reg as Reg
 import sirf.STIR as STIR
 from SIRF_data_preparation.data_QC import plot_image
 from SIRF_data_preparation.data_utilities import the_data_path, the_orgdata_path
+from SIRF_data_preparation.registration_utilities import STIR_to_nii, STIR_to_nii_hv, register_it, resample_STIR
 
 # %%
 __version__ = "0.1.0"
@@ -80,70 +80,6 @@ def connected_component(arr: npt.NDArray[bool], order=0) -> npt.NDArray[bool]:
     return labels == idx
 
 
-# %% Functions to convert between STIR and Reg.ImageData (via writing to file)
-# converters directly via the objects are not all implemented, and we want them on file anyway
-# WARNING: There will be flips in nii_to_STIR if the STIR.ImageData contains PatientPosition
-# information AND it is not in HFS. This is because .nii cannot store this information
-# and STIR reads .nii as HFS.
-# WARNING: filename_prefix etc should not contain dots (.)
-def nii_to_STIR(image_nii: Reg.ImageData, filename_prefix: str) -> STIR.ImageData:
-    """Write Reg.ImageData as nii, read as STIR.ImageData and write that as well"""
-    image_nii.write(filename_prefix + ".nii")
-    image = STIR.ImageData(filename_prefix + ".nii")
-    image.write(filename_prefix + ".hv")
-    return image
-
-
-def STIR_to_nii(image: STIR.ImageData, filename_nii: str, filename_hv: str | None = None) -> Reg.ImageData:
-    """Write STIR.ImageData object as Nifti and read as Reg.ImageData"""
-    image.write_par(
-        filename_nii,
-        os.path.join(
-            STIR.get_STIR_examples_dir(),
-            "samples",
-            "stir_math_ITK_output_file_format.par",
-        ),
-    )
-    nii_image = Reg.ImageData(filename_nii)
-    if filename_hv is not None:
-        image.write(filename_hv)
-    return nii_image
-
-
-def STIR_to_nii_hv(image: STIR.ImageData, filename_prefix: str) -> Reg.ImageData:
-    """Call STIR_to_nii by appendix .nii and .hv"""
-    return STIR_to_nii(image, filename_prefix + ".nii", filename_prefix + ".hv")
-
-
-# %% Function to do rigid registration and return resampler
-def register_it(reference_image: Reg.ImageData,
-                floating_image: Reg.ImageData) -> typing.Tuple[Reg.ImageData, Reg.NiftyResampler]:
-    """
-    Use rigid registration of floating to reference image
-    Return both the registered image and the resampler.
-    Currently the resampler is set to use NN interpolation (such that it can be used for VOIs)
-    Note that `registered = resampler.forward(floating_image)` up to interpolation choices
-    """
-    reg = Reg.NiftyAladinSym()
-    reg.set_parameter("SetPerformRigid", "1")
-    reg.set_parameter("SetPerformAffine", "0")
-    reg.set_reference_image(reference_image)
-    reg.add_floating_image(floating_image)
-    reg.process()
-    resampler = Reg.NiftyResampler()
-    resampler.add_transformation(reg.get_deformation_field_forward())
-    resampler.set_interpolation_type_to_nearest_neighbour()
-    resampler.set_reference_image(reference_image)
-    resampler.set_floating_image(floating_image)
-    return (reg.get_output(0), resampler)
-
-
-# %% resample Reg.ImageData and write to file
-def resample_STIR(resampler: Reg.NiftyResampler, image_nii: Reg.ImageData, filename_prefix: str) -> STIR.ImageData:
-    res_nii = resampler.forward(image_nii)
-    return nii_to_STIR(res_nii, filename_prefix)
-
-
 # %%%%%%%%%%%%% Read Hoffman data and find VOIs
 
 
@@ -244,7 +180,7 @@ def create_Hoffman_VOIs(Hoffman: STIR.ImageData,) -> typing.Tuple[typing.List[ST
 orgGT = VOI_GM*5 + VOI_WM*1
 orgGT_nii = STIR_to_nii_hv(orgGT, os.path.join(intermediate_data_path, "orgGT"))
 
-regGT_nii, resampler = register_it(OSEM_image_nii, orgGT_nii)
+regGT_nii, resampler, _ = register_it(OSEM_image_nii, orgGT_nii)
 regGT = resample_STIR(resampler, orgGT_nii, os.path.join(intermediate_data_path, "regGT"))
 
 # %% check registration

diff --git a/SIRF_data_preparation/dataset_settings.py b/SIRF_data_preparation/dataset_settings.py
@@ -5,7 +5,8 @@
 
 DATA_SUBSETS = {
     'Siemens_mMR_NEMA_IQ': 7, 'Siemens_mMR_NEMA_IQ_lowcounts': 7, 'Siemens_mMR_ACR': 7, 'NeuroLF_Hoffman_Dataset': 16,
-    'Mediso_NEMA_IQ': 12, 'Siemens_Vision600_thorax': 5, 'GE_DMI3_Torso': 8, 'Siemens_Vision600_Hoffman': 5}
+    'Mediso_NEMA_IQ': 12, 'Siemens_Vision600_thorax': 5, 'GE_DMI3_Torso': 8, 'Siemens_Vision600_Hoffman': 5,
+    'NeuroLF_Esser_Dataset': 8}
 
 
 @dataclass

diff --git a/SIRF_data_preparation/registration_utilities.py b/SIRF_data_preparation/registration_utilities.py
@@ -0,0 +1,77 @@
+# Copyright 2024 University College London
+# Licence: Apache-2.0
+
+import os
+
+# %% imports
+import typing
+
+import sirf.Reg as Reg
+import sirf.STIR as STIR
+
+# %%
+__version__ = "0.1.0"
+
+
+# %% Functions to convert between STIR and Reg.ImageData (via writing to file)
+# converters directly via the objects are not all implemented, and we want them on file anyway
+# WARNING: There will be flips in nii_to_STIR if the STIR.ImageData contains PatientPosition
+# information AND it is not in HFS. This is because .nii cannot store this information
+# and STIR reads .nii as HFS.
+# WARNING: filename_prefix etc should not contain dots (.)
+def nii_to_STIR(image_nii: Reg.ImageData, filename_prefix: str) -> STIR.ImageData:
+    """Write Reg.ImageData as nii, read as STIR.ImageData and write that as well"""
+    image_nii.write(filename_prefix + ".nii")
+    image = STIR.ImageData(filename_prefix + ".nii")
+    image.write(filename_prefix + ".hv")
+    return image
+
+
+def STIR_to_nii(image: STIR.ImageData, filename_nii: str, filename_hv: str | None = None) -> Reg.ImageData:
+    """Write STIR.ImageData object as Nifti and read as Reg.ImageData"""
+    image.write_par(
+        filename_nii,
+        os.path.join(
+            STIR.get_STIR_examples_dir(),
+            "samples",
+            "stir_math_ITK_output_file_format.par",
+        ),
+    )
+    nii_image = Reg.ImageData(filename_nii)
+    if filename_hv is not None:
+        image.write(filename_hv)
+    return nii_image
+
+
+def STIR_to_nii_hv(image: STIR.ImageData, filename_prefix: str) -> Reg.ImageData:
+    """Call STIR_to_nii by appendix .nii and .hv"""
+    return STIR_to_nii(image, filename_prefix + ".nii", filename_prefix + ".hv")
+
+
+# %% Function to do rigid registration and return resampler
+def register_it(reference_image: Reg.ImageData,
+                floating_image: Reg.ImageData) -> typing.Tuple[Reg.ImageData, Reg.NiftyResampler, Reg.NiftyAladinSym]:
+    """
+    Use rigid registration of floating to reference image
+    Return both the registered image and the resampler.
+    Currently the resampler is set to use NN interpolation (such that it can be used for VOIs)
+    Note that `registered = resampler.forward(floating_image)` up to interpolation choices
+    """
+    reg = Reg.NiftyAladinSym()
+    reg.set_parameter("SetPerformRigid", "1")
+    reg.set_parameter("SetPerformAffine", "0")
+    reg.set_reference_image(reference_image)
+    reg.add_floating_image(floating_image)
+    reg.process()
+    resampler = Reg.NiftyResampler()
+    resampler.add_transformation(reg.get_deformation_field_forward())
+    resampler.set_interpolation_type_to_nearest_neighbour()
+    resampler.set_reference_image(reference_image)
+    resampler.set_floating_image(floating_image)
+    return (reg.get_output(0), resampler, reg)
+
+
+# %% resample Reg.ImageData and write to file
+def resample_STIR(resampler: Reg.NiftyResampler, image_nii: Reg.ImageData, filename_prefix: str) -> STIR.ImageData:
+    res_nii = resampler.forward(image_nii)
+    return nii_to_STIR(res_nii, filename_prefix)
diff --git a/petric.py b/petric.py
@@ -293,7 +293,7 @@ def get_image(fname):
     'Siemens_mMR_NEMA_IQ_lowcounts': {'transverse_slice': 72, 'coronal_slice': 109, 'sagittal_slice': 89},
     'Siemens_mMR_ACR': {'transverse_slice': 99}, 'NeuroLF_Hoffman_Dataset': {'transverse_slice': 72},
     'Mediso_NEMA_IQ': {'transverse_slice': 22, 'coronal_slice': 89, 'sagittal_slice': 66},
-    'Siemens_Vision600_thorax': {}, 'GE_DMI3_Torso': {}, 'Siemens_Vision600_Hoffman': {}}
+    'Siemens_Vision600_thorax': {}, 'GE_DMI3_Torso': {}, 'Siemens_Vision600_Hoffman': {}, 'NeuroLF_Esser_Dataset': {}}
 
 if SRCDIR.is_dir():
     # create list of existing data
@@ -314,7 +314,9 @@ def get_image(fname):
         (SRCDIR / "GE_DMI3_Torso", OUTDIR / "DMI3_Torso",
          [MetricsWithTimeout(outdir=OUTDIR / "DMI3_Torso", **DATA_SLICES['GE_DMI3_Torso'])]),
         (SRCDIR / "Siemens_Vision600_Hoffman", OUTDIR / "Vision600_Hoffman",
-         [MetricsWithTimeout(outdir=OUTDIR / "Vision600_Hoffman", **DATA_SLICES['Siemens_Vision600_Hoffman'])])]
+         [MetricsWithTimeout(outdir=OUTDIR / "Vision600_Hoffman", **DATA_SLICES['Siemens_Vision600_Hoffman'])]),
+        (SRCDIR / "NeuroLF_Esser_Dataset", OUTDIR / "NeuroLF_Esser",
+         [MetricsWithTimeout(outdir=OUTDIR / "NeuroLF_Esser", **DATA_SLICES['NeuroLF_Esser_Dataset'])])]
 else:
     log.warning("Source directory does not exist: %s", SRCDIR)
     data_dirs_metrics = [(None, None, [])] # type: ignore