add QC for reference image and VOIs

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()