Siemens Vision600 Hoffman

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,28 @@
     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"] is not None) and 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(
+    f"srcdir: {srcdir}\ndownloaddir: {downloaddir}\nprocessingdir: {intermediate_data_path}\nwrite_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 +234,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 +266,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

@@ -40,7 +40,8 @@
         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,8 @@
 # logging.basicConfig(level=logging.INFO)
 
 scanID = args['<data_set>']
+num_updates = args['--updates']
+print(f"num_updates {num_updates}")
 
 if not all((SRCDIR.is_dir(), OUTDIR.is_dir())):
     PETRICDIR = Path('~/devel/PETRIC').expanduser()
@@ -53,7 +58,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, outdir=outdir, seconds=3600 * 100)])
 # %%
 fig = plt.figure()
 data_QC.plot_image(algo.get_output(), **settings.slices)

petric.py

@@ -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