From 7d26561c0c888853d1c556c8ce14831fc6600afd Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Wed, 21 Aug 2024 14:17:41 +0100 Subject: [PATCH 1/9] update Siemens_mMR_ACR README with additional steps Using Pawel's registered mu-map now. --- .../Siemens_mMR_ACR/README.md | 21 +++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/SIRF_data_preparation/Siemens_mMR_ACR/README.md b/SIRF_data_preparation/Siemens_mMR_ACR/README.md index e7c6771..249e836 100644 --- a/SIRF_data_preparation/Siemens_mMR_ACR/README.md +++ b/SIRF_data_preparation/Siemens_mMR_ACR/README.md @@ -4,9 +4,22 @@ 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) +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 data/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 + gunzip data/Siemens_mMR_ACR/output/acr-mumap-complete.nii.gz + stir_math data/Siemens_mMR_ACR/processing/reg_mumap.hv data/Siemens_mMR_ACR/output/acr-mumap-complete.nii + ``` + e. uncompress .gz VOI files as STIR isn't happy with them + ```sh + for f in data/Siemens_mMR_ACR/output/sampling_masks/*gz; do gunzip $f; done + ``` 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 +7. `python SIRF_data_preparation/create_initial_images.py -s 200 data/Siemens_mMR_ACR/final --template_image=data/Siemens_mMR_ACR/output/sampling_masks/acr-all-sampling-0-2mm_dipy.nii` +8. ` python ~/devel/PETRIC/SIRF_data_preparation/data_QC.py --srcdir='data/Siemens_mMR_ACR/final' --transverse_slice=99` From df947d9cc4e5a2dbe40736d67ba477bcc3b8e634 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Wed, 21 Aug 2024 14:18:46 +0100 Subject: [PATCH 2/9] initial BSREM for Siemens_mMR_ACR --- .../Siemens_mMR_ACR/BSREM.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) create mode 100644 SIRF_data_preparation/Siemens_mMR_ACR/BSREM.py diff --git a/SIRF_data_preparation/Siemens_mMR_ACR/BSREM.py b/SIRF_data_preparation/Siemens_mMR_ACR/BSREM.py new file mode 100644 index 0000000..10e2cda --- /dev/null +++ b/SIRF_data_preparation/Siemens_mMR_ACR/BSREM.py @@ -0,0 +1,19 @@ +from petric import MetricsWithTimeout, get_data +from sirf.contrib.BSREM.BSREM import BSREM1 +from sirf.contrib.partitioner import partitioner + +dataset = 'Siemens_mMR_ACR' +outdir = "./output/BSREM_" + dataset +data = get_data(srcdir="./data/" + dataset + "/final/", outdir=outdir) +num_subsets = 7 +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) + +algo = BSREM1(data_sub, obj_funs, initial=data.OSEM_image, initial_step_size=.3, relaxation_eta=.01, + update_objective_interval=10) +algo.run(15000, callbacks=[MetricsWithTimeout(seconds=3600*30, outdir=outdir, transverse_slice=99)]) From 194030208ad7138c1910dc8196b7103b3b4bd889 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Thu, 22 Aug 2024 09:46:49 +0100 Subject: [PATCH 3/9] clarify path of template_image --- SIRF_data_preparation/create_initial_images.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SIRF_data_preparation/create_initial_images.py b/SIRF_data_preparation/create_initial_images.py index 8af574c..aac29c0 100644 --- a/SIRF_data_preparation/create_initial_images.py +++ b/SIRF_data_preparation/create_initial_images.py @@ -7,7 +7,7 @@ path to data files Options: - -t , --template_image= filename of image to use for data sizes [default: VOI.hv] + -t , --template_image= filename (relative to ) of image to use for data sizes [default: PETRIC/VOI_whole_phantom.hv] -s , --xy-size= force xy-size (do not use when using VOIs as init) [default: -1] -S , --subsets= number of subsets [default: 2] -i , --subiterations= number of sub-iterations [default: 14] From cc83b458d0e657aa381b32e972c2067eddd061fc Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Thu, 22 Aug 2024 09:47:28 +0100 Subject: [PATCH 4/9] add the_orgdata_path for downloaded data etc --- SIRF_data_preparation/data_utilities.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/SIRF_data_preparation/data_utilities.py b/SIRF_data_preparation/data_utilities.py index 54b9a7c..817c6b1 100644 --- a/SIRF_data_preparation/data_utilities.py +++ b/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: From ac3531fe1f347808d5c7799c673d1e3b6f03111f Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Thu, 22 Aug 2024 09:56:30 +0100 Subject: [PATCH 5/9] add Siemens_mMR_ACR --- SIRF_data_preparation/dataset_settings.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/SIRF_data_preparation/dataset_settings.py b/SIRF_data_preparation/dataset_settings.py index 664a963..4fb0e55 100644 --- a/SIRF_data_preparation/dataset_settings.py +++ b/SIRF_data_preparation/dataset_settings.py @@ -11,6 +11,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 From 260a847c93ab63658ecaa10f05b4713b70336dd6 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Thu, 22 Aug 2024 09:58:19 +0100 Subject: [PATCH 6/9] Siemens_mMR_ACR processing updates - change to using orgdata for download/intermediate output - add hardware mu-map --- .gitignore | 1 + .../Siemens_mMR_ACR/README.md | 22 ++++++++++++------- .../Siemens_mMR_ACR/download.py | 10 +++------ .../Siemens_mMR_ACR/prepare.py | 14 +++++------- .../Siemens_mMR_ACR/register_mumap.py | 13 ++++------- SIRF_data_preparation/run_BSREM.py | 8 +++---- 6 files changed, 31 insertions(+), 37 deletions(-) diff --git a/.gitignore b/.gitignore index 0fb5791..48e71f3 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ data/ +orgdata/ output/ tmp*/ err*.txt diff --git a/SIRF_data_preparation/Siemens_mMR_ACR/README.md b/SIRF_data_preparation/Siemens_mMR_ACR/README.md index 249e836..4a90e54 100644 --- a/SIRF_data_preparation/Siemens_mMR_ACR/README.md +++ b/SIRF_data_preparation/Siemens_mMR_ACR/README.md @@ -6,20 +6,26 @@ Steps to follow: 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 and ideally would be renamed) -5. `python SIRF_data_preparation/Siemens_mMR_ACR/register_mumap.py` (output is data/Siemens_mMR_ACR/processing/reg_mumap.hv etc) +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 - gunzip data/Siemens_mMR_ACR/output/acr-mumap-complete.nii.gz - stir_math data/Siemens_mMR_ACR/processing/reg_mumap.hv data/Siemens_mMR_ACR/output/acr-mumap-complete.nii + 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 + e. uncompress .gz VOI files as STIR isn't happy with them (probably not necessary) ```sh - for f in data/Siemens_mMR_ACR/output/sampling_masks/*gz; do gunzip $f; done + for f in orgdata/Siemens_mMR_ACR/output/sampling_masks/*gz; do gunzip $f; done ``` -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 --template_image=data/Siemens_mMR_ACR/output/sampling_masks/acr-all-sampling-0-2mm_dipy.nii` -8. ` python ~/devel/PETRIC/SIRF_data_preparation/data_QC.py --srcdir='data/Siemens_mMR_ACR/final' --transverse_slice=99` +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. `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) +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 radionuclide and duration info which got lost. +11. `python SIRF_data_preparation/run_OSEM.py Siemens_mMR_ACR` +12. `python SIRF_data_preparation/run_BSREM.py Siemens_mMR_ACR` diff --git a/SIRF_data_preparation/Siemens_mMR_ACR/download.py b/SIRF_data_preparation/Siemens_mMR_ACR/download.py index 5569930..d073039 100644 --- a/SIRF_data_preparation/Siemens_mMR_ACR/download.py +++ b/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') diff --git a/SIRF_data_preparation/Siemens_mMR_ACR/prepare.py b/SIRF_data_preparation/Siemens_mMR_ACR/prepare.py index 3bd4207..7cce5c5 100644 --- a/SIRF_data_preparation/Siemens_mMR_ACR/prepare.py +++ b/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) diff --git a/SIRF_data_preparation/Siemens_mMR_ACR/register_mumap.py b/SIRF_data_preparation/Siemens_mMR_ACR/register_mumap.py index 02a1766..8735d9d 100644 --- a/SIRF_data_preparation/Siemens_mMR_ACR/register_mumap.py +++ b/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') diff --git a/SIRF_data_preparation/run_BSREM.py b/SIRF_data_preparation/run_BSREM.py index cc48211..b85938a 100644 --- a/SIRF_data_preparation/run_BSREM.py +++ b/SIRF_data_preparation/run_BSREM.py @@ -4,13 +4,13 @@ run_BSREM.py [--help | options] Arguments: - path to data files as well as prefix to use + 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' import logging import math @@ -38,7 +38,7 @@ SRCDIR = PETRICDIR / 'data' OUTDIR = PETRICDIR / 'output' -outdir = OUTDIR/ scanID / "BSREM" +outdir = OUTDIR / scanID / "BSREM" srcdir = SRCDIR / scanID #log.info("Finding files in %s", srcdir) @@ -57,7 +57,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(30000, callbacks=[MetricsWithTimeout(**settings.slices, outdir=outdir, seconds=3600*100)]) #%% fig = plt.figure() data_QC.plot_image(algo.get_output(), **settings.slices) From c4819abf6eff2f514b0b52bd87988363677be50e Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Thu, 22 Aug 2024 12:01:55 +0100 Subject: [PATCH 7/9] remove remnant VOI_checks (now in data_QC) --- .../Siemens_mMR_NEMA_VOIs.py | 46 +------------------ 1 file changed, 1 insertion(+), 45 deletions(-) diff --git a/SIRF_data_preparation/Siemens_mMR_NEMA_IQ/Siemens_mMR_NEMA_VOIs.py b/SIRF_data_preparation/Siemens_mMR_NEMA_IQ/Siemens_mMR_NEMA_VOIs.py index 54cb49b..a2f5611 100644 --- a/SIRF_data_preparation/Siemens_mMR_NEMA_IQ/Siemens_mMR_NEMA_VOIs.py +++ b/SIRF_data_preparation/Siemens_mMR_NEMA_IQ/Siemens_mMR_NEMA_VOIs.py @@ -78,59 +78,15 @@ plt.figure() data_QC.plot_image(reference_image, **slices, alpha=allVOIs) #%% -from data_QC import plot_image_if_exists +from data_QC import plot_image_if_exists, VOI_checks from data_QC import plot_image from data_QC import VOI_mean from scipy import ndimage -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] From 92d40a724f6da17f9b5ffaf7552a9c96a9857dcc Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Thu, 22 Aug 2024 12:02:28 +0100 Subject: [PATCH 8/9] add VOIs for the mMR ACR --- .../Siemens_mMR_ACR/README.md | 9 +- .../Siemens_mMR_ACR/VOI_prep.py | 91 +++++++++++++++++++ 2 files changed, 96 insertions(+), 4 deletions(-) create mode 100644 SIRF_data_preparation/Siemens_mMR_ACR/VOI_prep.py diff --git a/SIRF_data_preparation/Siemens_mMR_ACR/README.md b/SIRF_data_preparation/Siemens_mMR_ACR/README.md index 4a90e54..8d1562d 100644 --- a/SIRF_data_preparation/Siemens_mMR_ACR/README.md +++ b/SIRF_data_preparation/Siemens_mMR_ACR/README.md @@ -23,9 +23,10 @@ Steps to follow: ```sh stir_math --accumulate orgdata/Siemens_mMR_ACR/processing/reg_mumap.hv orgdata/Siemens_mMR_ACR/output/ACR_hardware-to-STIR.nii.gz ``` -7. `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. 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 radionuclide and duration info which got lost. -11. `python SIRF_data_preparation/run_OSEM.py Siemens_mMR_ACR` -12. `python SIRF_data_preparation/run_BSREM.py Siemens_mMR_ACR` +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` diff --git a/SIRF_data_preparation/Siemens_mMR_ACR/VOI_prep.py b/SIRF_data_preparation/Siemens_mMR_ACR/VOI_prep.py new file mode 100644 index 0000000..a712b27 --- /dev/null +++ b/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) From cba1219ec6a7196d126388c9f7bbfe2506a73de5 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Thu, 22 Aug 2024 12:25:29 +0100 Subject: [PATCH 9/9] removed obsolete file --- .../Siemens_mMR_ACR/BSREM.py | 19 ------------------- 1 file changed, 19 deletions(-) delete mode 100644 SIRF_data_preparation/Siemens_mMR_ACR/BSREM.py diff --git a/SIRF_data_preparation/Siemens_mMR_ACR/BSREM.py b/SIRF_data_preparation/Siemens_mMR_ACR/BSREM.py deleted file mode 100644 index 10e2cda..0000000 --- a/SIRF_data_preparation/Siemens_mMR_ACR/BSREM.py +++ /dev/null @@ -1,19 +0,0 @@ -from petric import MetricsWithTimeout, get_data -from sirf.contrib.BSREM.BSREM import BSREM1 -from sirf.contrib.partitioner import partitioner - -dataset = 'Siemens_mMR_ACR' -outdir = "./output/BSREM_" + dataset -data = get_data(srcdir="./data/" + dataset + "/final/", outdir=outdir) -num_subsets = 7 -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) - -algo = BSREM1(data_sub, obj_funs, initial=data.OSEM_image, initial_step_size=.3, relaxation_eta=.01, - update_objective_interval=10) -algo.run(15000, callbacks=[MetricsWithTimeout(seconds=3600*30, outdir=outdir, transverse_slice=99)])