diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 9879945..98c256d 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,3 +1,6 @@ +ci: + skip: [todo] + default_language_version: python: python3 repos: diff --git a/SIRF_data_preparation/create_Hoffman_VOIs.py b/SIRF_data_preparation/create_Hoffman_VOIs.py index 5bf8398..19c297d 100644 --- a/SIRF_data_preparation/create_Hoffman_VOIs.py +++ b/SIRF_data_preparation/create_Hoffman_VOIs.py @@ -8,6 +8,7 @@ Options: -h, --help --dataset= dataset name (required) + --srcdir= pathname. Will default to current directory unless dataset is set -s, --skip_write_PETRIC_VOIs do not write in data//PETRIC """ @@ -30,9 +31,8 @@ import sirf.Reg as Reg import sirf.STIR as STIR -from petric import OUTDIR, SRCDIR from SIRF_data_preparation.data_QC import plot_image -from SIRF_data_preparation.data_utilities import the_orgdata_path +from SIRF_data_preparation.data_utilities import the_data_path, the_orgdata_path # %% __version__ = "0.1.0" @@ -47,24 +47,29 @@ if scanID is None: print("Need to set the --dataset argument") exit(1) - if args["--skip_write_PETRIC_VOIs"] is not None: + if args["--skip_write_PETRIC_VOIs"]: write_PETRIC_VOIs = False else: # set it by hand, e.g. scanID = "NeuroLF_Hoffman_Dataset" write_PETRIC_VOIs = False -# %% standard PETRIC directories -if not all((SRCDIR.is_dir(), OUTDIR.is_dir())): - PETRICDIR = Path("~/devel/PETRIC").expanduser() - SRCDIR = PETRICDIR / "data" - OUTDIR = PETRICDIR / "output" +srcdir = args['--srcdir'] +# %% standard PETRIC directories +if srcdir is None: + srcdir = the_data_path(scanID) +srcdir = Path(srcdir) downloaddir = Path(the_orgdata_path("downloads")) intermediate_data_path = the_orgdata_path(scanID, "processing") os.makedirs(downloaddir, exist_ok=True) os.makedirs(intermediate_data_path, exist_ok=True) +print("srcdir:", srcdir) +print("downloaddir:", downloaddir) +print("processingdir:", intermediate_data_path) +print("write_VOIs:", write_PETRIC_VOIs) + # %% Function to find the n-th connected component (counting by size) def connected_component(arr: npt.NDArray[bool], order=0) -> npt.NDArray[bool]: @@ -230,7 +235,6 @@ def create_Hoffman_VOIs(Hoffman: STIR.ImageData,) -> typing.Tuple[typing.List[ST # %%%%%%%%%%%%% Register VOIs to OSEM_image # %% Now register this to the reconstructed image -srcdir = SRCDIR / scanID OSEM_image = STIR.ImageData(str(srcdir / "OSEM_image.hv")) OSEM_image_nii = STIR_to_nii(OSEM_image, os.path.join(intermediate_data_path, "OSEM_image.nii")) @@ -263,6 +267,6 @@ def create_Hoffman_VOIs(Hoffman: STIR.ImageData,) -> typing.Tuple[typing.List[ST print(f"Writing registered VOIs to {datadir}") os.makedirs(datadir, exist_ok=True) for VOI, n in zip(regVOIs, VOInames): - VOI.write(str(datadir / n)) + VOI.write(str(datadir / ("VOI_"+n))) # %% plt.show() diff --git a/SIRF_data_preparation/data_QC.py b/SIRF_data_preparation/data_QC.py index 52d801c..ebee999 100644 --- a/SIRF_data_preparation/data_QC.py +++ b/SIRF_data_preparation/data_QC.py @@ -26,6 +26,7 @@ import matplotlib.pyplot as plt import numpy as np +import numpy.typing as npt from docopt import docopt from scipy import ndimage @@ -36,6 +37,15 @@ STIR.AcquisitionData.set_storage_scheme('memory') +def check_values_non_negative(arr: npt.NDArray[np.float32], desc: str): + min = np.min(arr) + max = np.max(arr) + if np.isnan(min) or min < 0: + raise ValueError(f"{desc}: minimum should be non-negative but is {min} (max={max})") + if not np.isfinite(max): + raise ValueError(f"{desc}: maximum should be finite") + + def plot_sinogram_profile(prompts, background, sumaxis=(0, 1), select=0, srcdir='./'): """ Plot a profile through sirf.STIR.AcquisitionData @@ -102,6 +112,13 @@ def plot_image_if_exists(prefix, **kwargs): return None +def check_and_plot_image_if_exists(prefix, **kwargs): + im = plot_image_if_exists(prefix, **kwargs) + if im is not None: + check_values_non_negative(im.as_array(), prefix) + return im + + def VOI_mean(image, VOI): return float((image * VOI).sum() / VOI.sum()) @@ -123,6 +140,7 @@ def VOI_checks(allVOInames, OSEM_image=None, reference_image=None, srcdir='.', * continue VOI = STIR.ImageData(filename) VOI_arr = VOI.as_array() + check_values_non_negative(VOI_arr, VOIname) COM = np.rint(ndimage.center_of_mass(VOI_arr)) num_voxels = VOI_arr.sum() print(f"VOI: {VOIname}: COM (in indices): {COM} voxels {num_voxels} = {num_voxels * np.prod(VOI.spacing)} mm^3") @@ -180,10 +198,14 @@ def main(argv=None): mult_factors = STIR.AcquisitionData(os.path.join(srcdir, 'mult_factors.hs')) background = additive_term * mult_factors plot_sinogram_profile(acquired_data, background, srcdir=srcdir) - - OSEM_image = plot_image_if_exists(os.path.join(srcdir, 'OSEM_image'), **slices) - plot_image_if_exists(os.path.join(srcdir, 'kappa'), **slices) - reference_image = plot_image_if_exists(os.path.join(srcdir, 'PETRIC/reference_image'), **slices) + check_values_non_negative(acquired_data.as_array(), "prompts") + check_values_non_negative(additive_term.as_array(), "additive_term") + check_values_non_negative(mult_factors.as_array(), "mult_factors") + check_values_non_negative(background.as_array(), "background") + + OSEM_image = check_and_plot_image_if_exists(os.path.join(srcdir, 'OSEM_image'), **slices) + check_and_plot_image_if_exists(os.path.join(srcdir, 'kappa'), **slices) + reference_image = check_and_plot_image_if_exists(os.path.join(srcdir, 'PETRIC/reference_image'), **slices) VOIdir = os.path.join(srcdir, 'PETRIC') allVOInames = [os.path.basename(str(voi)[:-3]) for voi in Path(VOIdir).glob("VOI_*.hv")] diff --git a/SIRF_data_preparation/dataset_settings.py b/SIRF_data_preparation/dataset_settings.py index 1f8f827..74485cf 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, 'GE_DMI3_Torso': 8} + 'Mediso_NEMA_IQ': 12, 'Siemens_Vision600_thorax': 5, 'GE_DMI3_Torso': 8, 'Siemens_Vision600_Hoffman': 5} @dataclass diff --git a/SIRF_data_preparation/get_penalisation_factor.py b/SIRF_data_preparation/get_penalisation_factor.py index a6ce158..030d4ed 100755 --- a/SIRF_data_preparation/get_penalisation_factor.py +++ b/SIRF_data_preparation/get_penalisation_factor.py @@ -39,8 +39,8 @@ if dataset is None or ref_dataset is None: print("Need to set the --dataset arguments") exit(1) - if args["--write_penalisation_factor"] is not None: - write_penalisation_factor = False + write_penalisation_factor = args["--write_penalisation_factor"] + else: # set it by hand, e.g. ref_dataset = "NeuroLF_Hoffman_Dataset" dataset = "Siemens_mMR_NEMA_IQ" diff --git a/SIRF_data_preparation/run_BSREM.py b/SIRF_data_preparation/run_BSREM.py index 9181681..d6a00ec 100644 --- a/SIRF_data_preparation/run_BSREM.py +++ b/SIRF_data_preparation/run_BSREM.py @@ -4,7 +4,10 @@ run_BSREM.py [--help | options] Arguments: - path to data files as well as prefix to use (e.g. Siemens_mMR_NEMA_EQ) + path to data files as well as prefix to use (e.g. Siemens_mMR_NEMA_EQ) + +Options: + --updates= number of updates to run [default: 15000] """ # Copyright 2024 Rutherford Appleton Laboratory STFC @@ -28,6 +31,7 @@ # logging.basicConfig(level=logging.INFO) scanID = args[''] +num_updates = int(args['--updates']) if not all((SRCDIR.is_dir(), OUTDIR.is_dir())): PETRICDIR = Path('~/devel/PETRIC').expanduser() @@ -41,6 +45,9 @@ settings = get_settings(scanID) data = get_data(srcdir=srcdir, outdir=outdir) +print("Penalisation factor:", data.prior.get_penalisation_factor()) +print("num_subsets:", settings.num_subsets) +print("num_updates:", num_updates) data_sub, acq_models, obj_funs = partitioner.data_partition(data.acquired_data, data.additive_term, data.mult_factors, settings.num_subsets, mode="staggered", initial_image=data.OSEM_image) @@ -53,7 +60,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(num_updates, callbacks=[MetricsWithTimeout(**settings.slices, interval=80, outdir=outdir, seconds=3600 * 100)]) # %% fig = plt.figure() data_QC.plot_image(algo.get_output(), **settings.slices) diff --git a/SIRF_data_preparation/run_OSEM.py b/SIRF_data_preparation/run_OSEM.py index 00885ee..d628dad 100644 --- a/SIRF_data_preparation/run_OSEM.py +++ b/SIRF_data_preparation/run_OSEM.py @@ -6,12 +6,13 @@ Arguments: path to data files as well as prefix to use +Options: + --updates= number of updates to run [default: 400] """ # 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 @@ -27,6 +28,7 @@ # logging.basicConfig(level=logging.INFO) scanID = args[''] +num_updates = int(args['--updates']) if not all((SRCDIR.is_dir(), OUTDIR.is_dir())): PETRICDIR = Path('~/devel/PETRIC').expanduser() @@ -40,9 +42,11 @@ settings = get_settings(scanID) data = get_data(srcdir=srcdir, outdir=outdir) +print("num_subsets:", settings.num_subsets) +print("num_updates:", num_updates) # %% algo = main_OSEM.Submission(data, settings.num_subsets, update_objective_interval=20) -algo.run(400, callbacks=[MetricsWithTimeout(**settings.slices, seconds=5000, outdir=outdir)]) +algo.run(num_updates, callbacks=[MetricsWithTimeout(**settings.slices, seconds=5000, interval=20, outdir=outdir)]) # %% fig = plt.figure() data_QC.plot_image(algo.get_output(), **settings.slices) diff --git a/petric.py b/petric.py index e0c5ff5..d86ce77 100755 --- a/petric.py +++ b/petric.py @@ -178,7 +178,7 @@ def pass_index(metrics: np.ndarray, thresh: Iterable, window: int = 10) -> int: return np.where(res)[0][0] -class MetricsWithTimeout(cil_callbacks.Callback): +class MetricsWithTimeout(Callback): """Stops the algorithm after `seconds`""" def __init__(self, seconds=3600, outdir=OUTDIR, transverse_slice=None, coronal_slice=None, sagittal_slice=None, **kwargs): @@ -186,9 +186,9 @@ def __init__(self, seconds=3600, outdir=OUTDIR, transverse_slice=None, coronal_s self._seconds = seconds self.callbacks = [ cil_callbacks.ProgressCallback(), - SaveIters(outdir=outdir), + SaveIters(outdir=outdir, **kwargs), (tb_cbk := StatsLog(logdir=outdir, transverse_slice=transverse_slice, coronal_slice=coronal_slice, - sagittal_slice=sagittal_slice))] + sagittal_slice=sagittal_slice, **kwargs))] self.tb = tb_cbk.tb # convenient access to the underlying SummaryWriter self.reset() @@ -292,8 +292,8 @@ 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': {}} + 'Mediso_NEMA_IQ': {'transverse_slice': 22, 'coronal_slice': 89, 'sagittal_slice': 66}, + 'Siemens_Vision600_thorax': {}, 'GE_DMI3_Torso': {}, 'Siemens_Vision600_Hoffman': {}} if SRCDIR.is_dir(): # create list of existing data