diff --git a/SIRF_data_preparation/data_QC.py b/SIRF_data_preparation/data_QC.py index 83338f1..e184341 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,24 @@ def plot_image(image, save_name=None, transverse_slice=-1, coronal_slice=-1, sag vmax = image.max() arr = image.as_array() - plt.figure() - 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[:, :, sagittal_slice], vmin=vmin, vmax=vmax) + 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] + + ax = plt.subplot(131) + plt.imshow(arr[transverse_slice, :, :], vmin=vmin, vmax=vmax, alpha=alpha_trans, **kwargs) + ax.set_title(f"T={transverse_slice}") + ax = plt.subplot(132) + plt.imshow(arr[:, coronal_slice, :], vmin=vmin, vmax=vmax, alpha=alpha_cor, **kwargs) + 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') @@ -79,9 +92,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 +104,48 @@ def VOI_mean(image, VOI): return float((image * VOI).sum() / VOI.sum()) -def VOI_checks(allVOInames, OSEM_image, reference_image, srcdir='.'): +from scipy import ndimage + +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: + 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=prefix, 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 + 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, alpha=allVOIs, save_name="OSEM_image_and_VOIs", **kwargs) + # unformatted print of VOI values for now print(allVOInames) print(OSEM_VOI_values) @@ -111,27 +156,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()