SRT2D repository: SRT 2D PET and SPECT algorithms (#1420)

* Added SRT2DPETReconstruction and SRT2DSPECTReconstruction
* Added Wiener and Gamma filters
* Expanded recon_test_pack to run tests for SPECT data with SRT and others. Therefore,
renamed simulate_PET_data_for_tests.sh to simulate_data_for_tests.sh

---------

Co-authored-by: Kris Thielemans <[email protected]>

.github/workflows/build-test.yml

@@ -234,13 +234,12 @@ jobs:
               # brew install openblas
               # export OPENBLAS=$(brew --prefix openblas)
               #python -m pip install --no-cache-dir --no-binary  numpy numpy # avoid the cached .whl!
-              python -m pip install numpy
+              python -m pip install numpy pytest
               ;;
             (*)
-              python -m pip install numpy
+              python -m pip install numpy pytest
               ;;
           esac
-          python -m pip install pytest
 
           if test "${{matrix.parallelproj}}XX" == "ONXX"; then
              git clone --depth 1 --branch v1.7.3 https://github.com/gschramm/parallelproj

documentation/release_6.3.htm

@@ -19,6 +19,21 @@ <h2> Summary for end users (also to be read by developers)</h2>
 
 <h3>New functionality</h3>
 <ul>
+  <li>
+    The analytic Spline Reconstruction Technique (SRT) algorithm has been added in 2 different versions: for PET
+    (inverting the 2D Radon transform) and for SPECT (inverting the 2D attenuated Radon transform). The latter
+    allows quantitatively correct analytic reconstruction of SPECT data (after scatter correction).
+
+    The reference for the implemented algorithms is:<br>
+    Fokas, A. S., A. Iserles, and V. Marinakis.
+    <cite>Reconstruction algorithm for single photon emission computed tomography and its numerical implementation.</cite>
+    Journal of the Royal Society Interface* 3.6 (2006): 45-54.
+    <br>
+    The reference for the implementation is:<br>
+    Dimitra Kyriakopoulou, <cite>Analytical and Numerical Aspects of Tomography</cite>, UCL PhD thesis, (not yet publically accessible)
+    <br>
+    <a href=https://github.com/UCL/STIR/pull/1420>PR #1420</a>
+  </li>
   <li>
     <code>ScatterSimulation</code> can now downsample the scanner transaxially (crystals per ring) for <code>BlocksOnCylindrical</code>,
     scanners, which speeds up <code>ScatterEstimation</code> considerably. By default, downsampling the detectors per reading

examples/samples/SRT2D.par

@@ -0,0 +1,18 @@
+SRT2DParameters :=
+input file := sino_hof_no_noise.hs
+
+xy output image size (in pixels) := -1  
+z output image size (in pixels) := -1
+zoom := 1 
+
+;post-filter type := Wiener
+;Wiener Filter Parameters :=
+;End Wiener Filter Parameters :=
+
+;post-filter type := Gamma
+;Gamma Filter Parameters :=
+;End Gamma Filter Parameters :=
+
+output filename prefix := SRT2D
+
+END :=

examples/samples/SRT2DSPECT.par

@@ -0,0 +1,20 @@
+SRT2DSPECTParameters :=
+
+;SRT2DSPECT receives as input the forward projection of the attenuation image (attensino.hs) and the measured emission sinogram (sino.hs).
+input file := sino.hs
+attenuation projection filename := attensino.hs
+
+xy output image size (in pixels) := -1
+
+;post-filter type := Wiener
+;Wiener Filter Parameters :=
+;End Wiener Filter Parameters :=
+
+;post-filter type := Gamma
+;Gamma Filter Parameters :=
+;End Gamma Filter Parameters :=
+
+output filename prefix := srt_recon
+
+
+end:=

recon_test_pack/SPECT_test_Interfile_header.hs

@@ -0,0 +1,62 @@
+!INTERFILE  :=
+; This is a sample minimal header for SPECT tomographic data
+; The format is as per the 3.3 Interfile standard (aside from time frame info)
+
+!imaging modality := nucmed
+
+; name of file with binary data
+name of data file := SPECT_test_Interfile_header.s
+
+!version of keys := 3.3
+!GENERAL DATA :=
+!GENERAL IMAGE DATA :=
+!type of data := Tomographic
+
+; optional keywords specifying patient position
+; patient rotation := prone
+; patient orientation := feet_in
+
+imagedata byte order := LITTLEENDIAN
+
+number of energy windows:=1
+energy window lower level[1]:=120
+energy window upper level[1]:=160
+
+!SPECT STUDY (General) :=
+; specify how the data are stored on disk
+; here given as "single-precision float" (you could have "unsigned integer" data instead)
+!number format := float
+!number of bytes per pixel := 4
+!number of projections := 120
+; total rotation (or coverage) angle (in degrees)
+!extent of rotation := 360
+process status := acquired
+!SPECT STUDY (acquired data):=
+; rotation info (e.g. clock-wise or counter-clock wise)
+!direction of rotation := CW
+start angle := 180
+
+; Orbit definition
+orbit := Circular
+; radius in mm
+Radius := 166.5
+; or
+; orbit := Non-circular
+; give a list of "radii", one for every position
+; Radii := {150, 151, 153, ....}
+
+; pixel sizes in the acquired data, first in "transverse" direction, then in "axial" direction
+!matrix size [1] := 111
+!scaling factor (mm/pixel) [1] := 3
+!matrix size [2] := 47
+!scaling factor (mm/pixel) [2] := 3.27
+
+; optional keywords specifying frame duration etc
+; These are not according to the Interfile 3.3 specification
+; Currently only useful in STIR for dynamic applications
+; (but a "time frame" is considered to be all projections acquired at the same time)
+;number of time frames := 1
+;image duration (sec)[1] := 0
+;image relative start time (sec)[1] := 0
+
+!END OF INTERFILE :=

recon_test_pack/SPECT_test_Interfile_header.s

@@ -0,0 +1 @@
+

recon_test_pack/SRT2DSPECT_test_sim.par

@@ -0,0 +1,20 @@
+SRT2DSPECTParameters :=
+; test file for SRT2DSPECT
+input file := my_sino${suffix}.hs  
+attenuation projection filename := my_attenuation_sino${suffix}.hs
+
+xy output image size (in pixels) := 91
+zoom := 0.5
+
+;post-filter type := Wiener
+;Wiener Filter Parameters :=
+;End Wiener Filter Parameters :=
+
+;post-filter type := Gamma
+;Gamma Filter Parameters :=
+;End Gamma Filter Parameters :=
+
+output filename prefix := my_test_sim_image_SRT2DSPECT
+
+
+end:=

recon_test_pack/SRT2D_test_sim.par

@@ -0,0 +1,19 @@
+SRT2DParameters :=
+; test file for SRT2D
+input file := my_precorrected_sino${suffix}.hs
+
+xy output image size (in pixels) := 91
+zoom := .5
+
+;post-filter type := Wiener
+;Wiener Filter Parameters :=
+;End Wiener Filter Parameters :=
+
+;post-filter type := Gamma
+;Gamma Filter Parameters :=
+;End Gamma Filter Parameters :=
+
+output filename prefix := my_test_sim_image_SRT2D
+
+
+end:=

recon_test_pack/forward_project_SPECTUB.par

@@ -0,0 +1,35 @@
+Forward Projector parameters:=
+; example par file for specifying the forward projector for e.g. fwdtest
+type:=Matrix
+  Forward Projector Using Matrix Parameters :=
+
+     Matrix type := SPECT UB
+     Projection Matrix By Bin SPECT UB Parameters:=
+        ; same parameters as used in the OSMAPOSL_osem_SPECT.par file
+
+	psf type:= Geometrical
+
+        ;psf type:= 2D
+	;maximum number of sigmas:= 2.0
+        ;; sigma_at_depth = collimator_slope * depth_in_cm + collimator sigma 0(cm)
+        ;collimator slope := 0.0163
+        ;collimator sigma 0(cm) := 0.1466
+
+        ;Attenuation correction { Simple // Full // No }
+       	;attenuation type := No
+	;attenuation type := Simple
+        ;Values in attenuation map in cm-1 (float file)                         
+        ;attenuation map := attenuation.hv
+        attenuation map := my_atten_image_SPECT_modified.hv
+
+        ;Mask properties { Cylinder // Attenuation Map // Explicit Mask // No}
+        ;mask type := Attenuation Map
+        mask type := No
+
+     End Projection Matrix By Bin SPECT UB Parameters:=
+
+  End Forward Projector Using Matrix Parameters :=
+end:=
+
+
+

recon_test_pack/generate_atten_cylinder_SPECT.par

@@ -0,0 +1,24 @@
+generate_image Parameters :=
+output filename:=my_atten_image_SPECT
+X output image size (in pixels):=111
+Y output image size (in pixels):=111
+Z output image size (in pixels):=47
+X voxel size (in mm):= 3
+Y voxel size (in mm):= 3
+Z voxel size (in mm) :=3.27
+
+  Z number of samples to take per voxel := 1
+  Y number of samples to take per voxel := 1
+  X number of samples to take per voxel := 1
+
+shape type:= ellipsoidal cylinder
+Ellipsoidal Cylinder Parameters:=
+   radius-x (in mm):=100
+   radius-y (in mm):=100
+   length-z (in mm):=110
+   origin (in mm):={70,10,-20}
+;{70,10,-20}
+   END:=
+value := 0.096
+
+END:=

recon_test_pack/run_scatter_tests.sh

@@ -37,7 +37,7 @@ examples_dir=`stir_config --examples-dir`
 scatter_pardir="$examples_dir/samples/scatter_estimation_par_files"
 echo "Using scatter parameter files from $scatter_pardir"
 
-./simulate_PET_data_for_tests.sh
+./simulate_data_for_tests.sh
 if [ $? -ne 0 ]; then
   echo "Error running simulation"
   exit 1

recon_test_pack/run_test_simulate_and_recon.sh

@@ -9,9 +9,9 @@
 #  SPDX-License-Identifier: Apache-2.0
 #
 #  See STIR/LICENSE.txt for details
-#      
+#       
 # Author Kris Thielemans
-# 
+# Author Dimitra Kyriakopoulou
 
 echo This script should work with STIR version 6.2. If you have
 echo a later version, you might have to update your test pack.
@@ -58,27 +58,37 @@ fi
 echo "Using `command -v OSMAPOSL`"
 echo "Using `command -v OSSPS`"
 echo "Using `command -v FBP2D`"
+echo "Using `command -v FBP3DRP`"
+echo "Using `command -v SRT2D`"
+echo "Using `command -v SRT2DSPECT`"
 
 # first need to set this to the C locale, as this is what the STIR utilities use
 # otherwise, awk might interpret floating point numbers incorrectly
 LC_ALL=C
 export LC_ALL
 
-./simulate_PET_data_for_tests.sh
+./simulate_data_for_tests.sh
 if [ $? -ne 0 ]; then
   echo "Error running simulation"
   exit 1
 fi
 # need to repeat with zero-offset now as FBP doesn't support it
 zero_view_suffix=_force_zero_view_offset
-./simulate_PET_data_for_tests.sh --force_zero_view_offset --suffix $zero_view_suffix
+./simulate_data_for_tests.sh --force_zero_view_offset --suffix $zero_view_suffix
 if [ $? -ne 0 ]; then
   echo "Error running simulation with zero view offset"
   exit 1
 fi
 ## TOF data
 TOF_suffix=_TOF
-./simulate_PET_data_for_tests.sh --TOF --suffix "$TOF_suffix"
+./simulate_data_for_tests.sh --TOF --suffix "$TOF_suffix"
+if [ $? -ne 0 ]; then
+  echo "Error running simulation"
+  exit 1
+fi
+## SPECT data
+SPECT_suffix=_SPECT 
+./simulate_data_for_tests.sh --SPECT --suffix "$SPECT_suffix"
 if [ $? -ne 0 ]; then
   echo "Error running simulation"
   exit 1
@@ -96,7 +106,7 @@ input_ROI_mean=`awk 'NR>2 {print $2}' ${input_image}.roistats`
 # warning: currently OSMAPOSL needs to be run before OSSPS as 
 # the OSSPS par file uses an OSMAPOSL result as initial image
 # and reuses its subset sensitivities
-for recon in FBP2D FBP3DRP OSMAPOSL OSSPS; do
+for recon in FBP2D FBP3DRP SRT2D SRT2DSPECT OSMAPOSL OSSPS ; do  
   echo "========== Testing `command -v ${recon}`"
   # Check if we have CUDA code and parallelproj.
   # If so, check for test files in CUDA/*
@@ -118,22 +128,31 @@ for recon in FBP2D FBP3DRP OSMAPOSL OSSPS; do
     for dataSuffix in "" "$TOF_suffix"; do
       echo "===== data suffix: \"$dataSuffix\""
       # test first if analytic reconstruction and if so, run pre-correction
-      isFBP=0
+      is_analytic=0
       if expr "$recon" : FBP > /dev/null; then
-        if expr "$dataSuffix" : '.*TOF.*' > /dev/null; then
-          echo "Skipping TOF as not yet supported for FBP"
-          break
-        fi
-        isFBP=1
-        suffix=$zero_view_suffix
-        export suffix
-        echo "Running precorrection"
-        correct_projdata correct_projdata_simulation.par > my_correct_projdata_simulation.log 2>&1
-        if [ $? -ne 0 ]; then
-            echo "Error running precorrection. CHECK my_correct_projdata_simulation.log"
-            error_log_files="${error_log_files} my_correct_projdata_simulation.log"
+        is_analytic=1
+      elif expr "$recon" : SRT > /dev/null; then
+        is_analytic=1
+      fi
+      if [ $is_analytic = 1 ]; then
+          if expr "$dataSuffix" : '.*TOF.*' > /dev/null; then
+            echo "Skipping TOF as not yet supported for FBP and SRT"
             break
-        fi
+          fi
+	  if expr "$recon" : SRT2DSPECT > /dev/null; then
+	    suffix=$SPECT_suffix
+	    export suffix
+          else   
+            suffix=$zero_view_suffix
+            export suffix
+            echo "Running precorrection"
+	    correct_projdata correct_projdata_simulation.par > my_correct_projdata_simulation.log 2>&1
+	    if [ $? -ne 0 ]; then
+              echo "Error running precorrection. CHECK my_correct_projdata_simulation.log"
+	      error_log_files="${error_log_files} my_correct_projdata_simulation.log"
+	      break
+	    fi
+          fi
       else
           suffix="$dataSuffix"
           export suffix
@@ -160,7 +179,7 @@ for recon in FBP2D FBP3DRP OSMAPOSL OSSPS; do
       output_filename=`awk -F':='  '/output[ _]*filename[ _]*prefix/ { value=$2;gsub(/[ \t]/, "", value); printf("%s", value) }' "$parfile"`
       # substitute env variables (e.g. to fill in suffix)
       output_filename=`eval echo "${output_filename}"`
-      if [ ${isFBP} -eq 0 ]; then
+      if [ ${is_analytic} -eq 0 ]; then
           # iterative algorithm, so we need to append the num_subiterations
           num_subiterations=`awk -F':='  '/number[ _]*of[ _]*subiterations/ { value=$2;gsub(/[ \t]/, "", value); printf("%s", value) }' ${parfile}`
           output_filename=${output_filename}_${num_subiterations}
@@ -202,4 +221,3 @@ else
  tail ${error_log_files}
  exit 1
 fi
-