From e3aaf993fd90245895cbca991621933363f65d01 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Sun, 20 Oct 2024 21:35:26 +0000 Subject: [PATCH 1/8] added GE_D690_NEMA_IQ --- SIRF_data_preparation/dataset_settings.py | 2 +- petric.py | 10 ++++++---- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/SIRF_data_preparation/dataset_settings.py b/SIRF_data_preparation/dataset_settings.py index dc17ab6..a21b9fe 100644 --- a/SIRF_data_preparation/dataset_settings.py +++ b/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, 'Siemens_Vision600_ZrNEMAIQ': 5} + 'NeuroLF_Esser_Dataset': 8, 'Siemens_Vision600_ZrNEMAIQ': 5, 'GE_D690_NEMA_IQ': 16} @dataclass diff --git a/petric.py b/petric.py index b269fe8..4c26bf7 100755 --- a/petric.py +++ b/petric.py @@ -292,9 +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': {}, 'Siemens_Vision600_ZrNEMAIQ': {'transverse_slice': 60}} + '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}, 'GE_D690_NEMA_IQ': {'transverse_slice': 23}} if SRCDIR.is_dir(): # create list of existing data @@ -319,7 +319,9 @@ def get_image(fname): (SRCDIR / "NeuroLF_Esser_Dataset", OUTDIR / "NeuroLF_Esser", [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'])])] + [MetricsWithTimeout(outdir=OUTDIR / "Vision600_ZrNEMAIQ", **DATA_SLICES['Siemens_Vision600_ZrNEMAIQ'])]), + (SRCDIR / "GE_D690_NEMA_IQ", OUTDIR / "D690_NEMA", + [MetricsWithTimeout(outdir=OUTDIR / "D690_NEMA", **DATA_SLICES['GE_D690_NEMA_IQ'])])] else: log.warning("Source directory does not exist: %s", SRCDIR) data_dirs_metrics = [(None, None, [])] # type: ignore From b8e246f2ef2a0ba7575f2da3bfc7caf7eca4874a Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Sun, 20 Oct 2024 21:38:16 +0000 Subject: [PATCH 2/8] minor fix --- SIRF_data_preparation/create_Hoffman_VOIs.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/SIRF_data_preparation/create_Hoffman_VOIs.py b/SIRF_data_preparation/create_Hoffman_VOIs.py index 9d6da62..799d7fc 100644 --- a/SIRF_data_preparation/create_Hoffman_VOIs.py +++ b/SIRF_data_preparation/create_Hoffman_VOIs.py @@ -35,7 +35,7 @@ from SIRF_data_preparation.registration_utilities import STIR_to_nii, STIR_to_nii_hv, register_it, resample_STIR # %% -__version__ = "0.1.0" +__version__ = "0.1.1" write_PETRIC_VOIs = True if "ipykernel" not in sys.argv[0]: # clunky way to be able to set variables from within jupyter/VScode without docopt @@ -44,6 +44,7 @@ # logging.basicConfig(level=logging.INFO) scanID = args["--dataset"] + srcdir = args['--srcdir'] if scanID is None: print("Need to set the --dataset argument") exit(1) @@ -52,10 +53,9 @@ else: # set it by hand, e.g. scanID = "NeuroLF_Hoffman_Dataset" + srcdir = None write_PETRIC_VOIs = False -srcdir = args['--srcdir'] - # %% standard PETRIC directories if srcdir is None: srcdir = the_data_path(scanID) From dc04f67e060fd7cde606f4c092a483df5b8afdd4 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Sun, 20 Oct 2024 21:40:59 +0000 Subject: [PATCH 3/8] add function to create NEMA_IQ VOIs building on previous attempts for the ZrNEMAIQ, but hopefully more generic. Needs https://github.com/SyneRBI/SIRF-Contribs/pull/22 --- SIRF_data_preparation/create_NEMA_IQ_VOIs.py | 165 +++++++++++++++++++ 1 file changed, 165 insertions(+) create mode 100644 SIRF_data_preparation/create_NEMA_IQ_VOIs.py diff --git a/SIRF_data_preparation/create_NEMA_IQ_VOIs.py b/SIRF_data_preparation/create_NEMA_IQ_VOIs.py new file mode 100644 index 0000000..a496887 --- /dev/null +++ b/SIRF_data_preparation/create_NEMA_IQ_VOIs.py @@ -0,0 +1,165 @@ +#!/usr/bin/env python +"""Create VOIs for a scan of the NEMA IQ phantom + +Usage: + create_NEMA_IQ_VOIs.py [--help | options] + + +Options: + -h, --help + --dataset= dataset name (required) + --srcdir= pathname. Will default to current directory unless dataset is set + --central_VOI= Use the normal background VOI (in the centre of the phantom) or a + VOI "below" the largest sphere [default: True] + --angle_smallest_sphere= angle (in degrees) of smallest sphere [default: 210] + --spheres= string to be parsed as tuple of sphere numbers (1 is smallest) [default: (1,3,5)] +""" + +# Copyright 2024 University College London +# Licence: Apache-2.0 + +# %% +# %load_ext autoreload +# %autoreload 2 +# %% imports +import ast +import os +import sys +import typing +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np +import numpy.typing as npt +import scipy.ndimage as ndimage +from docopt import docopt + +import sirf.contrib.NEMA.generate_nema_rois as generate_nema_rois +import sirf.STIR as STIR +from SIRF_data_preparation.data_QC import VOI_checks, plot_image +from SIRF_data_preparation.data_utilities import the_data_path, the_orgdata_path +from SIRF_data_preparation.dataset_settings import get_settings + +# %% +__version__ = "0.1.0" + +if "ipykernel" not in sys.argv[0]: # clunky way to be able to set variables from within jupyter/VScode without docopt + args = docopt(__doc__, argv=None, version=__version__) + + # logging.basicConfig(level=logging.INFO) + + scanID = args["--dataset"] + srcdir = args['--srcdir'] + if scanID is None: + print("Need to set the --dataset argument") + exit(1) + angle_smallest_sphere = args["--angle_smallest_sphere"] + central_VOI = args["--central_VOI"] + spheres = ast.literal_eval(args["--spheres"]) +else: + # set it by hand, e.g. + scanID = "GE_D690_NEMA_IQ" + srcdir = None + central_VOI = False + spheres = (2, 3, 5) + angle_smallest_sphere = 30 + +# %% standard PETRIC directories +if srcdir is None: + srcdir = the_data_path(scanID) +srcdir = Path(srcdir) +intermediate_data_path = Path(the_orgdata_path(scanID, "processing")) +os.makedirs(intermediate_data_path, exist_ok=True) +settings = get_settings(scanID) +slices = settings.slices + +print("srcdir:", srcdir) +print("processingdir:", intermediate_data_path) +print("angle_smallest_sphere:", angle_smallest_sphere) +print("central_VOI:", central_VOI) +print("spheres:", spheres) + + +# %% Function to find the n-th connected component (counting by size) +def connected_component(arr: npt.NDArray[bool], order=0) -> npt.NDArray[bool]: + """Return connected component of a given order in an array. order=0 returns the largest.""" + labels, num_labels = ndimage.label(arr) + num_elems = [(labels == i).sum() for i in range(num_labels)] + idx = np.argsort(num_elems)[-1 - order] + return labels == idx + + +def read_unregistered_VOIs(intermediate_datadir: Path) -> typing.List[STIR.ImageData]: + orgVOIs = [] + for i in range(1, 8): + orgVOIs.append(STIR.ImageData(str(intermediate_datadir / f"unregistered_sphere{i}.nii"))) + return orgVOIs + + +def read_orgVOIs(intermediate_datadir: Path) -> typing.List[STIR.ImageData]: + orgVOIs = [] + for i in range(1, 8): + orgVOIs.append(STIR.ImageData(str(intermediate_datadir / f"S{i}.nii"))) + return orgVOIs + + +def create_whole_object_mask(reference_image: STIR.ImageData) -> STIR.ImageData: + whole_object_mask = reference_image.clone() + ref_arr = reference_image.as_array() + whole_object_mask.fill(ndimage.binary_erosion(ref_arr > np.percentile(ref_arr, 99) / 20, iterations=2)) + # plot_image(whole_object_mask, **slices) + return whole_object_mask + + +def create_background_VOI(orgVOIs: typing.List[STIR.ImageData], central_VOI: bool) -> STIR.ImageData: + if central_VOI: + return orgVOIs[6] + # 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() + # COM = ndimage.center_of_mass(background_arr) + z_voxel_size = background_mask.voxel_sizes()[0] + # 37 mm diameter + z_shift = int(np.floor(37 * 1.2 / z_voxel_size)) + background_arr = np.concatenate(( + background_arr[0:z_shift, :, :] * 0, + background_arr[0:-z_shift, :, :], + )) + background_mask.fill(background_arr) + return background_mask + + +# %% +reference_image = STIR.ImageData(str(srcdir / 'OSEM_image.hv')) +plot_image(reference_image, **slices) +# %% +thresholded_image = reference_image.clone() +thresholded_arr = thresholded_image.as_array() +thresholded_arr[thresholded_arr < np.max(thresholded_arr) / 5] = 0 +thresholded_image.fill(thresholded_arr) +# %% +# currently needs trailing slash +generate_nema_rois.data_output_path = str(intermediate_data_path) + "/" +generate_nema_rois.generate_nema_rois(thresholded_image, angle_smallest=angle_smallest_sphere) +# %% +unregVOIs = read_unregistered_VOIs(intermediate_data_path) +orgVOIs = read_orgVOIs(intermediate_data_path) +# %% +whole_object_mask = create_whole_object_mask(reference_image) +background_mask = create_background_VOI(orgVOIs, central_VOI) +sphere_VOIs = [orgVOIs[i - 1] for i in spheres] +# %% +VOIs = [whole_object_mask, background_mask] + sphere_VOIs +VOInames = ["VOI_whole_object", "VOI_background"] + [f"VOI_sphere{s}" for s in spheres] + +# %% write PETRIC VOIs +datadir = Path(srcdir) / "PETRIC" +print(f"Writing VOIs to {datadir}") +os.makedirs(datadir, exist_ok=True) +for VOI, n in zip(VOIs, VOInames): + VOI.write(str(datadir / n)) +# %% +VOI_checks(VOInames, OSEM_image=reference_image, reference_image=None, srcdir=datadir) +# %% +plt.show() From 2a72ef16876e98f9bf5928a375bc51e9c9403915 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Sun, 20 Oct 2024 21:56:41 +0000 Subject: [PATCH 4/8] fix location of final .png --- SIRF_data_preparation/data_QC.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SIRF_data_preparation/data_QC.py b/SIRF_data_preparation/data_QC.py index 962f1fe..c12e526 100644 --- a/SIRF_data_preparation/data_QC.py +++ b/SIRF_data_preparation/data_QC.py @@ -167,7 +167,7 @@ def VOI_checks(allVOInames, OSEM_image=None, reference_image=None, srcdir='.', * if OSEM_image is not None: plt.figure() - plot_image(OSEM_image, alpha=allVOIs, save_name="OSEM_image_and_VOIs", **kwargs) + plot_image(OSEM_image, alpha=allVOIs, save_name=prefix + "OSEM_image_and_VOIs", **kwargs) # unformatted print of VOI values for now print(allVOInames) From cced554440b452fb2aea138f75822a97f09cd59f Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Sun, 20 Oct 2024 22:29:33 +0000 Subject: [PATCH 5/8] minor changes for plotting --- SIRF_data_preparation/plot_BSREM_metrics.py | 28 ++++++++++----------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/SIRF_data_preparation/plot_BSREM_metrics.py b/SIRF_data_preparation/plot_BSREM_metrics.py index 6a0ca06..e039218 100644 --- a/SIRF_data_preparation/plot_BSREM_metrics.py +++ b/SIRF_data_preparation/plot_BSREM_metrics.py @@ -26,29 +26,37 @@ # scanID = 'NeuroLF_Hoffman_Dataset' # scanID = 'Siemens_mMR_NEMA_IQ' # scanID = 'Mediso_NEMA_IQ' -scanID = 'Siemens_Vision600_Hoffman' +# scanID = 'Siemens_Vision600_Hoffman' # scanID = 'NeuroLF_Esser_Dataset' # scanID = 'Siemens_Vision600_ZrNEMAIQ' +scanID = 'GE_D690_NEMA_IQ' srcdir = the_data_path(scanID) outdir = OUTDIR / scanID OSEMdir = outdir / 'OSEM' datadir = outdir / 'BSREM' # we will check for images obtained after restarting BSREM (with new settings) datadir1 = outdir / 'BSREM_cont' - +# datadir = Path('/opt/runner/logs/Casper/BSREM/') / scanID OSEM_image = STIR.ImageData(str(datadir / 'iter_0000.hv')) settings = get_settings(scanID) slices = settings.slices -cmax = numpy.percentile(OSEM_image.as_array(), 99.995) +cmax = numpy.percentile(OSEM_image.as_array(), 99.9) +# %% +data = get_data(srcdir=srcdir, outdir=None, read_sinos=False) # %% -image = data_QC.plot_image_if_exists(str(datadir / 'iter_final'), **slices, vmax=cmax) +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')) # image2=STIR.ImageData(datadir+'iter_14000.hv') # %% image2 = OSEM_image -diff = image2 - image -print("relative l1-norm diff final-OSEM:", diff.abs().max() / image.max()) +diff = image2 - reference_image +print("relative l1-norm diff final-OSEM:", diff.abs().max() / reference_image.max()) data_QC.plot_image(diff, **slices, vmin=-cmax / 100, vmax=cmax / 100) # %% objs = read_objectives(datadir) @@ -69,14 +77,6 @@ fig.savefig(outdir / f'{scanID}_BSREM_objectives_last.png') # %% -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')) 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 From f3490e99dc10c47b1dff78f7931b50baf56eecc3 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Sun, 20 Oct 2024 22:54:42 +0000 Subject: [PATCH 6/8] allow setting of output dir --- SIRF_data_preparation/run_BSREM.py | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/SIRF_data_preparation/run_BSREM.py b/SIRF_data_preparation/run_BSREM.py index d962ebe..96a7c39 100644 --- a/SIRF_data_preparation/run_BSREM.py +++ b/SIRF_data_preparation/run_BSREM.py @@ -1,4 +1,6 @@ -"""Run BSREM for PETRIC +#!/usr/bin/env python +""" +Run BSREM for PETRIC Usage: run_BSREM.py [--help | options] @@ -9,17 +11,18 @@ Options: --updates= number of updates to run [default: 15000] --initial_image= optional initial image, normally the OSEM_image from get_data. - If specified, - output will be in BSREM_cont. + --num_subsets= number of subsets. If not specified, will use dataset_settings.get_settings. --initial_step_size= start stepsize [default: .3] --relaxation_eta= relaxation factor per epoch [default: .01] --interval= interval to save [default: 80] + --outreldir= optional relative path to override + (defaults to 'BSREM' or 'BSREM_cont' if initial_image is set) """ # Copyright 2024 Rutherford Appleton Laboratory STFC # Copyright 2024 University College London # Licence: Apache-2.0 -__version__ = '0.3.0' +__version__ = '0.4.0' from pathlib import Path @@ -44,6 +47,7 @@ initial_step_size = float(args['--initial_step_size']) relaxation_eta = float(args['--relaxation_eta']) interval = int(args['--interval']) +outreldir = args['--outreldir'] if not all((SRCDIR.is_dir(), OUTDIR.is_dir())): PETRICDIR = Path('~/devel/PETRIC').expanduser() @@ -61,11 +65,11 @@ if initial_image is None: initial_image_name = "OSEM" initial_image = data.OSEM_image - outdir = outdir / "BSREM" + outdir = outdir / ("BSREM" if outreldir is None else outreldir) else: initial_image_name = initial_image initial_image = STIR.ImageData(initial_image) - outdir = outdir / "BSREM_cont" + outdir = outdir / ("BSREM_cont" if outreldir is None else outreldir) print("Penalisation factor:", data.prior.get_penalisation_factor()) print("num_subsets:", num_subsets) From 75951d938a4665310b05dc6092fbb00ed199d094 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Sun, 20 Oct 2024 23:52:55 +0000 Subject: [PATCH 7/8] add steps for GE_D690_NEMA_IQ --- SIRF_data_preparation/GE_D690_NEMA_IQ/README.md | 11 +++++++++++ 1 file changed, 11 insertions(+) create mode 100644 SIRF_data_preparation/GE_D690_NEMA_IQ/README.md diff --git a/SIRF_data_preparation/GE_D690_NEMA_IQ/README.md b/SIRF_data_preparation/GE_D690_NEMA_IQ/README.md new file mode 100644 index 0000000..9e5c622 --- /dev/null +++ b/SIRF_data_preparation/GE_D690_NEMA_IQ/README.md @@ -0,0 +1,11 @@ +# Preparation steps for NEMA IQ data scanned on GE Discovery 690 +## Reference solution +```sh +python -m SIRF_data_preparation.run_BSREM GE_D690_NEMA_IQ --updates=500 --initial_step_size=1 --interval=10 --outreldir=BSREM1 +python -m SIRF_data_preparation.run_BSREM GE_D690_NEMA_IQ --updates=7520 --initial_step_size=.8 --relaxation_eta=.01 --num_subsets=6 --initial_image=output/GE_D690_NEMA_IQ/BSREM1/iter_final.hv --outreldir=BSREM2 +python -m SIRF_data_preparation.run_BSREM GE_D690_NEMA_IQ --updates=1000 --initial_step_size=.2 --relaxation_eta=.01 --num_subsets=1 --initial_image=output/GE_D690_NEMA_IQ/BSREM2/iter_7520.hv --outreldir=BSREM3 --interval=2 +``` +## VOIs +```sh +python -m SIRF_data_preparation.create_NEMA_IQ_VOIs --dataset=GE_D690_NEMA_IQ --central_VOI=False --angle_smallest_sphere=30 --spheres="(2,3,5)" +``` \ No newline at end of file From 4d65710e4f007dd40112aeff754657217069dbb8 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sun, 20 Oct 2024 23:53:23 +0000 Subject: [PATCH 8/8] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- SIRF_data_preparation/GE_D690_NEMA_IQ/README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SIRF_data_preparation/GE_D690_NEMA_IQ/README.md b/SIRF_data_preparation/GE_D690_NEMA_IQ/README.md index 9e5c622..487033a 100644 --- a/SIRF_data_preparation/GE_D690_NEMA_IQ/README.md +++ b/SIRF_data_preparation/GE_D690_NEMA_IQ/README.md @@ -8,4 +8,4 @@ python -m SIRF_data_preparation.run_BSREM GE_D690_NEMA_IQ --updates=1000 --initi ## VOIs ```sh python -m SIRF_data_preparation.create_NEMA_IQ_VOIs --dataset=GE_D690_NEMA_IQ --central_VOI=False --angle_smallest_sphere=30 --spheres="(2,3,5)" -``` \ No newline at end of file +```