Merge pull request #102 from KrisThielemans/Mediso_NEMA_IQ

Mediso NEMA IQ data

SIRF_data_preparation/Mediso_NEMA_IQ/README.md

@@ -0,0 +1,36 @@
+# Data from the Mediso AnyScan at NPL
+Initial processing steps
+```sh
+orgpath=/mnt/share/petric-wip/NEMA_Challenge
+#orgpath=~/devel/PETRIC/orgdata/
+# cd $orgpath
+#wget https://nplanywhere.npl.co.uk/userportal/?v=4.6.1#/shared/public/nrKe5vWynvyMsp7m/e9e3d4b8-8f83-47be-bf2d-17a96a7a8aac
+#unzip NEMA_Challenge.zip
+cd ~/devel/PETRIC/data
+mkdir -p  Mediso_NEMA_IQ
+cd Mediso_NEMA_IQ/
+# trim sinograms to avoid "corner" problems in mult_factors
+# TODO for next data: use 30 (as some problems remain)
+for f in additive_term.hs  mult_factors.hs  prompts.hs; do
+   SSRB -t 20 $f $orgpath/sinograms/$f; 
+done
+# alternative if we don't need to trim
+#cp -rp $orgpath/sinograms/* .
+# get rid of NaNs
+python prepare.py
+# now handle VOIs
+mkdir PETRIC
+cp -rp $orgpath/VOIs/* PETRIC
+# cp -rp $orgpath/README.md .
+cd PETRIC
+stir_math VOI_background.hv VOI_backgroung.hv
+rm VOI_backgroung.*
+rm *ahv
+cd ..
+python ../../SIRF_data_preparation/create_initial_images.py --template_image=PETRIC/VOI_whole_object.hv .
+ ```
+ I needed to fix the OSEM header (manual):
+ ```
+ !imaging modality := PT
+ ```
+python ../../SIRF_data_preparation/data_QC

SIRF_data_preparation/Mediso_NEMA_IQ/prepare.py

@@ -0,0 +1,10 @@
+# Set NaNs to zero in additive_term
+import sirf.STIR
+import numpy as np
+additive = sirf.STIR.AcquisitionData('additive_term.hs')
+add_arr = additive.as_array()
+add_arr = np.nan_to_num(add_arr)
+new_add = additive.clone()
+new_add.fill(add_arr)
+new_add.write('additive_term.hs')
+

SIRF_data_preparation/data_QC.py

@@ -11,12 +11,13 @@
   --transverse_slice=<i>  idx [default: -1]
   --coronal_slice=<c>    idx [default: -1]
   --sagittal_slice=<s>   idx [default: -1]
+  --dataset=<name>       dataset name. if set, it is used to override default slices
 
 Note that -1 one means to use middle of image
 """
 # Copyright 2024 University College London
 # Licence: Apache-2.0
-__version__ = '0.2.0'
+__version__ = '0.3.0'
 
 import os
 import os.path
@@ -29,6 +30,7 @@
 from scipy import ndimage
 
 import sirf.STIR as STIR
+from SIRF_data_preparation.dataset_settings import get_settings
 
 STIR.AcquisitionData.set_storage_scheme('memory')
 
@@ -119,7 +121,10 @@ def VOI_checks(allVOInames, OSEM_image=None, reference_image=None, srcdir='.', *
             print(f"VOI {VOIname} does not exist")
             continue
         VOI = STIR.ImageData(filename)
-        COM = np.rint(ndimage.center_of_mass(VOI.as_array()))
+        VOI_arr = VOI.as_array()
+        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")
         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]))
@@ -149,13 +154,21 @@ def VOI_checks(allVOInames, OSEM_image=None, reference_image=None, srcdir='.', *
 
 def main(argv=None):
     args = docopt(__doc__, argv=argv, version=__version__)
+    dataset = args['--dataset']
     srcdir = args['--srcdir']
     skip_sino_profiles = args['--skip_sino_profiles']
     slices = {}
     slices["transverse_slice"] = literal_eval(args['--transverse_slice'])
     slices["coronal_slice"] = literal_eval(args['--coronal_slice'])
     slices["sagittal_slice"] = literal_eval(args['--sagittal_slice'])
 
+    if (dataset):
+        settings = get_settings(dataset)
+        for key in slices.keys():
+            if slices[key] == -1 and key in settings.slices:
+                slices[key] = settings.slices[key]
+    print(slices)
+
     if not skip_sino_profiles:
         acquired_data = STIR.AcquisitionData(os.path.join(srcdir, 'prompts.hs'))
         additive_term = STIR.AcquisitionData(os.path.join(srcdir, 'additive_term.hs'))

SIRF_data_preparation/dataset_settings.py

@@ -19,6 +19,9 @@ def get_settings(scanID: str):
     elif scanID == 'NeuroLF_Hoffman_Dataset':
         slices = {'transverse_slice': 72}
         num_subsets = 16
+    elif scanID == 'Mediso_NEMA_IQ':
+        slices = {'transverse_slice': 22, 'coronal_slice': 89, 'sagittal_slice': 66}
+        num_subsets = 12
     else:                                                       # Vision
         slices = {}
         num_subsets = 5

SIRF_data_preparation/plot_BSREM_metrics.py

@@ -1,4 +1,5 @@
-"""Preliminary file to check evolution of metrics as well as pass_index"""
+#%%
+# """Preliminary file to check evolution of metrics as well as pass_index"""
 # %load_ext autoreload
 # %autoreload 2
 from pathlib import Path
@@ -21,7 +22,8 @@
 
 # scanID = 'Siemens_Vision600_thorax'
 # scanID = 'NeuroLF_Hoffman_Dataset'
-scanID = 'Siemens_mMR_NEMA_IQ'
+# scanID = 'Siemens_mMR_NEMA_IQ'
+scanID = 'Mediso_NEMA_IQ'
 
 srcdir = SRCDIR / scanID
 outdir = OUTDIR / scanID
@@ -69,16 +71,19 @@
     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
 last_iteration = int(objs[-1, 0] + .5)
-iteration_interval = int(objs[-1, 0] - objs[-2, 0] + .5)
+# find interval(don't use last value, as that interval can be smaller)
+iteration_interval = int(objs[-2, 0] - objs[-3, 0] + .5)
 if datadir1.is_dir():
     last_iteration = int(objs0[-1, 0] + .5)
     iteration_interval = int(objs0[-1, 0] - objs0[-2, 0] + .5) * 2
 # %%
 iters = range(0, last_iteration, iteration_interval)
 m = get_metrics(qm, iters, srcdir=datadir)
 # %%
-OSEMiters = range(0, 400, 20)
+OSEMobjs = read_objectives(OSEMdir)
+OSEMiters = OSEMobjs[:, 0].astype(int)
 OSEMm = get_metrics(qm, OSEMiters, srcdir=OSEMdir)
 # %%
 fig = plt.figure()
@@ -107,7 +112,7 @@
     fig.savefig(outdir / f'{scanID}_metrics_BSREM.png')
 
 # %%
-idx = pass_index(m, numpy.array([.01, .01, .005, .005, .005]), 10)
+idx = pass_index(m, numpy.array([.01, .01] + [.005 for i in range(len(data.voi_masks))]), 10)
 iter = iters[idx]
 print(iter)
 image = STIR.ImageData(str(datadir / f"iter_{iter:04d}.hv"))

SIRF_data_preparation/run_OSEM.py

@@ -42,7 +42,7 @@
 data = get_data(srcdir=srcdir, outdir=outdir)
 # %%
 algo = main_OSEM.Submission(data, settings.num_subsets, update_objective_interval=20)
-algo.run(20, callbacks=[MetricsWithTimeout(**settings.slices, seconds=5000, outdir=outdir)])
+algo.run(400, callbacks=[MetricsWithTimeout(**settings.slices, seconds=5000, outdir=outdir)])
 # %%
 fig = plt.figure()
 data_QC.plot_image(algo.get_output(), **settings.slices)

petric.py

@@ -141,7 +141,7 @@ def keys(self):
 
 class MetricsWithTimeout(cil_callbacks.Callback):
     """Stops the algorithm after `seconds`"""
-    def __init__(self, seconds=600, outdir=OUTDIR, transverse_slice=None, coronal_slice=None, **kwargs):
+    def __init__(self, seconds=600, outdir=OUTDIR, transverse_slice=None, coronal_slice=None, sagittal_slice=None, **kwargs):
         super().__init__(**kwargs)
         self._seconds = seconds
         self.callbacks = [