Merge pull request #88 from SyneRBI/Siemens_mMR_ACR2

Siemens mMR_ACR updates

.gitignore

@@ -1,4 +1,5 @@
 data/
+orgdata/
 output/
 tmp*/
 err*.txt

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,91 @@
+#%% file to prepare VOIs from the original ACR NEMA VOIs supplied by Pawel Markiewicz
+#%%
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+import sys
+import sirf.STIR as STIR
+from sirf.Utilities import examples_data_path
+from scipy.ndimage import binary_erosion
+from pathlib import Path
+from SIRF_data_preparation.data_utilities import the_data_path, the_orgdata_path
+from SIRF_data_preparation.dataset_settings import get_settings
+import SIRF_data_preparation.data_QC as data_QC
+#%%
+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)

SIRF_data_preparation/Siemens_mMR_ACR/download.py

@@ -7,20 +7,16 @@
 from zenodo_get import zenodo_get
 
 # %%
+from SIRF_data_preparation.data_utilities import the_data_path
 from sirf.Utilities import examples_data_path
 
 sirf_data_path = os.path.join(examples_data_path('PET'), 'mMR')
 
 # %% set paths
-this_directory = os.path.dirname(__file__)
-repo_directory = os.path.dirname(os.path.dirname(this_directory))
-challenge_data_path = os.path.join(repo_directory, 'data')
-print('PETRIC data path: %s' % challenge_data_path)
-
-download_path = os.path.join(challenge_data_path, 'downloads')
+download_path = the_data_path('downloads')
 os.makedirs(download_path, exist_ok=True)
 
-dest_path = os.path.join(challenge_data_path, 'Siemens_mMR_ACR')
+dest_path = the_data_path('Siemens_mMR_ACR')
 os.makedirs(dest_path, exist_ok=True)
 print('dest path: %s' % dest_path)
 processing_path = os.path.join(dest_path, 'processing')

SIRF_data_preparation/Siemens_mMR_ACR/prepare.py

@@ -2,12 +2,9 @@
 import logging
 import os
 
-from data_utilities import prepare_challenge_Siemens_data, the_data_path
+from SIRF_data_preparation.data_utilities import prepare_challenge_Siemens_data, the_data_path, the_orgdata_path
 
-this_directory = os.path.dirname(__file__)
-repo_directory = os.path.dirname(this_directory)
-repo_directory = os.path.dirname(repo_directory)
-output_path = os.path.join(repo_directory, 'data')
+scanID = 'Siemens_mMR_ACR'
 
 if __name__ == '__main__':
     parser = argparse.ArgumentParser(description='SyneRBI PETRIC Siemens mMR ACR data preparation script.')
@@ -27,16 +24,15 @@
     end = args.end
 
     if args.raw_data_path is None:
-        data_path = the_data_path('Siemens_mMR_ACR', 'processing')
+        data_path = the_orgdata_path(scanID, 'processing')
     else:
         data_path = args.raw_data_path
 
     data_path = os.path.abspath(data_path)
     logging.debug(f"Raw data path: {data_path}")
 
-    output_path = os.path.join(output_path, 'Siemens_mMR_ACR')
-    intermediate_data_path = os.path.join(output_path, 'processing')
-    output_path = os.path.join(output_path, 'final')
+    intermediate_data_path = the_orgdata_path(scanID, 'processing')
+    output_path = the_data_path(scanID)
 
     os.makedirs(output_path, exist_ok=True)
     os.chdir(output_path)

SIRF_data_preparation/Siemens_mMR_ACR/register_mumap.py

@@ -2,16 +2,11 @@
 
 import sirf.Reg as Reg
 import sirf.STIR as STIR
-
-# %%
-this_directory = os.path.dirname(__file__)
-repo_directory = os.path.dirname(this_directory)
-repo_directory = os.path.dirname(repo_directory)
+from SIRF_data_preparation.data_utilities import the_data_path, the_orgdata_path
 # %% set paths filenames
-root_path = os.path.join(repo_directory, 'data', 'Siemens_mMR_ACR')
-intermediate_data_path = os.path.join(root_path, 'processing')
-output_path = os.path.join(root_path, 'final')
-mumap_filename = os.path.join(root_path, 'ACR_data_design/synth_mumap/acr-complete-umap.nii.gz')
+scanID = 'Siemens_mMR_ACR'
+intermediate_data_path = the_orgdata_path(scanID, 'processing')
+mumap_filename = the_orgdata_path(scanID, 'ACR_data_design/synth_mumap/acr-complete-umap.nii.gz')
 # %% write current OSEM image as Nifti
 NAC_image = STIR.ImageData(os.path.join(intermediate_data_path, 'OSEM_image.hv'))
 NAC_image_filename = os.path.join(intermediate_data_path, 'NAC_image.nii')

SIRF_data_preparation/Siemens_mMR_NEMA_IQ/Siemens_mMR_NEMA_VOIs.py

@@ -77,58 +77,10 @@
 data_QC.plot_image(reference_image, **slices)
 plt.figure()
 data_QC.plot_image(reference_image, **slices, alpha=allVOIs)
-# %%
-
-
-def VOI_checks(allVOInames, OSEM_image=None, reference_image=None, srcdir='.', **kwargs):
-    if len(allVOInames) == 0:
-        return
-    OSEM_VOI_values = []
-    ref_VOI_values = []
-    allVOIs = None
-    VOIkwargs = kwargs.copy()
-    VOIkwargs['vmax'] = 1
-    VOIkwargs['vmin'] = 0
-    for VOIname in allVOInames:
-        filename = os.path.join(srcdir, VOIname + '.hv')
-        if not os.path.isfile(filename):
-            print(f"VOI {VOIname} does not exist")
-            continue
-        VOI = STIR.ImageData(filename)
-        COM = numpy.rint(ndimage.measurements.center_of_mass(VOI.as_array()))
-        print(COM)
-        plt.figure()
-        plot_image(VOI, vmin=0, vmax=1, transverse_slice=int(COM[0]), coronal_slice=int(COM[1]),
-                   sagittal_slice=int(COM[2]))
-
-        # construct transparency image
-        if VOIname == 'VOI_whole_object':
-            VOI /= 2
-        if allVOIs is None:
-            allVOIs = VOI.clone()
-        else:
-            allVOIs += VOI
-        if OSEM_image is not None:
-            OSEM_VOI_values.append(VOI_mean(OSEM_image, VOI))
-        if reference_image:
-            ref_VOI_values.append(VOI_mean(reference_image, VOI))
-    allVOIs /= allVOIs.max()
-
-    if OSEM_image is not None:
-        plt.figure()
-        plot_image(OSEM_image, **kwargs)
-        plt.figure()
-        plot_image(OSEM_image, alpha=allVOIs, **kwargs)
-
 
-# %%
+#%%
 VOI_checks(['VOI_whole_object', 'VOI_sphere5'], OSEM_image, srcdir=os.path.join(datadir, 'PETRIC'), **slices)
-# %%
-OSEM_image
-# %%
-[data_QC.VOI_mean(OSEM_image, VOI) for VOI in VOIs]
-# %%
-[data_QC.VOI_mean(reference_image, VOI) for VOI in VOIs]
-# %%
-data_QC.main(['--srcdir', datadir])
-# %%
+#%%
+[ VOI_mean(OSEM_image, VOI) for VOI in VOIs]
+#%%
+[ VOI_mean(reference_image, VOI) for VOI in VOIs]

SIRF_data_preparation/create_initial_images.py

@@ -7,7 +7,7 @@
   <data_path>  path to data files
 
 Options:
-  -t <template_image>, --template_image=<template_image>  filename of image to use for data sizes [default: VOI.hv]
+  -t <template_image>, --template_image=<template_image>  filename (relative to <data_path>) of image to use for data sizes [default: PETRIC/VOI_whole_phantom.hv]
   -s <xy-size>, --xy-size=<xy-size>  force xy-size (do not use when using VOIs as init) [default: -1]
   -S <subsets>, --subsets=<subsets>  number of subsets [default: 2]
   -i <subiterations>, --subiterations=<subiterations>     number of sub-iterations [default: 14]

SIRF_data_preparation/data_utilities.py

@@ -17,18 +17,25 @@
 # DATA_PATH = '/home/KrisThielemans/devel/PETRIC/data'
 this_directory = os.path.dirname(__file__)
 repo_directory = os.path.dirname(this_directory)
+ORG_DATA_PATH = os.path.join(repo_directory, 'orgdata')
 DATA_PATH = os.path.join(repo_directory, 'data')
 
 
-def the_data_path(*data_type):
+def the_data_path(*folders):
     '''
     Returns the path to data.
 
-    data_type: either 'PET', 'MR' or 'Synergistic', or use multiple arguments for
-    subdirectories like the_data_path('PET', 'mMR', 'NEMA_IQ').
+    data_type: subfolders like the_data_path('Siemens_mMR_ACR').
     '''
-    return os.path.join(DATA_PATH, *data_type)
+    return os.path.join(DATA_PATH, *folders)
 
+def the_orgdata_path(*folders):
+    '''
+    Returns the path to original data (for downloads/processing)
+
+    data_type: subfolders like the_orgdata_path('Siemens_mMR_ACR', 'processing').
+    '''
+    return os.path.join(ORG_DATA_PATH, *folders)
 
 def fix_siemens_norm_EOL(in_filename, out_filename):
     with open(in_filename, mode="rb") as f:

SIRF_data_preparation/dataset_settings.py

@@ -13,6 +13,9 @@ def get_settings(scanID: str):
     if scanID == 'Siemens_mMR_NEMA_IQ':
         slices = {'transverse_slice': 72, 'coronal_slice': 109} # , 'sagittal_slice': 89}
         num_subsets = 7
+    elif scanID == 'Siemens_mMR_ACR':
+        slices = { 'transverse_slice': 99}
+        num_subsets = 7
     elif scanID == 'NeuroLF_Hoffman_Dataset':
         slices = {'transverse_slice': 72}
         num_subsets = 16

SIRF_data_preparation/run_BSREM.py

@@ -4,13 +4,14 @@
   run_BSREM.py <data_set> [--help | options]
 
 Arguments:
-  <data_set>  path to data files as well as prefix to use
+  <data_set>  path to data files as well as prefix to use (e.g. Siemens_mMR_NEMA_EQ)
 
 """
 # Copyright 2024 Rutherford Appleton Laboratory STFC
 # Copyright 2024 University College London
 # Licence: Apache-2.0
-__version__ = '0.1.0'
+__version__ = '0.2.0'
+
 from pathlib import Path
 
 import matplotlib.pyplot as plt