Merge pull request #134 from KrisThielemans/ZrNEMAIQ

Siemens_Vision600_ZrNEMAIQ

SIRF_data_preparation/Siemens_Vision600_ZrNEMAIQ/prep_VOIs.py

@@ -0,0 +1,69 @@
+# %%
+# %load_ext autoreload
+# %%
+from pathlib import Path
+
+import numpy as np
+from scipy.ndimage import binary_erosion
+
+import sirf.STIR as STIR
+from SIRF_data_preparation import data_QC
+from SIRF_data_preparation.data_utilities import the_data_path, the_orgdata_path
+from SIRF_data_preparation.dataset_settings import get_settings
+
+# %%
+# %autoreload 2
+
+# %%
+# set-up redirection of STIR messages to files
+# _ = STIR.MessageRedirector('BSREM_info.txt', 'BSREM_warnings.txt', 'BSREM_errors.txt')
+# fewer messages from STIR and SIRF (set to higher value to diagnose have problems)
+STIR.set_verbosity(0)
+
+# %% read
+scanID = 'Siemens_Vision600_ZrNEMAIQ'
+datadir = Path(the_data_path(scanID))
+intermediate_datadir = Path(the_orgdata_path(scanID, 'processing'))
+settings = get_settings(scanID)
+slices = settings.slices
+# %%
+OSEM_image = STIR.ImageData(str(datadir / 'OSEM_image.hv'))
+cmax = np.percentile(OSEM_image.as_array(), 99) / .99
+# %% read in original VOIs
+orgVOIs = []
+for i in range(1, 8):
+    orgVOIs.append(STIR.ImageData(str(intermediate_datadir / f"S{i}.hv")))
+    data_QC.plot_image(orgVOIs[i - 1], **slices)
+# %%
+reference_image = STIR.ImageData(str(datadir / 'PETRIC' / 'reference_image.hv'))
+data_QC.plot_image(reference_image, **slices)
+# %% create PETRIC VOIs
+whole_object_mask = reference_image.clone()
+ref_arr = reference_image.as_array()
+whole_object_mask.fill(binary_erosion(ref_arr > np.percentile(ref_arr, 5), iterations=2))
+data_QC.plot_image(whole_object_mask, **slices)
+# %%
+background_mask = orgVOIs[6]
+VOI1_mask = orgVOIs[4]
+VOI2_mask = orgVOIs[2]
+VOI3_mask = orgVOIs[1]
+# %% write PETRIC VOIs
+whole_object_mask.write(str(datadir / 'PETRIC' / 'VOI_whole_object.hv'))
+VOI1_mask.write(str(datadir / 'PETRIC' / 'VOI_sphere5.hv'))
+VOI2_mask.write(str(datadir / 'PETRIC' / 'VOI_sphere3.hv'))
+VOI3_mask.write(str(datadir / 'PETRIC' / 'VOI_sphere2.hv'))
+# %% 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()
+z_size = background_arr.shape[0]
+background_arr = np.concatenate((
+    background_arr[0:35, :, :] * 0,
+    background_arr[0:-35, :, :],
+))
+background_mask.fill(background_arr)
+background_mask.write(str(datadir / 'PETRIC' / 'VOI_background.hv'))
+
+# %%
+VOInames = ('VOI_whole_object', 'VOI_background', 'VOI_sphere5', 'VOI_sphere3', 'VOI_sphere2')
+data_QC.VOI_checks(VOInames, OSEM_image=OSEM_image, reference_image=reference_image, srcdir=str(datadir / 'PETRIC'))

SIRF_data_preparation/data_QC.py

@@ -73,10 +73,10 @@ def plot_image(image, save_name=None, transverse_slice=-1, coronal_slice=-1, sag
         coronal_slice = image.dimensions()[1] // 2
     if sagittal_slice < 0:
         sagittal_slice = image.dimensions()[2] // 2
+    arr = image.as_array()
     if vmax is None:
-        vmax = image.max()
+        vmax = np.percentile(arr, 99.995)
 
-    arr = image.as_array()
     alpha_trans = None
     alpha_cor = None
     alpha_sag = None
@@ -147,6 +147,10 @@ def VOI_checks(allVOInames, OSEM_image=None, reference_image=None, srcdir='.', *
         plt.figure()
         plot_image(VOI, save_name=prefix, vmin=0, vmax=1, transverse_slice=int(COM[0]), coronal_slice=int(COM[1]),
                    sagittal_slice=int(COM[2]))
+        if OSEM_image is not None:
+            plt.figure()
+            plot_image(OSEM_image, alpha=(VOI+.5) / 1.5, save_name=prefix + "_and_OSEM", transverse_slice=int(COM[0]),
+                       coronal_slice=int(COM[1]), sagittal_slice=int(COM[2]))
 
         # construct transparency image
         if VOIname == 'VOI_whole_object':

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

SIRF_data_preparation/plot_BSREM_metrics.py

@@ -2,6 +2,7 @@
 # """Preliminary file to check evolution of metrics as well as pass_index"""
 # %load_ext autoreload
 # %autoreload 2
+# %%
 from pathlib import Path
 
 import matplotlib.pyplot as plt
@@ -14,6 +15,7 @@
 from SIRF_data_preparation.dataset_settings import get_settings
 from SIRF_data_preparation.evaluation_utilities import get_metrics, plot_metrics, read_objectives
 
+# %%
 if not all((SRCDIR.is_dir(), OUTDIR.is_dir())):
     PETRICDIR = Path('~/devel/PETRIC').expanduser()
     OUTDIR = PETRICDIR / 'output'
@@ -25,7 +27,8 @@
 # scanID = 'Siemens_mMR_NEMA_IQ'
 # scanID = 'Mediso_NEMA_IQ'
 scanID = 'Siemens_Vision600_Hoffman'
-
+# scanID = 'NeuroLF_Esser_Dataset'
+# scanID = 'Siemens_Vision600_ZrNEMAIQ'
 srcdir = the_data_path(scanID)
 outdir = OUTDIR / scanID
 OSEMdir = outdir / 'OSEM'
@@ -38,7 +41,7 @@
 settings = get_settings(scanID)
 slices = settings.slices
 
-cmax = OSEM_image.max()
+cmax = numpy.percentile(OSEM_image.as_array(), 99.995)
 # %%
 image = data_QC.plot_image_if_exists(str(datadir / 'iter_final'), **slices, vmax=cmax)
 # image2=STIR.ImageData(datadir+'iter_14000.hv')
@@ -66,8 +69,11 @@
 fig.savefig(outdir / f'{scanID}_BSREM_objectives_last.png')
 
 # %%
-data = get_data(srcdir=srcdir, outdir=None)
-if datadir1.is_dir():
+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'))
@@ -82,7 +88,9 @@
     iteration_interval = int(objs0[-1, 0] - objs0[-2, 0] + .5) * 2
 # %%
 iters = range(0, last_iteration, iteration_interval)
+print('GETMETRICS')
 m = get_metrics(qm, iters, srcdir=datadir)
+print('DONE')
 # %%
 OSEMobjs = read_objectives(OSEMdir)
 OSEMiters = OSEMobjs[:, 0].astype(int)

SIRF_data_preparation/run_BSREM.py

@@ -8,10 +8,13 @@
 
 Options:
   --updates=<u>               number of updates to run [default: 15000]
-  --initial_image=<filename>  optional initial image, normally the OSEM_image. If specified,
+  --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]
 """
 # Copyright 2024 Rutherford Appleton Laboratory STFC
 # Copyright 2024 University College London
@@ -37,8 +40,10 @@
 scanID = args['<data_set>']
 num_updates = int(args['--updates'])
 initial_image = args['--initial_image']
+num_subsets = args['--num_subsets']
 initial_step_size = float(args['--initial_step_size'])
 relaxation_eta = float(args['--relaxation_eta'])
+interval = int(args['--interval'])
 
 if not all((SRCDIR.is_dir(), OUTDIR.is_dir())):
     PETRICDIR = Path('~/devel/PETRIC').expanduser()
@@ -50,25 +55,29 @@
 # log.info("Finding files in %s", srcdir)
 
 settings = get_settings(scanID)
+num_subsets = settings.num_subsets if num_subsets is None else int(num_subsets)
 
 data = get_data(srcdir=srcdir, outdir=outdir)
 if initial_image is None:
+    initial_image_name = "OSEM"
     initial_image = data.OSEM_image
     outdir = outdir / "BSREM"
 else:
+    initial_image_name = initial_image
     initial_image = STIR.ImageData(initial_image)
     outdir = outdir / "BSREM_cont"
 
 print("Penalisation factor:", data.prior.get_penalisation_factor())
-print("num_subsets:", settings.num_subsets)
+print("num_subsets:", num_subsets)
 print("num_updates:", num_updates)
-print("initial_image:", initial_image)
+print("initial_image:", initial_image_name)
 print("outdir:", outdir)
 print("initial_step_size:", initial_step_size)
 print("relaxation_eta:", relaxation_eta)
+print("interval:", interval)
 
 data_sub, acq_models, obj_funs = partitioner.data_partition(data.acquired_data, data.additive_term, data.mult_factors,
-                                                            settings.num_subsets, mode="staggered",
+                                                            num_subsets, mode="staggered",
                                                             initial_image=data.OSEM_image)
 # WARNING: modifies prior strength with 1/num_subsets (as currently needed for BSREM implementations)
 data.prior.set_penalisation_factor(data.prior.get_penalisation_factor() / len(obj_funs))
@@ -77,9 +86,10 @@
     f.set_prior(data.prior)
 
 algo = BSREM1(data_sub, obj_funs, initial=initial_image, initial_step_size=initial_step_size,
-              relaxation_eta=relaxation_eta, update_objective_interval=80)
+              relaxation_eta=relaxation_eta, update_objective_interval=interval)
 # %%
-algo.run(num_updates, callbacks=[MetricsWithTimeout(**settings.slices, interval=80, outdir=outdir, seconds=3600 * 100)])
+algo.run(num_updates,
+         callbacks=[MetricsWithTimeout(**settings.slices, interval=interval, outdir=outdir, seconds=3600 * 100)])
 # %%
 fig = plt.figure()
 data_QC.plot_image(algo.get_output(), **settings.slices)

petric.py

@@ -245,7 +245,7 @@ class Dataset:
     path: PurePath
 
 
-def get_data(srcdir=".", outdir=OUTDIR, sirf_verbosity=0):
+def get_data(srcdir=".", outdir=OUTDIR, sirf_verbosity=0, read_sinos=True):
     """
     Load data from `srcdir`, constructs prior and return as a `Dataset`.
     Also redirects sirf.STIR log output to `outdir`, unless that's set to None
@@ -257,9 +257,9 @@ def get_data(srcdir=".", outdir=OUTDIR, sirf_verbosity=0):
     if outdir is not None:
         outdir = Path(outdir)
         _ = STIR.MessageRedirector(str(outdir / 'info.txt'), str(outdir / 'warnings.txt'), str(outdir / 'errors.txt'))
-    acquired_data = STIR.AcquisitionData(str(srcdir / 'prompts.hs'))
-    additive_term = STIR.AcquisitionData(str(srcdir / 'additive_term.hs'))
-    mult_factors = STIR.AcquisitionData(str(srcdir / 'mult_factors.hs'))
+    acquired_data = STIR.AcquisitionData(str(srcdir / 'prompts.hs')) if read_sinos else None
+    additive_term = STIR.AcquisitionData(str(srcdir / 'additive_term.hs')) if read_sinos else None
+    mult_factors = STIR.AcquisitionData(str(srcdir / 'mult_factors.hs')) if read_sinos else None
     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.
@@ -292,8 +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': {}}
+    '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}}
 
 if SRCDIR.is_dir():
     # create list of existing data
@@ -316,7 +317,9 @@ def get_image(fname):
         (SRCDIR / "Siemens_Vision600_Hoffman", OUTDIR / "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'])])]
+         [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'])])]
 else:
     log.warning("Source directory does not exist: %s", SRCDIR)
     data_dirs_metrics = [(None, None, [])] # type: ignore