From e21ff315f9ac34547653c816a797efede36ce781 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Sat, 13 Jul 2024 10:48:16 +0100 Subject: [PATCH 1/3] data_QC: plot VOIs with transparancy as well also remove need for trialing slah in srcdir --- SIRF_data_preparation/data_QC.py | 81 ++++++++++++++++++++++---------- 1 file changed, 56 insertions(+), 25 deletions(-) diff --git a/SIRF_data_preparation/data_QC.py b/SIRF_data_preparation/data_QC.py index 83338f1..0349a0f 100644 --- a/SIRF_data_preparation/data_QC.py +++ b/SIRF_data_preparation/data_QC.py @@ -46,12 +46,14 @@ def plot_sinogram_profile(prompts, background, sumaxis=(0, 1), select=0, srcdir= plt.plot(np.sum(prompts.as_array(), axis=sumaxis)[select, :], label='prompts') plt.plot(np.sum(background.as_array(), axis=sumaxis)[select, :], label='background') ax.legend() - plt.savefig(srcdir + 'prompts_background_profiles.png') + plt.savefig(os.path.join(srcdir, 'prompts_background_profiles.png')) -def plot_image(image, save_name=None, transverse_slice=-1, coronal_slice=-1, sagittal_slice=-1, vmin=0, vmax=None): +def plot_image(image, save_name=None, + transverse_slice=-1, coronal_slice=-1, sagittal_slice=-1, + vmin=0, vmax=None, alpha=None, **kwargs): """ - Plot a profile through sirf.STIR.ImageData + Plot transverse/coronal/sagital slices through sirf.STIR.ImageData """ if transverse_slice < 0: transverse_slice = image.dimensions()[0] // 2 @@ -63,13 +65,21 @@ def plot_image(image, save_name=None, transverse_slice=-1, coronal_slice=-1, sag vmax = image.max() arr = image.as_array() - plt.figure() + alpha_trans = None + alpha_cor = None + alpha_sag = None + if alpha is not None: + alpha_arr = alpha.as_array() + alpha_trans = alpha_arr[transverse_slice, :, :] + alpha_cor = alpha_arr[:, coronal_slice, :] + alpha_sag = alpha_arr[:, :, sagittal_slice] + plt.subplot(131) - plt.imshow(arr[transverse_slice, :, :], vmin=vmin, vmax=vmax) + plt.imshow(arr[transverse_slice, :, :], vmin=vmin, vmax=vmax, alpha=alpha_trans, **kwargs) plt.subplot(132) - plt.imshow(arr[:, coronal_slice, :], vmin=vmin, vmax=vmax) + plt.imshow(arr[:, coronal_slice, :], vmin=vmin, vmax=vmax, alpha=alpha_cor, **kwargs) plt.subplot(133) - plt.imshow(arr[:, :, sagittal_slice], vmin=vmin, vmax=vmax) + plt.imshow(arr[:, :, sagittal_slice], vmin=vmin, vmax=vmax, alpha=alpha_sag, **kwargs) plt.colorbar(shrink=.6) if save_name is not None: plt.savefig(save_name + '_slices.png') @@ -79,9 +89,11 @@ def plot_image(image, save_name=None, transverse_slice=-1, coronal_slice=-1, sag def plot_image_if_exists(prefix, **kwargs): if os.path.isfile(prefix + '.hv'): im = STIR.ImageData(prefix + '.hv') + plt.figure() plot_image(im, prefix, **kwargs) return im else: + print(f"Image {prefix}.hv does not exist") return None @@ -89,18 +101,39 @@ def VOI_mean(image, VOI): return float((image * VOI).sum() / VOI.sum()) -def VOI_checks(allVOInames, OSEM_image, reference_image, srcdir='.'): +def VOI_checks(allVOInames, OSEM_image=None, reference_image=None, srcdir='.', **kwargs): if len(allVOInames) == 0: return OSEM_VOI_values = [] ref_VOI_values = [] + allVOIs = None + VOIkwargs = kwargs.copy() + VOIkwargs['vmax'] = 1 + VOIkwargs['vmin'] = 0 for VOIname in allVOInames: - VOI = plot_image_if_exists(os.path.join(srcdir, VOIname), transverse_slice=transverse_slice, - coronal_slice=coronal_slice, sagittal_slice=sagittal_slice) - if OSEM_image: + VOI = plot_image_if_exists(os.path.join(srcdir, VOIname), **VOIkwargs) + if VOI is None: + print(f"VOI {VOIname} does not exist") + continue + # construct transparency image + if VOIname == 'VOI_whole_object': + VOI /= 2 + if allVOIs is None: + allVOIs = VOI.clone() + else: + allVOIs += VOI + if OSEM_image is not None: OSEM_VOI_values.append(VOI_mean(OSEM_image, VOI)) if reference_image: ref_VOI_values.append(VOI_mean(reference_image, VOI)) + allVOIs /= allVOIs.max() + + if OSEM_image is not None: + plt.figure() + plot_image(OSEM_image, **kwargs) + plt.figure() + plot_image(OSEM_image, alpha=allVOIs, **kwargs) + # unformatted print of VOI values for now print(allVOInames) print(OSEM_VOI_values) @@ -111,27 +144,25 @@ def main(argv=None): args = docopt(__doc__, argv=argv, version=__version__) srcdir = args['--srcdir'] skip_sino_profiles = args['--skip_sino_profiles'] - transverse_slice = literal_eval(args['--transverse_slice']) - coronal_slice = literal_eval(args['--coronal_slice']) - sagittal_slice = literal_eval(args['--sagittal_slice']) + slices = {} + slices["transverse_slice"] = literal_eval(args['--transverse_slice']) + slices["coronal_slice"] = literal_eval(args['--coronal_slice']) + slices["sagittal_slice"] = literal_eval(args['--sagittal_slice']) if not skip_sino_profiles: - acquired_data = STIR.AcquisitionData(srcdir + 'prompts.hs') - additive_term = STIR.AcquisitionData(srcdir + 'additive_term.hs') - mult_factors = STIR.AcquisitionData(srcdir + 'mult_factors.hs') + acquired_data = STIR.AcquisitionData(os.path.join(srcdir, 'prompts.hs')) + additive_term = STIR.AcquisitionData(os.path.join(srcdir, 'additive_term.hs')) + 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('OSEM_image', transverse_slice=transverse_slice, coronal_slice=coronal_slice, - sagittal_slice=sagittal_slice) - plot_image_if_exists('kappa', transverse_slice=transverse_slice, coronal_slice=coronal_slice, - sagittal_slice=sagittal_slice) - reference_image = plot_image_if_exists(srcdir + 'PETRIC/reference_image', transverse_slice=transverse_slice, - coronal_slice=coronal_slice, sagittal_slice=sagittal_slice) + 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) - VOIdir = os.path.join(srcdir, srcdir + "PETRIC") + 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) + VOI_checks(allVOInames, OSEM_image, reference_image, srcdir=VOIdir, **slices) plt.show() From 2a741e9d19149695545bbfe64f1d421253b74dfe Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Sat, 13 Jul 2024 11:04:09 +0100 Subject: [PATCH 2/3] data_QC: centre VOI images and add COM to plots --- SIRF_data_preparation/data_QC.py | 27 +++++++++++++++++++-------- 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/SIRF_data_preparation/data_QC.py b/SIRF_data_preparation/data_QC.py index 0349a0f..61626f0 100644 --- a/SIRF_data_preparation/data_QC.py +++ b/SIRF_data_preparation/data_QC.py @@ -74,12 +74,15 @@ def plot_image(image, save_name=None, alpha_cor = alpha_arr[:, coronal_slice, :] alpha_sag = alpha_arr[:, :, sagittal_slice] - plt.subplot(131) + ax = plt.subplot(131) plt.imshow(arr[transverse_slice, :, :], vmin=vmin, vmax=vmax, alpha=alpha_trans, **kwargs) - plt.subplot(132) + ax.set_title(f"T={transverse_slice}") + ax = plt.subplot(132) plt.imshow(arr[:, coronal_slice, :], vmin=vmin, vmax=vmax, alpha=alpha_cor, **kwargs) - plt.subplot(133) + ax.set_title(f"C={coronal_slice}") + ax = plt.subplot(133) plt.imshow(arr[:, :, sagittal_slice], vmin=vmin, vmax=vmax, alpha=alpha_sag, **kwargs) + ax.set_title(f"S={sagittal_slice}") plt.colorbar(shrink=.6) if save_name is not None: plt.savefig(save_name + '_slices.png') @@ -101,6 +104,8 @@ def VOI_mean(image, VOI): return float((image * VOI).sum() / VOI.sum()) +from scipy import ndimage + def VOI_checks(allVOInames, OSEM_image=None, reference_image=None, srcdir='.', **kwargs): if len(allVOInames) == 0: return @@ -111,10 +116,18 @@ def VOI_checks(allVOInames, OSEM_image=None, reference_image=None, srcdir='.', * VOIkwargs['vmax'] = 1 VOIkwargs['vmin'] = 0 for VOIname in allVOInames: - VOI = plot_image_if_exists(os.path.join(srcdir, VOIname), **VOIkwargs) - if VOI is None: + filename = os.path.join(srcdir, VOIname + '.hv') + if not os.path.isfile(filename): print(f"VOI {VOIname} does not exist") continue + VOI = STIR.ImageData(filename) + COM = np.rint(ndimage.center_of_mass(VOI.as_array())) + plt.figure() + plot_image(VOI, save_name=VOIname, vmin=0, vmax=1, + transverse_slice = int(COM[0]), + coronal_slice = int(COM[1]), + sagittal_slice = int(COM[2])) + # construct transparency image if VOIname == 'VOI_whole_object': VOI /= 2 @@ -130,9 +143,7 @@ def VOI_checks(allVOInames, OSEM_image=None, reference_image=None, srcdir='.', * if OSEM_image is not None: plt.figure() - plot_image(OSEM_image, **kwargs) - plt.figure() - plot_image(OSEM_image, alpha=allVOIs, **kwargs) + plot_image(OSEM_image, alpha=allVOIs, save_name="OSEM_image and VOIs", **kwargs) # unformatted print of VOI values for now print(allVOInames) From c5d554cfaabc22ac13f5df9c5d80a2b0ac5eed48 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Sat, 13 Jul 2024 11:17:12 +0100 Subject: [PATCH 3/3] data_QC: fix location of VOI*.png --- SIRF_data_preparation/data_QC.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/SIRF_data_preparation/data_QC.py b/SIRF_data_preparation/data_QC.py index 61626f0..e184341 100644 --- a/SIRF_data_preparation/data_QC.py +++ b/SIRF_data_preparation/data_QC.py @@ -116,14 +116,15 @@ def VOI_checks(allVOInames, OSEM_image=None, reference_image=None, srcdir='.', * VOIkwargs['vmax'] = 1 VOIkwargs['vmin'] = 0 for VOIname in allVOInames: - filename = os.path.join(srcdir, VOIname + '.hv') + prefix = os.path.join(srcdir, VOIname) + filename = prefix + '.hv' if not os.path.isfile(filename): print(f"VOI {VOIname} does not exist") continue VOI = STIR.ImageData(filename) COM = np.rint(ndimage.center_of_mass(VOI.as_array())) plt.figure() - plot_image(VOI, save_name=VOIname, vmin=0, vmax=1, + 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])) @@ -143,7 +144,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="OSEM_image_and_VOIs", **kwargs) # unformatted print of VOI values for now print(allVOInames)