Allow TOF mashing in SSRB

documentation/release_6.2.htm

@@ -78,6 +78,11 @@ <h3>Changed functionality</h3>
     <br>
     <a href=https://github.com/UCL/STIR/pull/1458>PR #1458</a>
     </li>
+    <li>
+    SSRB now allows to mash TOF bins.
+    <br>
+    <a href=https://github.com/UCL/STIR/pull/1464>PR #1464</a>
+    </li>
 </ul>
 
 
@@ -164,6 +169,17 @@ <h4>C++ tests</h4>
   </li>
 </ul>
 
+<h4>recon_test_pack</h4>
+<ul>
+  <li>
+    The output of <tt>simulate_PET_data_for_tests.sh</tt> can now varied by setting environment variables, e.g.
+    <code>max_rd</code>.
+  </li>
+  <li>
+    New test <tt>run_test_SSRB.sh</tt>
+  </li>
+</ul>
+
 </body>
 
 </html>

recon_test_pack/README.txt

@@ -1,5 +1,5 @@
 
-README file for STIR recon_test_pack version 6.1 (and later versions?)
+README file for STIR recon_test_pack version 6.2 (and later versions?)
 ----------------------------------------------------------------------
 
 This test pack runs some simple tests to check if various STIR reconstruction
@@ -85,6 +85,9 @@ Similarly, a test for motion corrected reconstruction of gated data is run as
 
 sh run_test_simulate_and_recon_with_motion.sh  [ --mpicmd cmd] [optional_install_path]
 
+Finally, there are some tests on using SSRB to rebin the data to lower dimensions
+
+sh run_test_SSRB.sh  [ --mpicmd cmd] [optional_install_path]
 
 
 Testing SPECT reconstructions

recon_test_pack/forward_projector_proj_matrix_ray_tracing.par

@@ -4,6 +4,7 @@ type:=Matrix
   Forward Projector Using Matrix Parameters :=
   Matrix type := Ray Tracing
   Ray tracing matrix parameters :=
+   disable caching:= 1
    number of rays in tangential direction to trace for each bin := 10
   End Ray tracing matrix parameters :=
   End Forward Projector Using Matrix Parameters :=

recon_test_pack/run_test_SSRB.sh

@@ -0,0 +1,197 @@
+#! /bin/sh
+# A script to check to see if reconstruction of simulated data gives the expected result.
+#
+#  Copyright (C) 2011 - 2011-01-14, Hammersmith Imanet Ltd
+#  Copyright (C) 2011-07-01 - 2011, Kris Thielemans
+#  Copyright (C) 2014, 2022, 2024 University College London
+#  This file is part of STIR.
+#
+#  SPDX-License-Identifier: Apache-2.0
+#
+#  See STIR/LICENSE.txt for details
+#      
+# Author Kris Thielemans
+# 
+
+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.
+echo Please check the web site.
+echo
+
+#
+# Options
+#
+MPIRUN=""
+#
+# Parse option arguments (--)
+# Note that the -- is required to suppress interpretation of $1 as options 
+# to expr
+#
+while test `expr -- "$1" : "--.*"` -gt 0
+do
+
+  if test "$1" = "--mpicmd"
+  then
+    MPIRUN="$2"
+    shift 1
+  elif test "$1" = "--help"
+  then
+    echo "Usage: `basename $0` [--mpicmd somecmd] [install_dir]"
+    echo "(where [] means that an argument is optional)"
+    echo "See README.txt for more info."
+    exit 1
+  else
+    echo Warning: Unknown option "$1"
+    echo rerun with --help for more info.
+    exit 1
+  fi
+
+  shift 1
+
+done 
+
+if [ $# -eq 1 ]; then
+  echo "Prepending $1 to your PATH for the duration of this script."
+  PATH=$1:$PATH
+fi
+
+echo "Using `command -v SSRB`"
+echo "Using `command -v OSMAPOSL`"
+
+# 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
+
+echo "=== Creating single view data and test TOF rebinning ==="
+logfile=SSRB_sim1.1og
+max_rd=9 span=3 view_mash=288 tof_mash=1 ./simulate_PET_data_for_tests.sh --suffix _TOF_vm288_rd9_span3 --TOF > "$logfile" 2>&1
+if [ $? -ne 0 ]; then
+  echo "Error running simulation, Check $logfile"
+  exit 1
+fi
+logfile=SSRB_sim2.1og
+max_rd=9 span=3 view_mash=288 tof_mash=11 ./simulate_PET_data_for_tests.sh --keep_images --suffix _TOF_vm288_rd9_span3_tm11 --TOF > "$logfile" 2>&1
+if [ $? -ne 0 ]; then
+  echo "Error running simulation, Check $logfile"
+  exit 1
+fi
+logfile=SSRB_run_SSRB_sim1.log
+SSRB --template my_line_integrals_TOF_vm288_rd9_span3_tm11.hs my_line_integrals_TOF_vm288_rd9_span3_TOFSSRB11.hs my_line_integrals_TOF_vm288_rd9_span3.hs 0 > "$logfile" 2>&1
+if [ $? -ne 0 ]; then
+  echo "Error running SSRB, Check $logfile"
+  exit 1
+fi
+compare_projdata my_line_integrals_TOF_vm288_rd9_span3_tm11.hs  my_line_integrals_TOF_vm288_rd9_span3_TOFSSRB11.hs
+if [ $? -ne 0 ]; then
+  echo "Error comparing TOF mash simulation with TOF-SSRB"
+  exit 1
+fi
+
+
+echo ""
+echo "=== Creating 3D data and test axial+TOF rebinning via reconstruction ==="
+echo "We will simulate some data at span=1 and TOF-mashing=1, rebin it, then reconstruct."
+echo "The image will be checked by using an ROI calculation."
+
+echo "Running simulation (will take a while)"
+max_rd=5
+export max_rd
+view_mash=8
+export view_mash
+suffix="_TOF_vm${view_mash}_rd${max_rd}_span1"
+logfile=SSRB_sim3.1og
+span=1 tof_mash=1 background_value=.1 ./simulate_PET_data_for_tests.sh --keep_images --suffix "$suffix" --TOF > "$logfile" 2>&1
+if [ $? -ne 0 ]; then
+  echo "Error running simulation, Check $logfile"
+  exit 1
+fi
+echo "Rebinning data"
+tof_mash=11
+span=3
+suffixSSRB="_TOF_vm${view_mash}_rd${max_rd}_span${span}_TOFSSRB${tof_mash}"
+logfile="SSRB_run_SSRB_line_integrals$suffixSSRB.1og"
+SSRB "my_line_integrals${suffixSSRB}.hs" "my_line_integrals${suffix}.hs" "$span" 1 1 "$max_rd" "$tof_mash" > "$logfile" 2>&1
+if [ $? -ne 0 ]; then
+  echo "Error running SSRB, Check $logfile"
+  exit 1
+fi
+logfile="SSRB_run_SSRB_acfs$suffixSSRB.1og"
+SSRB "my_acfs${suffixSSRB}.hs" "my_acfs${suffix}.hs" "$span" 1 1 "$max_rd" "$tof_mash" > "$logfile" 2>&1
+if [ $? -ne 0 ]; then
+  echo "Error running SSRB, Check $logfile"
+  exit 1
+fi
+logfile="SSRB_run_SSRB_prompts$suffixSSRB.1og"
+SSRB "my_prompts${suffixSSRB}.hs" "my_prompts${suffix}.hs" "$span" 1 1 "$max_rd" "$tof_mash" > "$logfile" 2>&1
+if [ $? -ne 0 ]; then
+  echo "Error running SSRB, Check $logfile"
+  exit 1
+fi
+logfile="SSRB_run_SSRB_additive_sinograms$suffixSSRB.1og"
+SSRB "my_additive_sinogram${suffixSSRB}.hs" "my_additive_sinogram${suffix}.hs" "$span" 1 1 "$max_rd" "$tof_mash" > "$logfile" 2>&1
+if [ $? -ne 0 ]; then
+  echo "Error running SSRB, Check $logfile"
+  exit 1
+fi
+
+
+# Below code is a copy of relevant lines in run_test_simulate_and_recon.sh.
+# It would be far better to make those into functions/scripts...
+
+recon=OSMAPOSL
+isFBP=0
+parfile=${recon}_test_sim_PM.par
+suffix=${suffixSSRB}
+export suffix
+
+num_subsets=2
+export num_subsets
+
+input_image=my_uniform_cylinder.hv
+input_voxel_size_x=`stir_print_voxel_sizes.sh ${input_image}|awk '{print $3}'`
+ROI=ROI_uniform_cylinder.par
+list_ROI_values ${input_image}.roistats ${input_image} ${ROI} 0 > /dev/null 2>&1
+input_ROI_mean=`awk 'NR>2 {print $2}' ${input_image}.roistats`
+
+# Indentation is here because straight copy of the file.
+
+      # run actual reconstruction
+      echo "Running ${recon} ${parfile}"
+      logfile="my_${parfile}${suffix}.log"
+      ${MPIRUN} ${recon} ${parfile} > "$logfile" 2>&1
+      if [ $? -ne 0 ]; then
+          echo "Error running reconstruction. CHECK RECONSTRUCTION LOG \"$logfile\""
+          exit 1
+      fi
+
+      # find filename of (last) image from ${parfile}
+      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
+          # 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}
+      fi
+      output_image=${output_filename}.hv
+
+      # compute ROI value
+      list_ROI_values ${output_image}.roistats ${output_image} ${ROI} 0  > ${output_image}.roistats.log 2>&1
+      if [ $? -ne 0 ]; then
+          echo "Error running list_ROI_values. CHECK LOG ${output_image}.roistats.log"
+          exit 1
+      fi
+
+      # compare ROI value
+      output_voxel_size_x=`stir_print_voxel_sizes.sh ${output_image}|awk '{print $3}'`
+      output_ROI_mean=`awk "NR>2 {print \\$2*${input_voxel_size_x}/${output_voxel_size_x}}" ${output_image}.roistats`
+      echo "Input ROI mean: $input_ROI_mean"
+      echo "Output ROI mean: $output_ROI_mean"
+      error_bigger_than_1percent=`echo $input_ROI_mean $output_ROI_mean| awk '{ print(($2/$1 - 1)*($2/$1 - 1)>0.0001) }'`
+      if [ ${error_bigger_than_1percent} -eq 1 ]; then
+          echo "DIFFERENCE IN ROI VALUES IS TOO LARGE. CHECK RECONSTRUCTION LOG "$logfile""
+      else
+          echo "This seems fine."
+      fi
+

recon_test_pack/simulate_PET_data_for_tests.sh

@@ -1,13 +1,13 @@
-#! /bin/sh
+#! /bin/sh -vx
 # A script to simulate some data used by other tests.
 # Careful: run_scatter_tests.sh compares with previously generated data
-# so you cannot simply modify this one with adjusting that as well.
+# so you cannot simply modify this one without adjusting that as well.
 #
 # This script is not intended to be run on its own!
 #
 #  Copyright (C) 2011 - 2011-01-14, Hammersmith Imanet Ltd
 #  Copyright (C) 2011-07-01 - 2011, Kris Thielemans
-#  Copyright (C) 2014, 2020, 2022 University College London
+#  Copyright (C) 2014, 2020, 2022, 2024 University College London
 #  This file is part of STIR.
 #
 #  SPDX-License-Identifier: Apache-2.0
@@ -25,6 +25,7 @@ echo
 command -v generate_image >/dev/null 2>&1 || { echo "generate_image not found or not executable. Aborting." >&2; exit 1; }
 echo "Using `command -v generate_image`"
 
+generate_images=1
 force_zero_view_offset=0
 TOF=0
 suffix=""
@@ -46,10 +47,19 @@ do
   then
     suffix="$2"
     shift 1
+  elif test "$1" = "--keep_images"
+  then
+    generate_images=0
   elif test "$1" = "--help"
   then
-    echo "Usage: `basename $0` [--force_zero_view_offset]  [--suffix sometext] [install_dir]"
+    echo "Usage: `basename $0` [--force_zero_view_offset]  [--suffix sometext] [--keep_images] [install_dir]"
     echo "(where [] means that an argument is optional)"
+    echo "The suffix is appended to filenames (before adding the extensions .hs/.s)."
+    echo "--keep-images can be used to avoid regenerating the images (if you know they are correct)."
+    echo "You can choose data sizes different from the defaults by setting the following environment variables:"
+    echo "   view_mash, span, max_rd, tof_mash"
+    echo "Randoms+scatter are set to a constant proj_data. you can set that value via the environment variable:"
+    echo "   background_value"
     exit 1
   else
     echo Warning: Unknown option "$1"
@@ -71,38 +81,46 @@ fi
 LC_ALL=C
 export LC_ALL
 
+if [ "$generate_images" -eq 1 ]; then
+  echo "===  make emission image"
+  generate_image  generate_uniform_cylinder.par
+  echo "===  make attenuation image"
+  generate_image  generate_atten_cylinder.par
+fi
 
-echo "===  make emission image"
-generate_image  generate_uniform_cylinder.par
-echo "===  make attenuation image"
-generate_image  generate_atten_cylinder.par
 if [ "$TOF" -eq 0 ]; then
-  echo "===  create template sinogram (DSTE in 3D with max ring diff 2 to save time)"
-  template_sino=my_DSTE_3D_rd3_template.hs
+  : ${view_mash:=1}
+  : ${span:=2}
+  : ${max_rd:=3}
+  echo "===  create template sinogram (DSTE with view_mash=${view_mash}, max_ring_diff=${max_rd})"
+  template_sino=my_DSTE_3D_vm${view_mash}_span${span}_rd${max_rd}_template.hs
   cat > my_input.txt <<EOF
 Discovery STE
 
-1
+$view_mash
 n
 
-0
-2
+$span
+$max_rd
 EOF
 else
-  echo "===  create template sinogram (D690 in 3D with view-mash =2, TOF-mash=11, max ring diff 3 to save time)"
-  template_sino=my_D690_3D_rd2_template.hs
+  : ${view_mash:=2}
+  : ${span:=2}
+  : ${max_rd:=3}
+  : ${tof_mash:=11}
+  echo "===  create template sinogram (D690 with view_mash=${view_mash}, TOF_mash=${tof_mash}, max_ring_diff ${max_rd})"
+  template_sino=my_D690_3D_vm${view_mash}_span${span}_rd${max_rd}_TOFmash${tof_mash}_template.hs
   cat > my_input.txt <<EOF
 Discovery 690
 
-2
-11
+$view_mash
+$tof_mash
 N
 
-2
-3
+$span
+$max_rd
 EOF
 fi
-
 create_projdata_template  ${template_sino} < my_input.txt > my_create_${template_sino}.log 2>&1
 if [ $? -ne 0 ]; then 
   echo "ERROR running create_projdata_template. Check my_create_${template_sino}.log"; exit 1; 
@@ -125,7 +143,8 @@ if [ $force_zero_view_offset -eq 1 ]; then
 fi
 
 # create sinograms
-./simulate_data.sh my_uniform_cylinder.hv my_atten_image.hv ${template_sino} 10 ${suffix}
+: ${background_value:=10}
+./simulate_data.sh my_uniform_cylinder.hv my_atten_image.hv ${template_sino} ${background_value} ${suffix}
 if [ $? -ne 0 ]; then
   echo "Error running simulation"
   exit 1