new data

.gitignore

@@ -16,4 +16,4 @@ __pycache__/
 *.v
 *.csv
 # Tensorboard files
-events.*
+events.*

SIRF_data_preparation/GE_DMI3_Torso/README.md

@@ -0,0 +1,25 @@
+# Torso phantom scanned on a GE Discovery MI (3 ring)
+
+This data is from a scan performed by GE of the
+Data Spectrum [Antropomorphic Torso Phantom](https://www.spect.com/pdf/anthropomorphic-torso-phantom.pdf)
+with additional lesion inserts of diameter 10mm.
+
+Activity ratio between "liver" and "torso" is approximately 2. Lungs
+are empty (air). Contrast of the lesions is unfortunately unknown.
+
+"Corrections" were obtained with GE Duetto and converted to STIR Interfile
+using internal scripts.
+
+Due to data size issues and count levels, the supplied data is **reduced** from the original scan
+by keeping only the "direct" sinograms (segment 0, corresponding to min/max ring differences -1/+1), e.g. on the command line
+```
+SSRB prompts.hs prompts_f1b1.hs 1 1 0 0 1
+```
+
+In addition, as current PETRIC examples do not cope with non-TOF multfactors for TOF data, the
+multfactors were expanded to TOF (by duplication).
+
+## VOIs
+
+VOIs were determined by Kris Thielemans on reconstructed images and approximately placed
+where some of the lesions are visible. They are not necessarily at the correct location.

SIRF_data_preparation/GE_DMI3_Torso/VOI_prep.py

@@ -0,0 +1,83 @@
+# %% file to prepare VOIs for the Torso phantom, using coordinates derived manually
+import os
+from pathlib import Path
+
+import numpy as np
+from scipy.ndimage import binary_fill_holes
+
+import sirf.STIR as STIR
+from SIRF_data_preparation.data_utilities import the_data_path
+
+# %%
+scanID = 'GE_DMI3_Torso'
+data_path = Path(the_data_path(scanID))
+output_path = Path(the_data_path(scanID, 'PETRIC'))
+os.makedirs(output_path, exist_ok=True)
+# %%
+image_filename = str(data_path / 'OSEM_image.hv')
+image = STIR.ImageData(image_filename)
+# %% find whole phantom VOI by thresholding/filling
+threshold = image.max() * .05
+im_arr = image.as_array() >= threshold
+# we want to use binary_fill_holes to fill the spine and lung. However,
+# that function does not fill holes connected to the boundary, so let's
+# fill the boundary planes first
+im_arr[0, :, :] = True
+im_arr[-1, :, :] = True
+im_arr = binary_fill_holes(im_arr)
+# now exclude boundary
+im_arr[0, :, :] = False
+im_arr[-1, :, :] = False
+VOI = image.allocate(0)
+VOI.fill(im_arr.astype(np.float32))
+VOI.show()
+VOI.write(str(output_path / 'VOI_whole_object.hv'))
+# %% find centre image as KT got the coordinates from Amide
+# note that geo.get_index_to_physical_point_matrix() is still unstable, so we won't use it.
+# Warning: the following will likely fail for future versions of SIRF
+zcentre = (image.shape[0] - 1) / 2 * image.voxel_sizes()[0]
+# %% background VOI (in approx middle of torso)
+shape = STIR.EllipticCylinder()
+shape.set_length(23)
+shape.set_radii((120, 120))
+shape.set_origin((zcentre + 0, 0, 0))
+VOI.fill(0)
+VOI.add_shape(shape, scale=1)
+VOI.write(str(output_path / 'VOI_background.hv'))
+# %% spine (approximately top half of it)
+shape = STIR.EllipticCylinder()
+shape.set_length(100)
+shape.set_radii((22, 22))
+shape.set_origin((zcentre + 42.6, 86.9, 3))
+VOI.fill(0)
+VOI.add_shape(shape, scale=1)
+VOI.write(str(output_path / 'VOI_spine.hv'))
+VOI.show()
+# %% lung lesion
+shape = STIR.Ellipsoid()
+shape.set_radius_z(5)
+shape.set_radius_y(5)
+shape.set_radius_x(5)
+shape.set_origin((zcentre + 25, -42.2, -105.3))
+VOI.fill(0)
+VOI.add_shape(shape, scale=1)
+VOI.write(str(output_path / 'VOI_lung_lesion.hv'))
+# %%
+# %% liver lesion
+shape = STIR.Ellipsoid()
+shape.set_radius_z(5)
+shape.set_radius_y(5)
+shape.set_radius_x(5)
+shape.set_origin((zcentre - 38.2, -38.9, -49.7))
+VOI.fill(0)
+VOI.add_shape(shape, scale=1)
+VOI.write(str(output_path / 'VOI_liver_lesion.hv'))
+# %% lung lesion
+shape = STIR.Ellipsoid()
+shape.set_radius_z(5)
+shape.set_radius_y(5)
+shape.set_radius_x(5)
+shape.set_origin((zcentre + 47.6, 12.5, 151.5))
+VOI.fill(0)
+VOI.add_shape(shape, scale=1)
+VOI.write(str(output_path / 'VOI_chest_lesion.hv'))

SIRF_data_preparation/README.md

@@ -27,7 +27,7 @@ for the Siemens mMR NEMA IQ data (on Zenodo):
 - `BSREM_*.py`: functions with specific settings for a particular data-set
 
 # Steps to follow to prepare data
-If starting from Siemens list-mode data and letting SIRF take care of scatter etc, check for instance [Siemens_mMR_ACR/README.md](steps for Siemens mMR ACR). If pre-prepared data are given, check that naming of all files is correct. KT normally puts all data
+If starting from Siemens mMR list-mode data and letting SIRF take care of scatter etc, check for instance [steps for Siemens mMR ACR](Siemens_mMR_ACR/README.md). If pre-prepared data are given, check that naming of all files is correct. KT normally puts all data
 in `~/devel/PETRIC/data/<datasetname>` with `datasetname` following convention of `scanner_phantom-name` as others (instructions below and indeed some scripts might assume this location). Change working directory to where data sits and add PETRIC to your python-path, e.g.
 ```
 PYTHONPATH=~/devel/PETRIC:$PYTHONPATH`

SIRF_data_preparation/Siemens_mMR_NEMA_IQ/prepare_mMR_NEMA_IQ_data.py

@@ -2,7 +2,7 @@
 import logging
 import os
 
-from ..data_utilities import prepare_challenge_Siemens_data, the_data_path, the_orgdata_path
+from ..data_utilities import prepare_challenge_Siemens_data, the_data_path
 
 this_directory = os.path.dirname(__file__)
 repo_directory = os.path.dirname(this_directory)

SIRF_data_preparation/Siemens_mMR_NEMA_IQ_lowcounts/README.md

@@ -5,4 +5,4 @@ This relies on the "normal" count level processsing to be done first, specifical
 Steps to follow:
 1. `prepare_mMR_NEMA_IQ_data.py`
 2. copy VOIs
-3. further steps as [normal](../README.md#steps-to-follow-to-prepare-data)
+3. further steps as [normal](../README.md#steps-to-follow-to-prepare-data)

SIRF_data_preparation/data_QC.py

@@ -30,8 +30,8 @@
 from scipy import ndimage
 
 import sirf.STIR as STIR
-from SIRF_data_preparation.dataset_settings import get_settings
 from SIRF_data_preparation.data_utilities import the_data_path
+from SIRF_data_preparation.dataset_settings import get_settings
 
 STIR.AcquisitionData.set_storage_scheme('memory')
 

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

petric.py

@@ -43,7 +43,7 @@ class Callback(cil_callbacks.Callback):
     CIL Callback but with `self.skip_iteration` checking `min(self.interval, algo.update_objective_interval)`.
     TODO: backport this class to CIL.
     """
-    def __init__(self, interval: int = 1 << 31, **kwargs):
+    def __init__(self, interval: int = 3, **kwargs):
         super().__init__(**kwargs)
         self.interval = interval
 
@@ -110,7 +110,7 @@ def __call__(self, algo: Algorithm):
 
 class QualityMetrics(ImageQualityCallback, Callback):
     """From https://github.com/SyneRBI/PETRIC/wiki#metrics-and-thresholds"""
-    def __init__(self, reference_image, whole_object_mask, background_mask, interval: int = 1 << 31, **kwargs):
+    def __init__(self, reference_image, whole_object_mask, background_mask, interval: int = 3, **kwargs):
         # TODO: drop multiple inheritance once `interval` included in CIL
         Callback.__init__(self, interval=interval)
         ImageQualityCallback.__init__(self, reference_image, **kwargs)
@@ -205,6 +205,7 @@ class Dataset:
     whole_object_mask: STIR.ImageData | None
     background_mask: STIR.ImageData | None
     voi_masks: dict[str, STIR.ImageData]
+    FOV_mask: STIR.ImageData
     path: PurePath
 
 
@@ -223,6 +224,10 @@ def get_data(srcdir=".", outdir=OUTDIR, sirf_verbosity=0):
     additive_term = STIR.AcquisitionData(str(srcdir / 'additive_term.hs'))
     mult_factors = STIR.AcquisitionData(str(srcdir / 'mult_factors.hs'))
     OSEM_image = STIR.ImageData(str(srcdir / 'OSEM_image.hv'))
+    # Find FOV mask
+    # WARNING: we are currently using Parralelproj with default settings, which uses a cylindrical FOV.
+    # The current code gives identical results to thresholding the sensitivity image (for those settings)
+    FOV_mask = STIR.TruncateToCylinderProcessor().process(OSEM_image.allocate(1))
     kappa = STIR.ImageData(str(srcdir / 'kappa.hv'))
     if (penalty_strength_file := (srcdir / 'penalisation_factor.txt')).is_file():
         penalty_strength = float(np.loadtxt(penalty_strength_file))
@@ -243,27 +248,34 @@ def get_image(fname):
         for voi in (srcdir / 'PETRIC').glob("VOI_*.hv") if voi.stem[4:] not in ('background', 'whole_object')}
 
     return Dataset(acquired_data, additive_term, mult_factors, OSEM_image, prior, kappa, reference_image,
-                   whole_object_mask, background_mask, voi_masks, srcdir.resolve())
+                   whole_object_mask, background_mask, voi_masks, FOV_mask, srcdir.resolve())
 
 
 DATA_SLICES = {
     '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': {}}
+    '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': {}}
 
 if SRCDIR.is_dir():
     # create list of existing data
     # NB: `MetricsWithTimeout` initialises `SaveIters` which creates `outdir`
     data_dirs_metrics = [
         (SRCDIR / "Siemens_mMR_NEMA_IQ", OUTDIR / "mMR_NEMA",
          [MetricsWithTimeout(outdir=OUTDIR / "mMR_NEMA", **DATA_SLICES['Siemens_mMR_NEMA_IQ'])]),
+        (SRCDIR / "Siemens_mMR_NEMA_IQ_lowcounts", OUTDIR / "mMR_NEMA_lowcounts",
+         [MetricsWithTimeout(outdir=OUTDIR / "mMR_NEMA_lowcounts", **DATA_SLICES['Siemens_mMR_NEMA_IQ_lowcounts'])]),
         (SRCDIR / "NeuroLF_Hoffman_Dataset", OUTDIR / "NeuroLF_Hoffman",
          [MetricsWithTimeout(outdir=OUTDIR / "NeuroLF_Hoffman", **DATA_SLICES['NeuroLF_Hoffman_Dataset'])]),
         (SRCDIR / "Siemens_Vision600_thorax", OUTDIR / "Vision600_thorax",
-         [MetricsWithTimeout(outdir=OUTDIR / "Vision600_thorax", **DATA_SLICES['Siemens_Vision600_thorax'])])]
+         [MetricsWithTimeout(outdir=OUTDIR / "Vision600_thorax", **DATA_SLICES['Siemens_Vision600_thorax'])]),
+        (SRCDIR / "Siemens_mMR_ACR", OUTDIR / "mMR_ACR",
+         [MetricsWithTimeout(outdir=OUTDIR / "mMR_ACR", **DATA_SLICES['Siemens_mMR_ACR'])]),
+        (SRCDIR / "Mediso_NEMA_IQ", OUTDIR / "Mediso_NEMA",
+         [MetricsWithTimeout(outdir=OUTDIR / "Mediso_NEMA", **DATA_SLICES['Mediso_NEMA_IQ'])]),
+        (SRCDIR / "GE_DMI3_Torso", OUTDIR / "DMI3_Torso",
+         [MetricsWithTimeout(outdir=OUTDIR / "DMI3_Torso", **DATA_SLICES['GE_DMI3_Torso'])])]
 else:
     log.warning("Source directory does not exist: %s", SRCDIR)
     data_dirs_metrics = [(None, None, [])] # type: ignore
@@ -294,7 +306,7 @@ def get_image(fname):
         metrics_with_timeout.reset() # timeout from now
         algo = Submission(data)
         try:
-            algo.run(np.inf, callbacks=metrics + submission_callbacks)
+            algo.run(np.inf, callbacks=metrics + submission_callbacks, update_objective_interval=np.inf)
         except Exception:
             print_exc(limit=2)
         finally: