Merge branch 'main' of [email protected]:SyneRBI/PETRIC-MaGeZ.git into s…

…imulation

.gitignore

@@ -1,4 +1,5 @@
 data/
+orgdata/
 output/
 output_test/
 output_magez/

SIRF_data_preparation/.gitignore

@@ -0,0 +1 @@
+runs/

SIRF_data_preparation/BSREM_NeuroLF_Hoffman.py → SIRF_data_preparation/NeuroLF_Hoffman_Dataset/BSREM_NeuroLF_Hoffman.py

@@ -1,9 +1,16 @@
-from petric import MetricsWithTimeout, get_data
+from petric import SRCDIR, MetricsWithTimeout, get_data
 from sirf.contrib.BSREM.BSREM import BSREM1
 from sirf.contrib.partitioner import partitioner
+from SIRF_data_preparation.dataset_settings import get_settings
 
-data = get_data(srcdir="./data/NeuroLF_Hoffman_Dataset", outdir="./output/BSREM_NeuroLF_Hoffman")
-num_subsets = 16
+scanID = "NeuroLF_Hoffman_Dataset"
+settings = get_settings(scanID)
+slices = settings.slices
+num_subsets = settings.num_subsets
+outdir = f"./output/{scanID}/BSREM"
+outdir1 = f"./output/{scanID}/BSREM_cont"
+
+data = get_data(srcdir=SRCDIR / scanID, outdir=outdir)
 data_sub, acq_models, obj_funs = partitioner.data_partition(data.acquired_data, data.additive_term, data.mult_factors,
                                                             num_subsets, initial_image=data.OSEM_image)
 # WARNING: modifies prior strength with 1/num_subsets (as currently needed for BSREM implementations)
@@ -13,5 +20,6 @@
     f.set_prior(data.prior)
 
 algo = BSREM1(data_sub, obj_funs, initial=data.OSEM_image, initial_step_size=.3, relaxation_eta=.01,
-              update_objective_interval=10)
-algo.run(5000, callbacks=[MetricsWithTimeout(transverse_slice=72)])
+              update_objective_interval=80)
+# %%
+algo.run(15000, callbacks=[MetricsWithTimeout(**slices, outdir=outdir, seconds=3600 * 100)])

SIRF_data_preparation/NeuroLF_Hoffman_Dataset/NeuroLF_VOIs.py

@@ -0,0 +1,97 @@
+# %%
+import os
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+from scipy.ndimage import binary_erosion
+
+import sirf.STIR as STIR
+from petric import SRCDIR
+from SIRF_data_preparation import data_QC
+
+# %%
+# %matplotlib ipympl
+# %%
+STIR.AcquisitionData.set_storage_scheme('memory')
+# 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)
+
+
+# %%
+def plot_image(image, vmin=0, vmax=None):
+    if vmax is None:
+        vmax = image.max()
+    plt.figure()
+    plt.subplot(131)
+    plt.imshow(image.as_array()[im_slice, :, :], vmin=vmin, vmax=vmax)
+    plt.subplot(132)
+    plt.imshow(image.as_array()[:, cor_slice, :], vmin=vmin, vmax=vmax)
+    plt.subplot(133)
+    plt.imshow(image.as_array()[:, :, sag_slice], vmin=vmin, vmax=vmax)
+    plt.show()
+
+
+# %% read
+if not SRCDIR.is_dir():
+    PETRICDIR = Path('~/devel/PETRIC').expanduser()
+    SRCDIR = PETRICDIR / 'data'
+datadir = str(SRCDIR / 'NeuroLF_Hoffman_Dataset') + '/'
+
+# %%
+OSEM_image = STIR.ImageData(datadir + 'OSEM_image.hv')
+im_slice = 72
+cor_slice = OSEM_image.dimensions()[1] // 2
+sag_slice = OSEM_image.dimensions()[2] // 2
+cmax = OSEM_image.max()
+
+# %%
+reference_image = STIR.ImageData(datadir + 'reference_image.hv')
+# %% read in original VOIs
+whole_object_mask = STIR.ImageData(datadir + 'whole_phantom.hv')
+# whole_object_mask = STIR.ImageData(datadir+'VOI_whole_object.hv')
+# from README.md
+# 1: ventricles, 2: white matter, 3: grey matter
+allVOIs = STIR.ImageData(datadir + 'vois_ventricles_white_grey.hv')
+# %% create PETRIC VOIs
+background_mask = whole_object_mask.clone()
+background_mask.fill(binary_erosion(allVOIs.as_array() == 2, iterations=2))
+VOI1_mask = whole_object_mask.clone()
+VOI1_mask.fill(allVOIs.as_array() == 1)
+VOI2_mask = whole_object_mask.clone()
+VOI2_mask.fill(binary_erosion(allVOIs.as_array() == 2, iterations=1))
+VOI3_mask = whole_object_mask.clone()
+VOI3_mask.fill(allVOIs.as_array() == 3)
+# %% write PETRIC VOIs
+whole_object_mask.write(datadir + 'PETRIC/VOI_whole_object.hv')
+background_mask.write(datadir + 'PETRIC/VOI_background.hv')
+VOI1_mask.write(datadir + 'PETRIC/VOI_ventricles.hv')
+VOI2_mask.write(datadir + 'PETRIC/VOI_WM.hv')
+VOI3_mask.write(datadir + 'PETRIC/VOI_GM.hv')
+# %%
+VOIs = (whole_object_mask, background_mask, VOI1_mask, VOI2_mask, VOI3_mask)
+# %%
+[plot_image(VOI) for VOI in VOIs]
+# %%
+[plot_image(VOI) for VOI in VOIs]
+
+
+# %%
+def VOI_mean(image, VOI):
+    return float((image * VOI).sum() / VOI.sum())
+
+
+# %%
+[VOI_mean(OSEM_image, VOI) for VOI in VOIs]
+# %%
+[VOI_mean(reference_image, VOI) for VOI in VOIs]
+
+# %%
+os.chdir(datadir)
+# %%
+data_QC.main()
+# %%
+voidir = Path(datadir) / 'PETRIC'
+[str(voi)[:-3] for voi in voidir.glob("VOI_*.hv")]
+# %%

SIRF_data_preparation/README.md

@@ -1,23 +1,27 @@
 # Functions/classes for PETRIC
 
-## Functions to create images that are provided to participants
-
-- `create_initial_images.py`: functions+script to run OSEM and compute the "kappa" image from existing data
-- `BSREM_common.py`: functions to run BSREM with various callbacks
-- `BSREM_*.py`: functions with specific settings for a particular data-set
-
 ## Utility functions to prepare data for the Challenge
 
 Participants should never have to use these (unless you want to create your own dataset).
 
-- `data_utilities.py`: functions to use sirf.STIR to output prompts/mult_factors and additive_term
-  and handle Siemens data
+- `create_initial_images.py`: functions+script to run OSEM and compute the "kappa" image from existing data
 - `data_QC.py`: generates plots for QC
-- Siemens mMR NEMA IQ data (on Zenodo)
-  - `download_Siemens_mMR_NEMA_IQ.py`: download and extract
-  - `prepare_mMR_NEMA_IQ_data.py`: prepare the data (prompts etc)
+- `plot_BSREM_metrics.py`: plot objective functions/metrics after a BSREM run
+- `run_BSREM.py` and `run_OSEM.py`: scripts to run these algorithms for a dataset
 
 ## Helpers
 
+- `data_utilities.py`: functions to use sirf.STIR to output prompts/mult_factors and additive_term
+  and handle Siemens data
 - `evaluation_utilities.py`: reading/plotting helpers for values of the objective function and metrics
 - `PET_plot_functions.py`: plotting helpers
+- `dataset_settings.py`: settings for display of good slices, subsets etc
+
+## Sub-folders per data-set
+
+These contain files specific to the data-set, e.g. for downloading, VOI conversion, settings used for recon, etc.
+Warning: this will be messy, and might be specific to whoever prepared this data-set. For instance,
+for the Siemens mMR NEMA IQ data (on Zenodo):
+- `download_Siemens_mMR_NEMA_IQ.py`: download and extract
+- `prepare_mMR_NEMA_IQ_data.py`: prepare the data (prompts etc)
+- `BSREM_*.py`: functions with specific settings for a particular data-set

SIRF_data_preparation/Siemens_Vision600_thorax/BSREM_Vision600_thorax.py

@@ -0,0 +1,39 @@
+import os
+
+import sirf.STIR as STIR
+from petric import SRCDIR, MetricsWithTimeout, get_data
+from sirf.contrib.BSREM.BSREM import BSREM1
+from sirf.contrib.partitioner import partitioner
+from SIRF_data_preparation.dataset_settings import get_settings
+
+scanID = 'Siemens_Vision600_thorax'
+settings = get_settings(scanID)
+slices = settings.slices
+num_subsets = settings.num_subsets
+outdir = f"./output/{scanID}/BSREM"
+outdir1 = f"./output/{scanID}/BSREM_cont"
+
+data = get_data(srcdir=SRCDIR / scanID, outdir=outdir)
+data_sub, acq_models, obj_funs = partitioner.data_partition(data.acquired_data, data.additive_term, data.mult_factors,
+                                                            num_subsets, 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))
+data.prior.set_up(data.OSEM_image)
+for f in obj_funs: # add prior evenly to every objective function
+    f.set_prior(data.prior)
+
+OSEM_image = data.OSEM_image
+# clean-up some data, now that we have the subsets
+del data
+try:
+    image1000 = STIR.ImageData(os.path.join(outdir, 'iter_1000.hv'))
+except Exception:
+    algo = BSREM1(data_sub, obj_funs, initial=OSEM_image, initial_step_size=.3, relaxation_eta=.01,
+                  update_objective_interval=5)
+    algo.run(1005, callbacks=[MetricsWithTimeout(seconds=3600 * 34, outdir=outdir)])
+
+# restart
+
+algo = BSREM1(data_sub, obj_funs, initial=image1000, initial_step_size=.3, relaxation_eta=.01,
+              update_objective_interval=5)
+algo.run(9000, callbacks=[MetricsWithTimeout(seconds=3600 * 34, outdir=outdir1)])

SIRF_data_preparation/Siemens_Vision600_thorax/Siemens_Vision600_thorax_VOI_prep.sh

@@ -0,0 +1,21 @@
+#! /bin/bash
+# get rid of offset in niftis, as they are currently not in the reconstructed images
+for f in *nii.gz; do
+   base=${f%%.nii.gz}
+   stir_math $base.hv $base.nii.gz
+   sed -i -e '/first pixel offset/ d' \
+     -e '/!INTERFILE/a !imaging modality := PT' \
+    $base.hv
+    # patient position is currently not in reconstructed images
+    # -e '/type of data/a patient orientation := feet_in' \
+    # -e '/type of data/a patient rotation := supine' \
+    invert_axis x ${base}_flip_x.hv ${base}.hv
+    invert_axis z ${base}.hv ${base}_flip_x.hv
+    rm ${base}_flip_x.*v
+done
+
+stir_math ../PETRIC/VOI_whole_object.hv phantom.hv
+stir_math ../PETRIC/VOI_background.hv background.hv
+for r in 1 4 7; do
+   stir_math ../PETRIC/VOI_${r}.hv VOI${r}_ds.hv
+done

SIRF_data_preparation/Siemens_mMR_ACR/README.md

@@ -4,9 +4,29 @@ ugly and temporary
 Steps to follow:
 1. `python SIRF_data_preparation/Siemens_mMR_ACR/download.py`
 2. `python SIRF_data_preparation/Siemens_mMR_ACR/prepare.py` (As there is no useful mumap, this will **not** do attenuation and scatter estimation)
-3. `python SIRF_data_preparation/create_initial_images.py data/Siemens_mMR_ACR/final --template-image=None` (Reconstruct at full FOV size)
-4. `mv data/Siemens_mMR_ACR/final/OSEM_image.* data/Siemens_mMR_ACR/processing` (this is really an NAC image)
-5. `python SIRF_data_preparation/Siemens_mMR_ACR/register_mumap.py` (output is data/Siemens_mMR_ACR/processing/reg_mumap.hv etc)
-6. `python SIRF_data_preparation/Siemens_mMR_ACR/prepare.py --end 200` (As there is now a useful mumap, this will now do attenuation and scatter estimation)
-7. `python SIRF_data_preparation/create_initial_images.py -s 200 data/Siemens_mMR_ACR/final` (ideally add `--template-image some_VOI`)
-However, registration failed, so this needs a manual intervention step. Output of that should be data/Siemens_mMR_ACR/processing/reg_mumap.hv
+3. `python SIRF_data_preparation/create_initial_images.py data/Siemens_mMR_ACR/final --template_image=None` (Reconstruct at full FOV size)
+4. `mv data/Siemens_mMR_ACR/final/OSEM_image.* data/Siemens_mMR_ACR/processing` (this is really an NAC image and ideally would be renamed)
+5. `python SIRF_data_preparation/Siemens_mMR_ACR/register_mumap.py` (output is orgdata/Siemens_mMR_ACR/processing/reg_mumap.hv etc)
+6. However, registration failed, so this needs a manual intervention step:
+   a. send data to Pawel who registers.
+   b. Download from https://drive.google.com/file/d/13eGUL3fwdpCiWLjoP1Lhkd8xQOyV_AAZ/view?usp=sharing
+   c. `cd  data/Siemens_mMR_ACR; unzip ../downloads/ACR_STIR_output.zip`
+   d. copy data and rename (currently using `stir_math` executable here)
+      ```sh
+      stir_math orgdata/Siemens_mMR_ACR/processing/reg_mumap.hv   orgdata/Siemens_mMR_ACR/output/acr-mumap-complete.nii.gz
+      ```
+   e. uncompress .gz VOI files as STIR isn't happy with them (probably not necessary)
+      ```sh
+      for f in orgdata/Siemens_mMR_ACR/output/sampling_masks/*gz; do gunzip $f; done
+      ```
+6. Add hardware-mumap
+   ```sh
+   stir_math --accumulate  orgdata/Siemens_mMR_ACR/processing/reg_mumap.hv  orgdata/Siemens_mMR_ACR/output/ACR_hardware-to-STIR.nii.gz
+   ```
+7. 11. `python SIRF_data_preparation/run_OSEM.py Siemens_mMR_ACR`--end 200` (As there is now a useful mumap, this will now do attenuation and scatter estimation)
+8. `python SIRF_data_preparation/create_initial_images.py data/Siemens_mMR_ACR --template_image=../../orgdata/Siemens_mMR_ACR/output/sampling_masks/acr-all-sampling-0-2mm_dipy.nii`
+9. `python SIRF_data_preparation/data_QC.py --srcdir='data/Siemens_mMR_ACR' --transverse_slice=99`
+10. edit `SIRF_data_prepatation/dataset_settings.py` for subsets etc. edit `OSEM_image.hv` to add modality, radionuclide and duration info which got lost.
+11. `python SIRF_data_preparation/Siemens_mMR_ACR/VOI_prep.py`
+12. `python SIRF_data_preparation/run_OSEM.py Siemens_mMR_ACR`
+13. `python SIRF_data_preparation/run_BSREM.py Siemens_mMR_ACR`

SIRF_data_preparation/Siemens_mMR_ACR/VOI_prep.py

@@ -0,0 +1,95 @@
+# %% file to prepare VOIs from the original ACR NEMA VOIs supplied by Pawel Markiewicz
+# %%
+import os
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import sirf.STIR as STIR
+import SIRF_data_preparation.data_QC as data_QC
+from SIRF_data_preparation.data_utilities import the_data_path, the_orgdata_path
+from SIRF_data_preparation.dataset_settings import get_settings
+
+# %%
+scanID = 'Siemens_mMR_ACR'
+org_VOI_path = the_orgdata_path(scanID, 'output', 'sampling_masks')
+data_path = the_data_path(scanID)
+output_path = the_data_path(scanID, 'PETRIC')
+os.makedirs(output_path, exist_ok=True)
+settings = get_settings(scanID)
+slices = settings.slices
+# %%
+OSEM_image = STIR.ImageData(os.path.join(data_path, 'OSEM_image.hv'))
+cmax = OSEM_image.max()
+# %% read in original VOIs
+orgVOIs = STIR.ImageData(os.path.join(org_VOI_path, 'acr-all-sampling-0-2mm_dipy.nii'))
+data_QC.plot_image(orgVOIs, **slices)
+# %%
+orgVOIs_arr = orgVOIs.as_array()
+print(np.unique(orgVOIs_arr))
+
+
+# %% output
+# [  0,  10-29,  40-81,  90- 101, 300- 317]
+# %%
+def plotMask(mask):
+    plt.figure()
+    plt.subplot(141)
+    plt.imshow(mask[:, 109, :])
+    plt.subplot(142)
+    plt.imshow(mask[85, :, :])
+    plt.subplot(143)
+    plt.imshow(mask[99, :, :])
+    plt.subplot(144)
+    plt.imshow(mask[40, :, :])
+
+
+# %% background mask: slices in uniform part (idx 300-317), taking slightly eroded mask
+mask = (orgVOIs_arr >= 300) * (orgVOIs_arr < 316)
+plotMask(mask)
+# %%
+background_mask = OSEM_image.clone()
+background_mask.fill(mask)
+# %% 6 cylinder masks
+mask = (orgVOIs_arr >= 10) * (orgVOIs_arr < 101)
+plotMask(mask)
+# %% cold cylinder mask
+mask = (orgVOIs_arr >= 90) * (orgVOIs_arr < 100)
+plotMask(mask)
+cold_cyl_mask = OSEM_image.clone()
+cold_cyl_mask.fill(mask)
+# %% faint hot cylinder mask
+mask = (orgVOIs_arr >= 20) * (orgVOIs_arr < 29)
+plotMask(mask)
+hot_cyl_mask = OSEM_image.clone()
+hot_cyl_mask.fill(mask)
+# %% Jasczcak part, not available so derive from "background mask"
+mask = (orgVOIs_arr >= 300) * (orgVOIs_arr < 316)
+# get single cylinder from a slice in the background
+oneslice = mask[85, :, :].copy()
+mask[35:45, :, :] = oneslice[:, :]
+mask[46:127, :, :] = False
+plotMask(mask)
+rods_mask = OSEM_image.clone()
+rods_mask.fill(mask)
+# %% whole phantom
+mask = OSEM_image.as_array() > cmax / 20
+plotMask(mask)
+whole_object_mask = OSEM_image.clone()
+whole_object_mask.fill(mask)
+# %%
+plotMask(OSEM_image.as_array())
+# %% write PETRIC VOIs
+whole_object_mask.write(os.path.join(output_path, 'VOI_whole_object.hv'))
+background_mask.write(os.path.join(output_path, 'VOI_background.hv'))
+cold_cyl_mask.write(os.path.join(output_path, 'VOI_cold_cylinder.hv'))
+hot_cyl_mask.write(os.path.join(output_path, 'VOI_hot_cylinder.hv'))
+rods_mask.write(os.path.join(output_path, 'VOI_rods.hv'))
+
+# %%
+VOIs = (whole_object_mask, background_mask, cold_cyl_mask, hot_cyl_mask, rods_mask)
+# %%
+[data_QC.plot_image(VOI, **slices) for VOI in VOIs]
+# %%
+data_QC.VOI_checks(['VOI_whole_object', 'VOI_background', 'VOI_cold_cylinder', 'VOI_hot_cylinder', 'VOI_rods'],
+                   OSEM_image, srcdir=output_path, **slices)