diff --git a/.gitignore b/.gitignore index e71f246..69eeb8d 100644 --- a/.gitignore +++ b/.gitignore @@ -16,4 +16,4 @@ __pycache__/ *.v *.csv # Tensorboard files -events.* \ No newline at end of file +events.* diff --git a/SIRF_data_preparation/GE_DMI3_Torso/README.md b/SIRF_data_preparation/GE_DMI3_Torso/README.md new file mode 100644 index 0000000..f0945d5 --- /dev/null +++ b/SIRF_data_preparation/GE_DMI3_Torso/README.md @@ -0,0 +1,25 @@ +# Torso phantom scanned on a GE Discovery MI (3 ring) + +This data is from a scan performed by GE of the +Data Spectrum [Antropomorphic Torso Phantom](https://www.spect.com/pdf/anthropomorphic-torso-phantom.pdf) +with additional lesion inserts of diameter 10mm. + +Activity ratio between "liver" and "torso" is approximately 2. Lungs +are empty (air). Contrast of the lesions is unfortunately unknown. + +"Corrections" were obtained with GE Duetto and converted to STIR Interfile +using internal scripts. + +Due to data size issues and count levels, the supplied data is **reduced** from the original scan +by keeping only the "direct" sinograms (segment 0, corresponding to min/max ring differences -1/+1), e.g. on the command line +``` +SSRB prompts.hs prompts_f1b1.hs 1 1 0 0 1 +``` + +In addition, as current PETRIC examples do not cope with non-TOF multfactors for TOF data, the +multfactors were expanded to TOF (by duplication). + +## VOIs + +VOIs were determined by Kris Thielemans on reconstructed images and approximately placed +where some of the lesions are visible. They are not necessarily at the correct location. diff --git a/SIRF_data_preparation/GE_DMI3_Torso/VOI_prep.py b/SIRF_data_preparation/GE_DMI3_Torso/VOI_prep.py new file mode 100644 index 0000000..2e3d31f --- /dev/null +++ b/SIRF_data_preparation/GE_DMI3_Torso/VOI_prep.py @@ -0,0 +1,83 @@ +# %% file to prepare VOIs for the Torso phantom, using coordinates derived manually +import os +from pathlib import Path + +import numpy as np +from scipy.ndimage import binary_fill_holes + +import sirf.STIR as STIR +from SIRF_data_preparation.data_utilities import the_data_path + +# %% +scanID = 'GE_DMI3_Torso' +data_path = Path(the_data_path(scanID)) +output_path = Path(the_data_path(scanID, 'PETRIC')) +os.makedirs(output_path, exist_ok=True) +# %% +image_filename = str(data_path / 'OSEM_image.hv') +image = STIR.ImageData(image_filename) +# %% find whole phantom VOI by thresholding/filling +threshold = image.max() * .05 +im_arr = image.as_array() >= threshold +# we want to use binary_fill_holes to fill the spine and lung. However, +# that function does not fill holes connected to the boundary, so let's +# fill the boundary planes first +im_arr[0, :, :] = True +im_arr[-1, :, :] = True +im_arr = binary_fill_holes(im_arr) +# now exclude boundary +im_arr[0, :, :] = False +im_arr[-1, :, :] = False +VOI = image.allocate(0) +VOI.fill(im_arr.astype(np.float32)) +VOI.show() +VOI.write(str(output_path / 'VOI_whole_object.hv')) +# %% find centre image as KT got the coordinates from Amide +# note that geo.get_index_to_physical_point_matrix() is still unstable, so we won't use it. +# Warning: the following will likely fail for future versions of SIRF +zcentre = (image.shape[0] - 1) / 2 * image.voxel_sizes()[0] +# %% background VOI (in approx middle of torso) +shape = STIR.EllipticCylinder() +shape.set_length(23) +shape.set_radii((120, 120)) +shape.set_origin((zcentre + 0, 0, 0)) +VOI.fill(0) +VOI.add_shape(shape, scale=1) +VOI.write(str(output_path / 'VOI_background.hv')) +# %% spine (approximately top half of it) +shape = STIR.EllipticCylinder() +shape.set_length(100) +shape.set_radii((22, 22)) +shape.set_origin((zcentre + 42.6, 86.9, 3)) +VOI.fill(0) +VOI.add_shape(shape, scale=1) +VOI.write(str(output_path / 'VOI_spine.hv')) +VOI.show() +# %% lung lesion +shape = STIR.Ellipsoid() +shape.set_radius_z(5) +shape.set_radius_y(5) +shape.set_radius_x(5) +shape.set_origin((zcentre + 25, -42.2, -105.3)) +VOI.fill(0) +VOI.add_shape(shape, scale=1) +VOI.write(str(output_path / 'VOI_lung_lesion.hv')) +# %% +# %% liver lesion +shape = STIR.Ellipsoid() +shape.set_radius_z(5) +shape.set_radius_y(5) +shape.set_radius_x(5) +shape.set_origin((zcentre - 38.2, -38.9, -49.7)) +VOI.fill(0) +VOI.add_shape(shape, scale=1) +VOI.write(str(output_path / 'VOI_liver_lesion.hv')) +# %% lung lesion +shape = STIR.Ellipsoid() +shape.set_radius_z(5) +shape.set_radius_y(5) +shape.set_radius_x(5) +shape.set_origin((zcentre + 47.6, 12.5, 151.5)) +VOI.fill(0) +VOI.add_shape(shape, scale=1) +VOI.write(str(output_path / 'VOI_chest_lesion.hv')) diff --git a/SIRF_data_preparation/README.md b/SIRF_data_preparation/README.md index 1666620..e6b0ab1 100644 --- a/SIRF_data_preparation/README.md +++ b/SIRF_data_preparation/README.md @@ -27,7 +27,7 @@ for the Siemens mMR NEMA IQ data (on Zenodo): - `BSREM_*.py`: functions with specific settings for a particular data-set # Steps to follow to prepare data -If starting from Siemens list-mode data and letting SIRF take care of scatter etc, check for instance [Siemens_mMR_ACR/README.md](steps for Siemens mMR ACR). If pre-prepared data are given, check that naming of all files is correct. KT normally puts all data +If starting from Siemens mMR list-mode data and letting SIRF take care of scatter etc, check for instance [steps for Siemens mMR ACR](Siemens_mMR_ACR/README.md). If pre-prepared data are given, check that naming of all files is correct. KT normally puts all data in `~/devel/PETRIC/data/` with `datasetname` following convention of `scanner_phantom-name` as others (instructions below and indeed some scripts might assume this location). Change working directory to where data sits and add PETRIC to your python-path, e.g. ``` PYTHONPATH=~/devel/PETRIC:$PYTHONPATH` diff --git a/SIRF_data_preparation/Siemens_mMR_NEMA_IQ/prepare_mMR_NEMA_IQ_data.py b/SIRF_data_preparation/Siemens_mMR_NEMA_IQ/prepare_mMR_NEMA_IQ_data.py index c88ec51..b7af0c5 100644 --- a/SIRF_data_preparation/Siemens_mMR_NEMA_IQ/prepare_mMR_NEMA_IQ_data.py +++ b/SIRF_data_preparation/Siemens_mMR_NEMA_IQ/prepare_mMR_NEMA_IQ_data.py @@ -2,7 +2,7 @@ import logging import os -from ..data_utilities import prepare_challenge_Siemens_data, the_data_path, the_orgdata_path +from ..data_utilities import prepare_challenge_Siemens_data, the_data_path this_directory = os.path.dirname(__file__) repo_directory = os.path.dirname(this_directory) diff --git a/SIRF_data_preparation/Siemens_mMR_NEMA_IQ_lowcounts/README.md b/SIRF_data_preparation/Siemens_mMR_NEMA_IQ_lowcounts/README.md index 8acaf9b..814052a 100644 --- a/SIRF_data_preparation/Siemens_mMR_NEMA_IQ_lowcounts/README.md +++ b/SIRF_data_preparation/Siemens_mMR_NEMA_IQ_lowcounts/README.md @@ -5,4 +5,4 @@ This relies on the "normal" count level processsing to be done first, specifical Steps to follow: 1. `prepare_mMR_NEMA_IQ_data.py` 2. copy VOIs -3. further steps as [normal](../README.md#steps-to-follow-to-prepare-data) \ No newline at end of file +3. further steps as [normal](../README.md#steps-to-follow-to-prepare-data) diff --git a/SIRF_data_preparation/data_QC.py b/SIRF_data_preparation/data_QC.py index 9176da6..52d801c 100644 --- a/SIRF_data_preparation/data_QC.py +++ b/SIRF_data_preparation/data_QC.py @@ -30,8 +30,8 @@ from scipy import ndimage import sirf.STIR as STIR -from SIRF_data_preparation.dataset_settings import get_settings from SIRF_data_preparation.data_utilities import the_data_path +from SIRF_data_preparation.dataset_settings import get_settings STIR.AcquisitionData.set_storage_scheme('memory') diff --git a/SIRF_data_preparation/dataset_settings.py b/SIRF_data_preparation/dataset_settings.py index a0a9628..1f8f827 100644 --- a/SIRF_data_preparation/dataset_settings.py +++ b/SIRF_data_preparation/dataset_settings.py @@ -5,7 +5,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} + 'Mediso_NEMA_IQ': 12, 'Siemens_Vision600_thorax': 5, 'GE_DMI3_Torso': 8} @dataclass diff --git a/petric.py b/petric.py index 005d7dd..31e2d20 100755 --- a/petric.py +++ b/petric.py @@ -43,7 +43,7 @@ class Callback(cil_callbacks.Callback): CIL Callback but with `self.skip_iteration` checking `min(self.interval, algo.update_objective_interval)`. TODO: backport this class to CIL. """ - def __init__(self, interval: int = 1 << 31, **kwargs): + def __init__(self, interval: int = 3, **kwargs): super().__init__(**kwargs) self.interval = interval @@ -110,7 +110,7 @@ def __call__(self, algo: Algorithm): class QualityMetrics(ImageQualityCallback, Callback): """From https://github.com/SyneRBI/PETRIC/wiki#metrics-and-thresholds""" - def __init__(self, reference_image, whole_object_mask, background_mask, interval: int = 1 << 31, **kwargs): + def __init__(self, reference_image, whole_object_mask, background_mask, interval: int = 3, **kwargs): # TODO: drop multiple inheritance once `interval` included in CIL Callback.__init__(self, interval=interval) ImageQualityCallback.__init__(self, reference_image, **kwargs) @@ -205,6 +205,7 @@ class Dataset: whole_object_mask: STIR.ImageData | None background_mask: STIR.ImageData | None voi_masks: dict[str, STIR.ImageData] + FOV_mask: STIR.ImageData path: PurePath @@ -223,6 +224,10 @@ def get_data(srcdir=".", outdir=OUTDIR, sirf_verbosity=0): additive_term = STIR.AcquisitionData(str(srcdir / 'additive_term.hs')) mult_factors = STIR.AcquisitionData(str(srcdir / 'mult_factors.hs')) OSEM_image = STIR.ImageData(str(srcdir / 'OSEM_image.hv')) + # Find FOV mask + # WARNING: we are currently using Parralelproj with default settings, which uses a cylindrical FOV. + # The current code gives identical results to thresholding the sensitivity image (for those settings) + FOV_mask = STIR.TruncateToCylinderProcessor().process(OSEM_image.allocate(1)) kappa = STIR.ImageData(str(srcdir / 'kappa.hv')) if (penalty_strength_file := (srcdir / 'penalisation_factor.txt')).is_file(): penalty_strength = float(np.loadtxt(penalty_strength_file)) @@ -243,16 +248,15 @@ def get_image(fname): for voi in (srcdir / 'PETRIC').glob("VOI_*.hv") if voi.stem[4:] not in ('background', 'whole_object')} return Dataset(acquired_data, additive_term, mult_factors, OSEM_image, prior, kappa, reference_image, - whole_object_mask, background_mask, voi_masks, srcdir.resolve()) + whole_object_mask, background_mask, voi_masks, FOV_mask, srcdir.resolve()) DATA_SLICES = { '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': {}} + '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': {}} if SRCDIR.is_dir(): # create list of existing data @@ -260,10 +264,18 @@ def get_image(fname): data_dirs_metrics = [ (SRCDIR / "Siemens_mMR_NEMA_IQ", OUTDIR / "mMR_NEMA", [MetricsWithTimeout(outdir=OUTDIR / "mMR_NEMA", **DATA_SLICES['Siemens_mMR_NEMA_IQ'])]), + (SRCDIR / "Siemens_mMR_NEMA_IQ_lowcounts", OUTDIR / "mMR_NEMA_lowcounts", + [MetricsWithTimeout(outdir=OUTDIR / "mMR_NEMA_lowcounts", **DATA_SLICES['Siemens_mMR_NEMA_IQ_lowcounts'])]), (SRCDIR / "NeuroLF_Hoffman_Dataset", OUTDIR / "NeuroLF_Hoffman", [MetricsWithTimeout(outdir=OUTDIR / "NeuroLF_Hoffman", **DATA_SLICES['NeuroLF_Hoffman_Dataset'])]), (SRCDIR / "Siemens_Vision600_thorax", OUTDIR / "Vision600_thorax", - [MetricsWithTimeout(outdir=OUTDIR / "Vision600_thorax", **DATA_SLICES['Siemens_Vision600_thorax'])])] + [MetricsWithTimeout(outdir=OUTDIR / "Vision600_thorax", **DATA_SLICES['Siemens_Vision600_thorax'])]), + (SRCDIR / "Siemens_mMR_ACR", OUTDIR / "mMR_ACR", + [MetricsWithTimeout(outdir=OUTDIR / "mMR_ACR", **DATA_SLICES['Siemens_mMR_ACR'])]), + (SRCDIR / "Mediso_NEMA_IQ", OUTDIR / "Mediso_NEMA", + [MetricsWithTimeout(outdir=OUTDIR / "Mediso_NEMA", **DATA_SLICES['Mediso_NEMA_IQ'])]), + (SRCDIR / "GE_DMI3_Torso", OUTDIR / "DMI3_Torso", + [MetricsWithTimeout(outdir=OUTDIR / "DMI3_Torso", **DATA_SLICES['GE_DMI3_Torso'])])] else: log.warning("Source directory does not exist: %s", SRCDIR) data_dirs_metrics = [(None, None, [])] # type: ignore @@ -294,7 +306,7 @@ def get_image(fname): metrics_with_timeout.reset() # timeout from now algo = Submission(data) try: - algo.run(np.inf, callbacks=metrics + submission_callbacks) + algo.run(np.inf, callbacks=metrics + submission_callbacks, update_objective_interval=np.inf) except Exception: print_exc(limit=2) finally: