diff --git a/SIRF_data_preparation/Mediso_NEMA_IQ/README.md b/SIRF_data_preparation/Mediso_NEMA_IQ/README.md new file mode 100644 index 0000000..670bce5 --- /dev/null +++ b/SIRF_data_preparation/Mediso_NEMA_IQ/README.md @@ -0,0 +1,36 @@ +# Data from the Mediso AnyScan at NPL +Initial processing steps +```sh +orgpath=/mnt/share/petric-wip/NEMA_Challenge +#orgpath=~/devel/PETRIC/orgdata/ +# cd $orgpath +#wget https://nplanywhere.npl.co.uk/userportal/?v=4.6.1#/shared/public/nrKe5vWynvyMsp7m/e9e3d4b8-8f83-47be-bf2d-17a96a7a8aac +#unzip NEMA_Challenge.zip +cd ~/devel/PETRIC/data +mkdir -p Mediso_NEMA_IQ +cd Mediso_NEMA_IQ/ +# trim sinograms to avoid "corner" problems in mult_factors +# TODO for next data: use 30 (as some problems remain) +for f in additive_term.hs mult_factors.hs prompts.hs; do + SSRB -t 20 $f $orgpath/sinograms/$f; +done +# alternative if we don't need to trim +#cp -rp $orgpath/sinograms/* . +# get rid of NaNs +python prepare.py +# now handle VOIs +mkdir PETRIC +cp -rp $orgpath/VOIs/* PETRIC +# cp -rp $orgpath/README.md . +cd PETRIC +stir_math VOI_background.hv VOI_backgroung.hv +rm VOI_backgroung.* +rm *ahv +cd .. +python ../../SIRF_data_preparation/create_initial_images.py --template_image=PETRIC/VOI_whole_object.hv . + ``` + I needed to fix the OSEM header (manual): + ``` + !imaging modality := PT + ``` +python ../../SIRF_data_preparation/data_QC diff --git a/SIRF_data_preparation/Mediso_NEMA_IQ/prepare.py b/SIRF_data_preparation/Mediso_NEMA_IQ/prepare.py new file mode 100644 index 0000000..68c2f42 --- /dev/null +++ b/SIRF_data_preparation/Mediso_NEMA_IQ/prepare.py @@ -0,0 +1,10 @@ +# Set NaNs to zero in additive_term +import sirf.STIR +import numpy as np +additive = sirf.STIR.AcquisitionData('additive_term.hs') +add_arr = additive.as_array() +add_arr = np.nan_to_num(add_arr) +new_add = additive.clone() +new_add.fill(add_arr) +new_add.write('additive_term.hs') + diff --git a/SIRF_data_preparation/data_QC.py b/SIRF_data_preparation/data_QC.py index 77401e1..cdffe03 100644 --- a/SIRF_data_preparation/data_QC.py +++ b/SIRF_data_preparation/data_QC.py @@ -11,12 +11,13 @@ --transverse_slice= idx [default: -1] --coronal_slice= idx [default: -1] --sagittal_slice= idx [default: -1] + --dataset= dataset name. if set, it is used to override default slices Note that -1 one means to use middle of image """ # Copyright 2024 University College London # Licence: Apache-2.0 -__version__ = '0.2.0' +__version__ = '0.3.0' import os import os.path @@ -29,6 +30,7 @@ from scipy import ndimage import sirf.STIR as STIR +from SIRF_data_preparation.dataset_settings import get_settings STIR.AcquisitionData.set_storage_scheme('memory') @@ -119,7 +121,10 @@ def VOI_checks(allVOInames, OSEM_image=None, reference_image=None, srcdir='.', * print(f"VOI {VOIname} does not exist") continue VOI = STIR.ImageData(filename) - COM = np.rint(ndimage.center_of_mass(VOI.as_array())) + VOI_arr = VOI.as_array() + 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") plt.figure() plot_image(VOI, save_name=prefix, vmin=0, vmax=1, transverse_slice=int(COM[0]), coronal_slice=int(COM[1]), sagittal_slice=int(COM[2])) @@ -149,6 +154,7 @@ def VOI_checks(allVOInames, OSEM_image=None, reference_image=None, srcdir='.', * def main(argv=None): args = docopt(__doc__, argv=argv, version=__version__) + dataset = args['--dataset'] srcdir = args['--srcdir'] skip_sino_profiles = args['--skip_sino_profiles'] slices = {} @@ -156,6 +162,13 @@ def main(argv=None): slices["coronal_slice"] = literal_eval(args['--coronal_slice']) slices["sagittal_slice"] = literal_eval(args['--sagittal_slice']) + if (dataset): + settings = get_settings(dataset) + for key in slices.keys(): + if slices[key] == -1 and key in settings.slices: + slices[key] = settings.slices[key] + print(slices) + if not skip_sino_profiles: acquired_data = STIR.AcquisitionData(os.path.join(srcdir, 'prompts.hs')) additive_term = STIR.AcquisitionData(os.path.join(srcdir, 'additive_term.hs')) diff --git a/SIRF_data_preparation/dataset_settings.py b/SIRF_data_preparation/dataset_settings.py index 0891e24..25b6ab8 100644 --- a/SIRF_data_preparation/dataset_settings.py +++ b/SIRF_data_preparation/dataset_settings.py @@ -19,6 +19,9 @@ def get_settings(scanID: str): elif scanID == 'NeuroLF_Hoffman_Dataset': slices = {'transverse_slice': 72} num_subsets = 16 + elif scanID == 'Mediso_NEMA_IQ': + slices = {'transverse_slice': 22, 'coronal_slice': 89, 'sagittal_slice': 66} + num_subsets = 12 else: # Vision slices = {} num_subsets = 5 diff --git a/SIRF_data_preparation/plot_BSREM_metrics.py b/SIRF_data_preparation/plot_BSREM_metrics.py index 6ff643f..a5d0606 100644 --- a/SIRF_data_preparation/plot_BSREM_metrics.py +++ b/SIRF_data_preparation/plot_BSREM_metrics.py @@ -1,4 +1,5 @@ -"""Preliminary file to check evolution of metrics as well as pass_index""" +#%% +# """Preliminary file to check evolution of metrics as well as pass_index""" # %load_ext autoreload # %autoreload 2 from pathlib import Path @@ -21,7 +22,8 @@ # scanID = 'Siemens_Vision600_thorax' # scanID = 'NeuroLF_Hoffman_Dataset' -scanID = 'Siemens_mMR_NEMA_IQ' +# scanID = 'Siemens_mMR_NEMA_IQ' +scanID = 'Mediso_NEMA_IQ' srcdir = SRCDIR / scanID outdir = OUTDIR / scanID @@ -69,8 +71,10 @@ 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 last_iteration = int(objs[-1, 0] + .5) -iteration_interval = int(objs[-1, 0] - objs[-2, 0] + .5) +# find interval(don't use last value, as that interval can be smaller) +iteration_interval = int(objs[-2, 0] - objs[-3, 0] + .5) if datadir1.is_dir(): last_iteration = int(objs0[-1, 0] + .5) iteration_interval = int(objs0[-1, 0] - objs0[-2, 0] + .5) * 2 @@ -78,7 +82,8 @@ iters = range(0, last_iteration, iteration_interval) m = get_metrics(qm, iters, srcdir=datadir) # %% -OSEMiters = range(0, 400, 20) +OSEMobjs = read_objectives(OSEMdir) +OSEMiters = OSEMobjs[:, 0].astype(int) OSEMm = get_metrics(qm, OSEMiters, srcdir=OSEMdir) # %% fig = plt.figure() @@ -107,7 +112,7 @@ fig.savefig(outdir / f'{scanID}_metrics_BSREM.png') # %% -idx = pass_index(m, numpy.array([.01, .01, .005, .005, .005]), 10) +idx = pass_index(m, numpy.array([.01, .01] + [.005 for i in range(len(data.voi_masks))]), 10) iter = iters[idx] print(iter) image = STIR.ImageData(str(datadir / f"iter_{iter:04d}.hv")) diff --git a/SIRF_data_preparation/run_OSEM.py b/SIRF_data_preparation/run_OSEM.py index 500e63e..00885ee 100644 --- a/SIRF_data_preparation/run_OSEM.py +++ b/SIRF_data_preparation/run_OSEM.py @@ -42,7 +42,7 @@ data = get_data(srcdir=srcdir, outdir=outdir) # %% algo = main_OSEM.Submission(data, settings.num_subsets, update_objective_interval=20) -algo.run(20, callbacks=[MetricsWithTimeout(**settings.slices, seconds=5000, outdir=outdir)]) +algo.run(400, callbacks=[MetricsWithTimeout(**settings.slices, seconds=5000, outdir=outdir)]) # %% fig = plt.figure() data_QC.plot_image(algo.get_output(), **settings.slices) diff --git a/petric.py b/petric.py index d23a230..18fb47a 100755 --- a/petric.py +++ b/petric.py @@ -141,7 +141,7 @@ def keys(self): class MetricsWithTimeout(cil_callbacks.Callback): """Stops the algorithm after `seconds`""" - def __init__(self, seconds=600, outdir=OUTDIR, transverse_slice=None, coronal_slice=None, **kwargs): + def __init__(self, seconds=600, outdir=OUTDIR, transverse_slice=None, coronal_slice=None, sagittal_slice=None, **kwargs): super().__init__(**kwargs) self._seconds = seconds self.callbacks = [