Merge pull request #61 from SyneRBI/data_QC_VOI

data_QC VOI improvements

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,28 +92,60 @@ 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
 
 
 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()