Merge pull request #135 from KrisThielemans/GE_D690_NEMA_IQ

GE_D690_NEMA_IQ

SIRF_data_preparation/GE_D690_NEMA_IQ/README.md

@@ -0,0 +1,11 @@
+# Preparation steps for NEMA IQ data scanned on GE Discovery 690
+## Reference solution
+```sh
+python -m SIRF_data_preparation.run_BSREM GE_D690_NEMA_IQ  --updates=500 --initial_step_size=1 --interval=10 --outreldir=BSREM1
+python -m SIRF_data_preparation.run_BSREM GE_D690_NEMA_IQ --updates=7520 --initial_step_size=.8 --relaxation_eta=.01  --num_subsets=6 --initial_image=output/GE_D690_NEMA_IQ/BSREM1/iter_final.hv --outreldir=BSREM2
+python -m SIRF_data_preparation.run_BSREM GE_D690_NEMA_IQ --updates=1000 --initial_step_size=.2 --relaxation_eta=.01  --num_subsets=1 --initial_image=output/GE_D690_NEMA_IQ/BSREM2/iter_7520.hv --outreldir=BSREM3 --interval=2
+```
+## VOIs
+```sh
+python -m SIRF_data_preparation.create_NEMA_IQ_VOIs --dataset=GE_D690_NEMA_IQ --central_VOI=False --angle_smallest_sphere=30 --spheres="(2,3,5)"
+```

SIRF_data_preparation/create_Hoffman_VOIs.py

@@ -35,7 +35,7 @@
 from SIRF_data_preparation.registration_utilities import STIR_to_nii, STIR_to_nii_hv, register_it, resample_STIR
 
 # %%
-__version__ = "0.1.0"
+__version__ = "0.1.1"
 
 write_PETRIC_VOIs = True
 if "ipykernel" not in sys.argv[0]: # clunky way to be able to set variables from within jupyter/VScode without docopt
@@ -44,6 +44,7 @@
     # logging.basicConfig(level=logging.INFO)
 
     scanID = args["--dataset"]
+    srcdir = args['--srcdir']
     if scanID is None:
         print("Need to set the --dataset argument")
         exit(1)
@@ -52,10 +53,9 @@
 else:
     # set it by hand, e.g.
     scanID = "NeuroLF_Hoffman_Dataset"
+    srcdir = None
     write_PETRIC_VOIs = False
 
-srcdir = args['--srcdir']
-
 # %% standard PETRIC directories
 if srcdir is None:
     srcdir = the_data_path(scanID)

SIRF_data_preparation/create_NEMA_IQ_VOIs.py

@@ -0,0 +1,165 @@
+#!/usr/bin/env python
+"""Create VOIs for a scan of the NEMA IQ phantom
+
+Usage:
+  create_NEMA_IQ_VOIs.py [--help | options]
+
+
+Options:
+  -h, --help
+  --dataset=<name>              dataset name (required)
+  --srcdir=<path>               pathname. Will default to current directory unless dataset is set
+  --central_VOI=<bool>          Use the normal background VOI (in the centre of the phantom) or a
+                                VOI "below" the largest sphere [default: True]
+  --angle_smallest_sphere=<a>   angle (in degrees) of smallest sphere [default: 210]
+  --spheres=<s>                 string to be parsed as tuple of sphere numbers (1 is smallest) [default: (1,3,5)]
+"""
+
+# Copyright 2024 University College London
+# Licence: Apache-2.0
+
+# %%
+# %load_ext autoreload
+# %autoreload 2
+# %% imports
+import ast
+import os
+import sys
+import typing
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+import numpy.typing as npt
+import scipy.ndimage as ndimage
+from docopt import docopt
+
+import sirf.contrib.NEMA.generate_nema_rois as generate_nema_rois
+import sirf.STIR as STIR
+from SIRF_data_preparation.data_QC import VOI_checks, plot_image
+from SIRF_data_preparation.data_utilities import the_data_path, the_orgdata_path
+from SIRF_data_preparation.dataset_settings import get_settings
+
+# %%
+__version__ = "0.1.0"
+
+if "ipykernel" not in sys.argv[0]: # clunky way to be able to set variables from within jupyter/VScode without docopt
+    args = docopt(__doc__, argv=None, version=__version__)
+
+    # logging.basicConfig(level=logging.INFO)
+
+    scanID = args["--dataset"]
+    srcdir = args['--srcdir']
+    if scanID is None:
+        print("Need to set the --dataset argument")
+        exit(1)
+    angle_smallest_sphere = args["--angle_smallest_sphere"]
+    central_VOI = args["--central_VOI"]
+    spheres = ast.literal_eval(args["--spheres"])
+else:
+    # set it by hand, e.g.
+    scanID = "GE_D690_NEMA_IQ"
+    srcdir = None
+    central_VOI = False
+    spheres = (2, 3, 5)
+    angle_smallest_sphere = 30
+
+# %% standard PETRIC directories
+if srcdir is None:
+    srcdir = the_data_path(scanID)
+srcdir = Path(srcdir)
+intermediate_data_path = Path(the_orgdata_path(scanID, "processing"))
+os.makedirs(intermediate_data_path, exist_ok=True)
+settings = get_settings(scanID)
+slices = settings.slices
+
+print("srcdir:", srcdir)
+print("processingdir:", intermediate_data_path)
+print("angle_smallest_sphere:", angle_smallest_sphere)
+print("central_VOI:", central_VOI)
+print("spheres:", spheres)
+
+
+# %% Function to find the n-th connected component (counting by size)
+def connected_component(arr: npt.NDArray[bool], order=0) -> npt.NDArray[bool]:
+    """Return connected component of a given order in an array. order=0 returns the largest."""
+    labels, num_labels = ndimage.label(arr)
+    num_elems = [(labels == i).sum() for i in range(num_labels)]
+    idx = np.argsort(num_elems)[-1 - order]
+    return labels == idx
+
+
+def read_unregistered_VOIs(intermediate_datadir: Path) -> typing.List[STIR.ImageData]:
+    orgVOIs = []
+    for i in range(1, 8):
+        orgVOIs.append(STIR.ImageData(str(intermediate_datadir / f"unregistered_sphere{i}.nii")))
+    return orgVOIs
+
+
+def read_orgVOIs(intermediate_datadir: Path) -> typing.List[STIR.ImageData]:
+    orgVOIs = []
+    for i in range(1, 8):
+        orgVOIs.append(STIR.ImageData(str(intermediate_datadir / f"S{i}.nii")))
+    return orgVOIs
+
+
+def create_whole_object_mask(reference_image: STIR.ImageData) -> STIR.ImageData:
+    whole_object_mask = reference_image.clone()
+    ref_arr = reference_image.as_array()
+    whole_object_mask.fill(ndimage.binary_erosion(ref_arr > np.percentile(ref_arr, 99) / 20, iterations=2))
+    # plot_image(whole_object_mask, **slices)
+    return whole_object_mask
+
+
+def create_background_VOI(orgVOIs: typing.List[STIR.ImageData], central_VOI: bool) -> STIR.ImageData:
+    if central_VOI:
+        return orgVOIs[6]
+    # This is the NEMA with "lung" so need to move background somewhere else
+    # will use the largest sphere VOI and shift it down
+    background_mask = orgVOIs[5].clone()
+    background_arr = background_mask.as_array()
+    # COM = ndimage.center_of_mass(background_arr)
+    z_voxel_size = background_mask.voxel_sizes()[0]
+    # 37 mm diameter
+    z_shift = int(np.floor(37 * 1.2 / z_voxel_size))
+    background_arr = np.concatenate((
+        background_arr[0:z_shift, :, :] * 0,
+        background_arr[0:-z_shift, :, :],
+    ))
+    background_mask.fill(background_arr)
+    return background_mask
+
+
+# %%
+reference_image = STIR.ImageData(str(srcdir / 'OSEM_image.hv'))
+plot_image(reference_image, **slices)
+# %%
+thresholded_image = reference_image.clone()
+thresholded_arr = thresholded_image.as_array()
+thresholded_arr[thresholded_arr < np.max(thresholded_arr) / 5] = 0
+thresholded_image.fill(thresholded_arr)
+# %%
+# currently needs trailing slash
+generate_nema_rois.data_output_path = str(intermediate_data_path) + "/"
+generate_nema_rois.generate_nema_rois(thresholded_image, angle_smallest=angle_smallest_sphere)
+# %%
+unregVOIs = read_unregistered_VOIs(intermediate_data_path)
+orgVOIs = read_orgVOIs(intermediate_data_path)
+# %%
+whole_object_mask = create_whole_object_mask(reference_image)
+background_mask = create_background_VOI(orgVOIs, central_VOI)
+sphere_VOIs = [orgVOIs[i - 1] for i in spheres]
+# %%
+VOIs = [whole_object_mask, background_mask] + sphere_VOIs
+VOInames = ["VOI_whole_object", "VOI_background"] + [f"VOI_sphere{s}" for s in spheres]
+
+# %% write PETRIC VOIs
+datadir = Path(srcdir) / "PETRIC"
+print(f"Writing VOIs to {datadir}")
+os.makedirs(datadir, exist_ok=True)
+for VOI, n in zip(VOIs, VOInames):
+    VOI.write(str(datadir / n))
+# %%
+VOI_checks(VOInames, OSEM_image=reference_image, reference_image=None, srcdir=datadir)
+# %%
+plt.show()

SIRF_data_preparation/data_QC.py

@@ -167,7 +167,7 @@ def VOI_checks(allVOInames, OSEM_image=None, reference_image=None, srcdir='.', *
 
     if OSEM_image is not None:
         plt.figure()
-        plot_image(OSEM_image, alpha=allVOIs, save_name="OSEM_image_and_VOIs", **kwargs)
+        plot_image(OSEM_image, alpha=allVOIs, save_name=prefix + "OSEM_image_and_VOIs", **kwargs)
 
     # unformatted print of VOI values for now
     print(allVOInames)

SIRF_data_preparation/dataset_settings.py

@@ -6,7 +6,7 @@
 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,
-    'NeuroLF_Esser_Dataset': 8, 'Siemens_Vision600_ZrNEMAIQ': 5}
+    'NeuroLF_Esser_Dataset': 8, 'Siemens_Vision600_ZrNEMAIQ': 5, 'GE_D690_NEMA_IQ': 16}
 
 
 @dataclass

SIRF_data_preparation/plot_BSREM_metrics.py

@@ -26,29 +26,37 @@
 # scanID = 'NeuroLF_Hoffman_Dataset'
 # scanID = 'Siemens_mMR_NEMA_IQ'
 # scanID = 'Mediso_NEMA_IQ'
-scanID = 'Siemens_Vision600_Hoffman'
+# scanID = 'Siemens_Vision600_Hoffman'
 # scanID = 'NeuroLF_Esser_Dataset'
 # scanID = 'Siemens_Vision600_ZrNEMAIQ'
+scanID = 'GE_D690_NEMA_IQ'
 srcdir = the_data_path(scanID)
 outdir = OUTDIR / scanID
 OSEMdir = outdir / 'OSEM'
 datadir = outdir / 'BSREM'
 # we will check for images obtained after restarting BSREM (with new settings)
 datadir1 = outdir / 'BSREM_cont'
-
+# datadir = Path('/opt/runner/logs/Casper/BSREM/') / scanID
 OSEM_image = STIR.ImageData(str(datadir / 'iter_0000.hv'))
 
 settings = get_settings(scanID)
 slices = settings.slices
 
-cmax = numpy.percentile(OSEM_image.as_array(), 99.995)
+cmax = numpy.percentile(OSEM_image.as_array(), 99.9)
+# %%
+data = get_data(srcdir=srcdir, outdir=None, read_sinos=False)
 # %%
-image = data_QC.plot_image_if_exists(str(datadir / 'iter_final'), **slices, vmax=cmax)
+if data.reference_image is not None:
+    reference_image = data.reference_image
+elif datadir1.is_dir():
+    reference_image = STIR.ImageData(str(datadir1 / 'iter_final.hv'))
+else:
+    reference_image = STIR.ImageData(str(datadir / 'iter_final.hv'))
 # image2=STIR.ImageData(datadir+'iter_14000.hv')
 # %%
 image2 = OSEM_image
-diff = image2 - image
-print("relative l1-norm diff final-OSEM:", diff.abs().max() / image.max())
+diff = image2 - reference_image
+print("relative l1-norm diff final-OSEM:", diff.abs().max() / reference_image.max())
 data_QC.plot_image(diff, **slices, vmin=-cmax / 100, vmax=cmax / 100)
 # %%
 objs = read_objectives(datadir)
@@ -69,14 +77,6 @@
 fig.savefig(outdir / f'{scanID}_BSREM_objectives_last.png')
 
 # %%
-data = get_data(srcdir=srcdir, outdir=None, read_sinos=False)
-# %%
-if data.reference_image is not None:
-    reference_image = data.reference_image
-elif datadir1.is_dir():
-    reference_image = STIR.ImageData(str(datadir1 / 'iter_final.hv'))
-else:
-    reference_image = STIR.ImageData(str(datadir / 'iter_final.hv'))
 qm = QualityMetrics(reference_image, data.whole_object_mask, data.background_mask, tb_summary_writer=None,
                     voi_mask_dict=data.voi_masks)
 # %% get update ("iteration") numbers from objective functions

SIRF_data_preparation/run_BSREM.py

@@ -1,4 +1,6 @@
-"""Run BSREM for PETRIC
+#!/usr/bin/env python
+"""
+Run BSREM for PETRIC
 
 Usage:
   run_BSREM.py <data_set> [--help | options]
@@ -9,17 +11,18 @@
 Options:
   --updates=<u>               number of updates to run [default: 15000]
   --initial_image=<filename>  optional initial image, normally the OSEM_image from get_data.
-                              If specified,
-                              output will be in BSREM_cont.
+
   --num_subsets=<n>           number of subsets. If not specified, will use dataset_settings.get_settings.
   --initial_step_size=<s>     start stepsize [default: .3]
   --relaxation_eta=<r>        relaxation factor per epoch [default: .01]
   --interval=<i>              interval to save [default: 80]
+  --outreldir=<relpath>       optional relative path to override
+                              (defaults to 'BSREM' or 'BSREM_cont' if initial_image is set)
 """
 # Copyright 2024 Rutherford Appleton Laboratory STFC
 # Copyright 2024 University College London
 # Licence: Apache-2.0
-__version__ = '0.3.0'
+__version__ = '0.4.0'
 
 from pathlib import Path
 
@@ -44,6 +47,7 @@
 initial_step_size = float(args['--initial_step_size'])
 relaxation_eta = float(args['--relaxation_eta'])
 interval = int(args['--interval'])
+outreldir = args['--outreldir']
 
 if not all((SRCDIR.is_dir(), OUTDIR.is_dir())):
     PETRICDIR = Path('~/devel/PETRIC').expanduser()
@@ -61,11 +65,11 @@
 if initial_image is None:
     initial_image_name = "OSEM"
     initial_image = data.OSEM_image
-    outdir = outdir / "BSREM"
+    outdir = outdir / ("BSREM" if outreldir is None else outreldir)
 else:
     initial_image_name = initial_image
     initial_image = STIR.ImageData(initial_image)
-    outdir = outdir / "BSREM_cont"
+    outdir = outdir / ("BSREM_cont" if outreldir is None else outreldir)
 
 print("Penalisation factor:", data.prior.get_penalisation_factor())
 print("num_subsets:", num_subsets)

petric.py

@@ -292,9 +292,9 @@ def get_image(fname):
     'Siemens_mMR_NEMA_IQ': {'transverse_slice': 72, 'coronal_slice': 109, 'sagittal_slice': 89},
     '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': {},
-    'NeuroLF_Esser_Dataset': {}, 'Siemens_Vision600_ZrNEMAIQ': {'transverse_slice': 60}}
+    'Mediso_NEMA_IQ': {'transverse_slice': 22, 'coronal_slice': 89, 'sagittal_slice': 66},
+    'Siemens_Vision600_thorax': {}, 'GE_DMI3_Torso': {}, 'Siemens_Vision600_Hoffman': {}, 'NeuroLF_Esser_Dataset': {},
+    'Siemens_Vision600_ZrNEMAIQ': {'transverse_slice': 60}, 'GE_D690_NEMA_IQ': {'transverse_slice': 23}}
 
 if SRCDIR.is_dir():
     # create list of existing data
@@ -319,7 +319,9 @@ def get_image(fname):
         (SRCDIR / "NeuroLF_Esser_Dataset", OUTDIR / "NeuroLF_Esser",
          [MetricsWithTimeout(outdir=OUTDIR / "NeuroLF_Esser", **DATA_SLICES['NeuroLF_Esser_Dataset'])]),
         (SRCDIR / "Siemens_Vision600_ZrNEMAIQ", OUTDIR / "Vision600_ZrNEMAIQ",
-         [MetricsWithTimeout(outdir=OUTDIR / "Vision600_ZrNEMAIQ", **DATA_SLICES['Siemens_Vision600_ZrNEMAIQ'])])]
+         [MetricsWithTimeout(outdir=OUTDIR / "Vision600_ZrNEMAIQ", **DATA_SLICES['Siemens_Vision600_ZrNEMAIQ'])]),
+        (SRCDIR / "GE_D690_NEMA_IQ", OUTDIR / "D690_NEMA",
+         [MetricsWithTimeout(outdir=OUTDIR / "D690_NEMA", **DATA_SLICES['GE_D690_NEMA_IQ'])])]
 else:
     log.warning("Source directory does not exist: %s", SRCDIR)
     data_dirs_metrics = [(None, None, [])] # type: ignore