From 479a16b51fe87296277f6d87eac94223a34ee8c0 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Wed, 10 Jul 2024 16:11:15 +0100 Subject: [PATCH] add QC for reference image and VOIs --- SIRF_data_preparation/data_QC.py | 68 ++++++++++++++++++++++++-------- 1 file changed, 52 insertions(+), 16 deletions(-) diff --git a/SIRF_data_preparation/data_QC.py b/SIRF_data_preparation/data_QC.py index 7e7eeee..adade4b 100644 --- a/SIRF_data_preparation/data_QC.py +++ b/SIRF_data_preparation/data_QC.py @@ -9,16 +9,18 @@ """ # Copyright 2024 University College London # Licence: Apache-2.0 -__version__ = '0.1.0' +__version__ = '0.2.0' import os import matplotlib.pyplot as plt import numpy as np from docopt import docopt - +from pathlib import Path +import os.path import sirf.STIR as STIR +STIR.AcquisitionData.set_storage_scheme('memory') def plot_sinogram_profile(prompts, background, sumaxis=(0, 1), select=0): """ @@ -36,48 +38,82 @@ def plot_sinogram_profile(prompts, background, sumaxis=(0, 1), select=0): plt.savefig('prompts_background_profiles.png') -def plot_image(image, save_name, transverse_slice=-1, coronal_slice=-1): +def plot_image(image, save_name, + transverse_slice=-1, + coronal_slice=-1, + sagital_slice=-1, + vmin=0, vmax=None): """ Plot a profile through sirf.STIR.ImageData """ - plt.figure() if transverse_slice < 0: transverse_slice = image.dimensions()[0] // 2 if coronal_slice < 0: coronal_slice = image.dimensions()[1] // 2 - cmax = image.max() + if sagital_slice < 0: + sagital_slice = image.dimensions()[2] // 2 + if vmax is None: + vmax = image.max() + arr = image.as_array() plt.figure() - plt.subplot(121) - plt.imshow(arr[transverse_slice, :, :], vmax=cmax) - plt.subplot(122) - plt.imshow(arr[:, coronal_slice, :], vmax=cmax) - plt.colorbar() + plt.subplot(131) + plt.imshow(arr[transverse_slice, :, :], vmin=vmin, vmax=vmax) + plt.subplot(132) + plt.imshow(arr[:, coronal_slice, :], vmin=vmin, vmax=vmax) + plt.subplot(133) + plt.imshow(arr[:,:, sagital_slice], vmin=vmin, vmax=vmax) + plt.colorbar(shrink=.6) plt.savefig(save_name + '_slices.png') + plt.suptitle(os.path.basename(save_name)) def plot_image_if_exists(prefix): if os.path.isfile(prefix + '.hv'): im = STIR.ImageData(prefix + '.hv') plot_image(im, prefix) - - -STIR.AcquisitionData.set_storage_scheme('memory') + return im + else: + return None + +def VOI_mean(image, VOI): + return float((image*VOI).sum() / VOI.sum()) + +def VOI_checks(allVOInames, OSEM_image, reference_image, srcdir='.'): + if len(allVOInames) == 0: + return + OSEM_VOI_values = [] + ref_VOI_values = [] + for VOIname in allVOInames: + VOI = plot_image_if_exists(os.path.join(srcdir, VOIname)) + if OSEM_image: + OSEM_VOI_values.append(VOI_mean(OSEM_image, VOI)) + if reference_image: + ref_VOI_values.append(VOI_mean(reference_image, VOI)) + # unformatted print of VOI values for now + print(allVOInames) + print(OSEM_VOI_values) + print(ref_VOI_values) def main(argv=None): args = docopt(__doc__, argv=argv, version=__version__) - acquired_data = STIR.AcquisitionData('prompts.hs') + srcdir = './' + acquired_data = STIR.AcquisitionData(srcdir+'prompts.hs') additive_term = STIR.AcquisitionData('additive_term.hs') mult_factors = STIR.AcquisitionData('mult_factors.hs') background = additive_term * mult_factors plot_sinogram_profile(acquired_data, background) - plot_image_if_exists('OSEM_image') + OSEM_image = plot_image_if_exists('OSEM_image') plot_image_if_exists('kappa') - plt.show() + reference_image = plot_image_if_exists('reference_image') + VOIdir = os.path.join(srcdir, "PETRIC") + allVOInames = [ os.path.basename(str(voi)[:-3]) for voi in Path(VOIdir).glob("VOI_*.hv")] + VOI_checks(allVOInames, OSEM_image, reference_image, srcdir=VOIdir) + plt.show() if __name__ == '__main__': main()