Merge pull request #127 from KrisThielemans/Siemens_Vision600_Hoffman

- Siemens Vision600 Hoffman
- various improvements in QC and generation utilities

.pre-commit-config.yaml

@@ -1,3 +1,6 @@
+ci:
+  skip: [todo]
+
 default_language_version:
   python: python3
 repos:

SIRF_data_preparation/create_Hoffman_VOIs.py

@@ -8,6 +8,7 @@
 Options:
   -h, --help
   --dataset=<name>              dataset name (required)
+  --srcdir=<path>               pathname. Will default to current directory unless dataset is set
   -s, --skip_write_PETRIC_VOIs  do not write in data/<dataset>/PETRIC
 """
 
@@ -30,9 +31,8 @@
 
 import sirf.Reg as Reg
 import sirf.STIR as STIR
-from petric import OUTDIR, SRCDIR
 from SIRF_data_preparation.data_QC import plot_image
-from SIRF_data_preparation.data_utilities import the_orgdata_path
+from SIRF_data_preparation.data_utilities import the_data_path, the_orgdata_path
 
 # %%
 __version__ = "0.1.0"
@@ -47,24 +47,29 @@
     if scanID is None:
         print("Need to set the --dataset argument")
         exit(1)
-    if args["--skip_write_PETRIC_VOIs"] is not None:
+    if args["--skip_write_PETRIC_VOIs"]:
         write_PETRIC_VOIs = False
 else:
     # set it by hand, e.g.
     scanID = "NeuroLF_Hoffman_Dataset"
     write_PETRIC_VOIs = False
 
-# %% standard PETRIC directories
-if not all((SRCDIR.is_dir(), OUTDIR.is_dir())):
-    PETRICDIR = Path("~/devel/PETRIC").expanduser()
-    SRCDIR = PETRICDIR / "data"
-    OUTDIR = PETRICDIR / "output"
+srcdir = args['--srcdir']
 
+# %% standard PETRIC directories
+if srcdir is None:
+    srcdir = the_data_path(scanID)
+srcdir = Path(srcdir)
 downloaddir = Path(the_orgdata_path("downloads"))
 intermediate_data_path = the_orgdata_path(scanID, "processing")
 os.makedirs(downloaddir, exist_ok=True)
 os.makedirs(intermediate_data_path, exist_ok=True)
 
+print("srcdir:", srcdir)
+print("downloaddir:", downloaddir)
+print("processingdir:", intermediate_data_path)
+print("write_VOIs:", write_PETRIC_VOIs)
+
 
 # %% Function to find the n-th connected component (counting by size)
 def connected_component(arr: npt.NDArray[bool], order=0) -> npt.NDArray[bool]:
@@ -230,7 +235,6 @@ def create_Hoffman_VOIs(Hoffman: STIR.ImageData,) -> typing.Tuple[typing.List[ST
 # %%%%%%%%%%%%% Register VOIs to OSEM_image
 
 # %% Now register this to the reconstructed image
-srcdir = SRCDIR / scanID
 OSEM_image = STIR.ImageData(str(srcdir / "OSEM_image.hv"))
 OSEM_image_nii = STIR_to_nii(OSEM_image, os.path.join(intermediate_data_path, "OSEM_image.nii"))
 
@@ -263,6 +267,6 @@ def create_Hoffman_VOIs(Hoffman: STIR.ImageData,) -> typing.Tuple[typing.List[ST
     print(f"Writing registered VOIs to {datadir}")
     os.makedirs(datadir, exist_ok=True)
     for VOI, n in zip(regVOIs, VOInames):
-        VOI.write(str(datadir / n))
+        VOI.write(str(datadir / ("VOI_"+n)))
 # %%
 plt.show()

SIRF_data_preparation/data_QC.py

@@ -26,6 +26,7 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
+import numpy.typing as npt
 from docopt import docopt
 from scipy import ndimage
 
@@ -36,6 +37,15 @@
 STIR.AcquisitionData.set_storage_scheme('memory')
 
 
+def check_values_non_negative(arr: npt.NDArray[np.float32], desc: str):
+    min = np.min(arr)
+    max = np.max(arr)
+    if np.isnan(min) or min < 0:
+        raise ValueError(f"{desc}: minimum should be non-negative but is {min} (max={max})")
+    if not np.isfinite(max):
+        raise ValueError(f"{desc}: maximum should be finite")
+
+
 def plot_sinogram_profile(prompts, background, sumaxis=(0, 1), select=0, srcdir='./'):
     """
     Plot a profile through sirf.STIR.AcquisitionData
@@ -102,6 +112,13 @@ def plot_image_if_exists(prefix, **kwargs):
         return None
 
 
+def check_and_plot_image_if_exists(prefix, **kwargs):
+    im = plot_image_if_exists(prefix, **kwargs)
+    if im is not None:
+        check_values_non_negative(im.as_array(), prefix)
+    return im
+
+
 def VOI_mean(image, VOI):
     return float((image * VOI).sum() / VOI.sum())
 
@@ -123,6 +140,7 @@ def VOI_checks(allVOInames, OSEM_image=None, reference_image=None, srcdir='.', *
             continue
         VOI = STIR.ImageData(filename)
         VOI_arr = VOI.as_array()
+        check_values_non_negative(VOI_arr, VOIname)
         COM = np.rint(ndimage.center_of_mass(VOI_arr))
         num_voxels = VOI_arr.sum()
         print(f"VOI: {VOIname}: COM (in indices): {COM} voxels {num_voxels} = {num_voxels * np.prod(VOI.spacing)} mm^3")
@@ -180,10 +198,14 @@ def main(argv=None):
         mult_factors = STIR.AcquisitionData(os.path.join(srcdir, 'mult_factors.hs'))
         background = additive_term * mult_factors
         plot_sinogram_profile(acquired_data, background, srcdir=srcdir)
-
-    OSEM_image = plot_image_if_exists(os.path.join(srcdir, 'OSEM_image'), **slices)
-    plot_image_if_exists(os.path.join(srcdir, 'kappa'), **slices)
-    reference_image = plot_image_if_exists(os.path.join(srcdir, 'PETRIC/reference_image'), **slices)
+        check_values_non_negative(acquired_data.as_array(), "prompts")
+        check_values_non_negative(additive_term.as_array(), "additive_term")
+        check_values_non_negative(mult_factors.as_array(), "mult_factors")
+        check_values_non_negative(background.as_array(), "background")
+
+    OSEM_image = check_and_plot_image_if_exists(os.path.join(srcdir, 'OSEM_image'), **slices)
+    check_and_plot_image_if_exists(os.path.join(srcdir, 'kappa'), **slices)
+    reference_image = check_and_plot_image_if_exists(os.path.join(srcdir, 'PETRIC/reference_image'), **slices)
 
     VOIdir = os.path.join(srcdir, 'PETRIC')
     allVOInames = [os.path.basename(str(voi)[:-3]) for voi in Path(VOIdir).glob("VOI_*.hv")]

SIRF_data_preparation/dataset_settings.py

@@ -5,7 +5,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}
+    'Mediso_NEMA_IQ': 12, 'Siemens_Vision600_thorax': 5, 'GE_DMI3_Torso': 8, 'Siemens_Vision600_Hoffman': 5}
 
 
 @dataclass

SIRF_data_preparation/get_penalisation_factor.py

@@ -39,8 +39,8 @@
     if dataset is None or ref_dataset is None:
         print("Need to set the --dataset arguments")
         exit(1)
-    if args["--write_penalisation_factor"] is not None:
-        write_penalisation_factor = False
+    write_penalisation_factor = args["--write_penalisation_factor"]
+
 else:         # set it by hand, e.g.
     ref_dataset = "NeuroLF_Hoffman_Dataset"
     dataset = "Siemens_mMR_NEMA_IQ"

SIRF_data_preparation/run_BSREM.py

@@ -4,7 +4,10 @@
   run_BSREM.py <data_set> [--help | options]
 
 Arguments:
-  <data_set>  path to data files as well as prefix to use (e.g. Siemens_mMR_NEMA_EQ)
+  <data_set>     path to data files as well as prefix to use (e.g. Siemens_mMR_NEMA_EQ)
+
+Options:
+  --updates=<u>  number of updates to run [default: 15000]
 
 """
 # Copyright 2024 Rutherford Appleton Laboratory STFC
@@ -28,6 +31,7 @@
 # logging.basicConfig(level=logging.INFO)
 
 scanID = args['<data_set>']
+num_updates = int(args['--updates'])
 
 if not all((SRCDIR.is_dir(), OUTDIR.is_dir())):
     PETRICDIR = Path('~/devel/PETRIC').expanduser()
@@ -41,6 +45,9 @@
 settings = get_settings(scanID)
 
 data = get_data(srcdir=srcdir, outdir=outdir)
+print("Penalisation factor:", data.prior.get_penalisation_factor())
+print("num_subsets:", settings.num_subsets)
+print("num_updates:", num_updates)
 data_sub, acq_models, obj_funs = partitioner.data_partition(data.acquired_data, data.additive_term, data.mult_factors,
                                                             settings.num_subsets, mode="staggered",
                                                             initial_image=data.OSEM_image)
@@ -53,7 +60,7 @@
 algo = BSREM1(data_sub, obj_funs, initial=data.OSEM_image, initial_step_size=.3, relaxation_eta=.01,
               update_objective_interval=80)
 # %%
-algo.run(15000, callbacks=[MetricsWithTimeout(**settings.slices, outdir=outdir, seconds=3600 * 100)])
+algo.run(num_updates, callbacks=[MetricsWithTimeout(**settings.slices, interval=80, outdir=outdir, seconds=3600 * 100)])
 # %%
 fig = plt.figure()
 data_QC.plot_image(algo.get_output(), **settings.slices)

SIRF_data_preparation/run_OSEM.py

@@ -6,12 +6,13 @@
 Arguments:
   <data_set>  path to data files as well as prefix to use
 
+Options:
+  --updates=<u>  number of updates to run [default: 400]
 """
 # Copyright 2024 Rutherford Appleton Laboratory STFC
 # Copyright 2024 University College London
 # Licence: Apache-2.0
-__version__ = '0.1.0'
-
+__version__ = '0.2.0'
 from pathlib import Path
 
 import matplotlib.pyplot as plt
@@ -27,6 +28,7 @@
 # logging.basicConfig(level=logging.INFO)
 
 scanID = args['<data_set>']
+num_updates = int(args['--updates'])
 
 if not all((SRCDIR.is_dir(), OUTDIR.is_dir())):
     PETRICDIR = Path('~/devel/PETRIC').expanduser()
@@ -40,9 +42,11 @@
 settings = get_settings(scanID)
 
 data = get_data(srcdir=srcdir, outdir=outdir)
+print("num_subsets:", settings.num_subsets)
+print("num_updates:", num_updates)
 # %%
 algo = main_OSEM.Submission(data, settings.num_subsets, update_objective_interval=20)
-algo.run(400, callbacks=[MetricsWithTimeout(**settings.slices, seconds=5000, outdir=outdir)])
+algo.run(num_updates, callbacks=[MetricsWithTimeout(**settings.slices, seconds=5000, interval=20, outdir=outdir)])
 # %%
 fig = plt.figure()
 data_QC.plot_image(algo.get_output(), **settings.slices)

petric.py

@@ -178,17 +178,17 @@ def pass_index(metrics: np.ndarray, thresh: Iterable, window: int = 10) -> int:
         return np.where(res)[0][0]
 
 
-class MetricsWithTimeout(cil_callbacks.Callback):
+class MetricsWithTimeout(Callback):
     """Stops the algorithm after `seconds`"""
     def __init__(self, seconds=3600, outdir=OUTDIR, transverse_slice=None, coronal_slice=None, sagittal_slice=None,
                  **kwargs):
         super().__init__(**kwargs)
         self._seconds = seconds
         self.callbacks = [
             cil_callbacks.ProgressCallback(),
-            SaveIters(outdir=outdir),
+            SaveIters(outdir=outdir, **kwargs),
             (tb_cbk := StatsLog(logdir=outdir, transverse_slice=transverse_slice, coronal_slice=coronal_slice,
-                                sagittal_slice=sagittal_slice))]
+                                sagittal_slice=sagittal_slice, **kwargs))]
         self.tb = tb_cbk.tb # convenient access to the underlying SummaryWriter
         self.reset()
 
@@ -292,8 +292,8 @@ 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': {}}
+    'Mediso_NEMA_IQ': {'transverse_slice': 22, 'coronal_slice': 89, 'sagittal_slice': 66},
+    'Siemens_Vision600_thorax': {}, 'GE_DMI3_Torso': {}, 'Siemens_Vision600_Hoffman': {}}
 
 if SRCDIR.is_dir():
     # create list of existing data