Merge pull request #181 from pwadhwa351/ViewOffset

View offset for PET data finally taken into account. WARNING: this PR changes image orientation when view-offset!=0 or when view-mashing is used.

CMakeLists.txt

@@ -150,8 +150,9 @@ if(NOT DISABLE_NLOHMANN_JSON)
     endif()
 endif()
 
-# Enable version 3 projectors
+# Legacy options
 option(STIR_PROJECTORS_AS_V3 "Enable STIR version 3 legacy code" OFF)
+option(STIR_LEGACY_IGNORE_VIEW_OFFSET "Ignore using scanner view-offset (or intrinsic azimuthal tilt), as in STIR 4 or earlier" OFF)
 
 # NiftyPET projector
 if(NOT DISABLE_NiftyPET_PROJECTOR OR NOT DISABLE_Parallelproj_PROJECTOR)

documentation/release_5.0.htm

@@ -50,8 +50,18 @@ <h2> Summary for end users (also to be read by developers)</h2>
 
 <h3>Changes breaking backwards compatibility from a user-perspective</h3>
 <ul>
-  <li> </li>
-  <ul>
+<li>
+View Offset Support enabled for the PET scanners.
+<ul>
+ <li>For the scanners that have non-zero intrinsic azimuthal tilt angle, reconstructed images will now get rotated.</li>
+ <li>If you use view mashing, in previous STIR versions images were rotated according to half the number of mashed views. This is now corrected.</li>
+</ul>
+<strong>WARNING</strong> This means that reconstructed PET images will <strong>not be identical</strong> when either the scanner has non-zero view-offset, or view-mashing is on.<br />
+To reflect this change, <code>Scanner::get_default_intrinsic_tilt()</code> has been renamed to <code>get_intrinsic_azimuthal_tilt()</code>.<br />
+Note: start angle was already supported for SPECT<br />
+Backward compatibility for reconstructed images can be achieved by setting the <tt>STIR_LEGACY_IGNORE_VIEW_OFFSET</tt> CMake option to <tt>ON</tt>. (However, copied sinogram data will then always have the offset set to 0, in contrast to earlier versions of STIR).
+</li>
+</ul>
 
 
 <h3>Bug fixes</h3>
@@ -115,6 +125,8 @@ <h3>New functionality</h3>
 <li>
 	SAFIR input file format now supports a Gaussian LOR randomization, which is applied on the endpoints of each LOR when sorting the listmode events into virtual scanner projection data.
 </li>
+<li><code>Scanner</code> has now a new <code>static</code> member <code>get_names_of_predefined_scanners</code> returning a list of names.
+</li>
 </ul>
 
 
@@ -245,6 +257,9 @@ <h3>Other code changes</h3>
 <ul>
     <li>store data in <code>ProjDataInMemory</code> in the same order as what is used by <code>copy_to</code> and <code>fill_from</code>. This enabled using straight-forward copy. (This change should not affect anyone, except if you relied on a specific order in the buffer before.)
     </li>
+<li><code>LOR</code> classes now no longer require phi at input to be between 0 and pi,
+  nor psi to be between 0 and 2 pi. They will standardise the input automatically.
+  </li>
 </ul>
 
 </body>

examples/python/print_scanner_info.py

@@ -0,0 +1,32 @@
+# Demo of how to use STIR from python to list some scanner properties
+
+# Copyright 2021 University College London
+# Author Kris Thielemans
+
+# This file is part of STIR.
+#
+# This file is free software; you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation; either version 2.1 of the License, or
+# (at your option) any later version.
+#
+# This file is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# See STIR/LICENSE.txt for details
+
+#%% Initial imports
+import stir
+import math
+
+#%% get the list of all scanners
+l=stir.Scanner.get_names_of_predefined_scanners()
+
+#%% go through the list to print some properties
+for s in l:
+  sc=stir.Scanner.get_scanner_from_name(s)
+  if (sc.get_intrinsic_azimuthal_tilt() != 0):
+    print("{0:30}: intrinsic tilt (degrees): {1:3.2f}".format(sc.get_name(), sc.get_intrinsic_azimuthal_tilt()*180/math.pi))
+

recon_test_pack/FBP2D_test_sim.par

@@ -1,6 +1,6 @@
 FBP2DParameters :=
 ; test file for FBP2D
-input file := my_precorrected_sino.hs
+input file := my_precorrected_sino${suffix}.hs
 xy output image size (in pixels) := 91
 zoom := .5
 

recon_test_pack/FBP3DRP_test_sim.par

@@ -1,6 +1,6 @@
 FBP3DRPParameters :=
 ; test file for FBP2D
-input file := my_precorrected_sino.hs
+input file := my_precorrected_sino${suffix}.hs
 xy output image size (in pixels) := 91
 zoom := .5
 

recon_test_pack/Utahscat600k_ca_seg4.hs

@@ -26,6 +26,10 @@ maximum ring difference per segment := { -4,-3,-2,-1,0,1,2,3,4}
 number of rings := 16
 number of detectors per ring := 384
 distance between rings (cm) := 0.675
+; This offset is constructed such that it compensates for the offset implied
+; in the view-mashing. This means that the first view will be "vertical"
+; as it was in older versions of STIR
+view offset (degrees) := -0.46875
 number of time frames := 1
 image duration (sec) [1]:= 60
 image relative start time (sec)[1] := 100

recon_test_pack/correct_projdata_simulation.par

@@ -1,19 +1,19 @@
 correct_projdata Parameters := 
 
-  input file := my_prompts.hs
-  output filename := my_precorrected_sino.hs
+  input file := my_prompts${suffix}.hs
+  output filename := my_precorrected_sino${suffix}.hs
   ; use data (1) or set to one (0) := 
   apply (1) or undo (0) correction := 1
-  randoms projdata filename := my_randoms.hs
+  randoms projdata filename := my_randoms${suffix}.hs
 
   Bin Normalisation type := from projdata
     Bin Normalisation From ProjData :=
     ; only attenuation here for this simulation
-    normalisation projdata filename:= my_acfs.hs
+    normalisation projdata filename:= my_acfs${suffix}.hs
     End Bin Normalisation From ProjData:=
 
   ; scatter term to be subtracted AFTER norm+atten correction
   ; defaults to 0
-  ;scatter projdata filename := scatter.hs
+  ;scatter projdata filename := scatter${suffix}.hs
 
 END:= 

recon_test_pack/run_test_simulate_and_recon.sh

@@ -82,6 +82,13 @@ 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
+if [ $? -ne 0 ]; then
+  echo "Error running simulation with zero view offset"
+  exit 1
+fi
 
 error_log_files=""
 
@@ -103,13 +110,18 @@ for recon in FBP2D FBP3DRP OSMAPOSL OSSPS; do
     isFBP=0
     if expr ${recon} : FBP > /dev/null; then
       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"
         break
       fi
+    else
+      suffix=""
+      export suffix
     fi
 
     # run actual reconstruction

recon_test_pack/run_tests_modelling.sh

@@ -153,17 +153,9 @@ done
 # Create the appropriate proj_data files
 extract_single_images_from_dynamic_image dyn_from_p0005-p5_img_f%dg1d0b0.hv dyn_from_p0005-p5.hv
 # if [ ! -r fwd_dyn_from_p0005-p5.S ]; then
-rm -f fwd.par
-cat <<EOF > fwd.par
-Forward Projector parameters:=
-type :=  ray tracing
-  forward  Projector Using Ray Tracing Parameters :=
-  End Forward Projector Using Ray Tracing Parameters :=
-end :=
-EOF
 
 for fr in `count 23 28`; do
-    forward_project fwd_dyn_from_p0005-p5_f${fr}g1d0b0  dyn_from_p0005-p5_img_f${fr}g1d0b0.hv ${INPUTDIR}ECAT_931_projdata_template.hs fwd.par
+    forward_project fwd_dyn_from_p0005-p5_f${fr}g1d0b0  dyn_from_p0005-p5_img_f${fr}g1d0b0.hv ${INPUTDIR}ECAT_931_projdata_template.hs
 done
 #fi
 

recon_test_pack/simulate_PET_data_for_tests.sh

@@ -7,7 +7,7 @@
 #
 #  Copyright (C) 2011 - 2011-01-14, Hammersmith Imanet Ltd
 #  Copyright (C) 2011-07-01 - 2011, Kris Thielemans
-#  Copyright (C) 2014, University College London
+#  Copyright (C) 2014,2020 University College London
 #  This file is part of STIR.
 #
 #  This file is free software; you can redistribute it and/or modify
@@ -25,14 +25,51 @@
 # Author Kris Thielemans
 # 
 
-echo This script should work with STIR version 3.x, 4.x. If you have
+echo This script should work with STIR version 5.x. If you have
 echo a later version, you might have to update your test pack.
 echo Please check the web site.
 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`"
 
+force_zero_view_offset=0
+suffix=""
+#
+# 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" = "--force_zero_view_offset"
+  then
+    force_zero_view_offset=1
+  elif test "$1" = "--suffix"
+  then
+    suffix="$2"
+    shift 1
+  elif test "$1" = "--help"
+  then
+    echo "Usage: `basename $0` [--force_zero_view_offset]  [--suffix sometext] [install_dir]"
+    echo "(where [] means that an argument is optional)"
+    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
+
 # 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
@@ -64,8 +101,14 @@ awk '/END OF INTERFILE/ { print "number of energy windows := 1\nenergy window lo
     ${template_sino} > tmp_header.hs
 mv tmp_header.hs ${template_sino}
 
+if [ $force_zero_view_offset -eq 1 ]; then
+  new_template_sino=my_DSTE_3D_rd2_template$suffix.hs
+  force_view_offset_to_zero.sh ${new_template_sino} ${template_sino}
+  template_sino=${new_template_sino}
+fi
+
 # create sinograms
-./simulate_data.sh my_uniform_cylinder.hv my_atten_image.hv ${template_sino}
+./simulate_data.sh my_uniform_cylinder.hv my_atten_image.hv ${template_sino} 10 ${suffix}
 if [ $? -ne 0 ]; then
   echo "Error running simulation"
   exit 1

recon_test_pack/simulate_data.sh

@@ -20,9 +20,9 @@
 # Author Kris Thielemans
 # 
 
-if [ $# -lt 3  -o $# -gt 4 ]; then
-  echo "Usage: `basename $0` emission_image attenuation_image template_sino  [ background_value ]"
-  echo "Creates my_prompts.hs my_randoms.hs my_acfs.hs my_additive_sinogram.hs"
+if [ $# -lt 3  -o $# -gt 5 ]; then
+  echo "Usage: `basename $0` emission_image attenuation_image template_sino  [ background_value [suffix] ]"
+  echo "Creates my_prompts$suffix.hs my_randoms$suffix.hs my_acfs$suffix.hs my_additive_sinogram$suffix.hs"
   exit 1
 fi
 
@@ -34,15 +34,17 @@ if [ $# -gt 3 ]; then
 else
   background_value=10
 fi
-
+if [ $# -gt 4 ]; then
+  suffix=$5
+fi
 echo "===  create ACFs"
-calculate_attenuation_coefficients --ACF my_acfs.hs ${atten_image} ${template_sino} > my_create_acfs.log 2>&1
+calculate_attenuation_coefficients --ACF my_acfs$suffix.hs ${atten_image} ${template_sino} > my_create_acfs.log 2>&1
 if [ $? -ne 0 ]; then 
   echo "ERROR running calculate_attenuation_coefficients. Check my_create_acfs.log"; exit 1; 
 fi
 
 echo "===  create line integrals"
-forward_project my_line_integrals.hs  ${emission_image} ${template_sino} forward_projector_proj_matrix_ray_tracing.par > my_create_line_integrals.log 2>&1
+forward_project my_line_integrals$suffix.hs  ${emission_image} ${template_sino} forward_projector_proj_matrix_ray_tracing.par > my_create_line_integrals.log 2>&1
 if [ $? -ne 0 ]; then 
   echo "ERROR running forward_project. Check my_create_line_integrals.log"; exit 1; 
 fi
@@ -51,17 +53,18 @@ fi
 echo "=== create constant randoms background"
 ${INSTALL_DIR}stir_math -s --including-first \
          --times-scalar 0 --add-scalar $background_value \
-         my_randoms my_line_integrals.hs
+         my_randoms$suffix my_line_integrals$suffix.hs
 if [ $? -ne 0 ]; then 
   echo "ERROR running stir_math"; exit 1; 
 fi
 
 echo "===  create norm factors"
 # currently just 1 as not used in rest of script yet.
 stir_math -s --including-first \
-          --times-scalar 0 --add-scalar 1 my_norm.hs my_line_integrals.hs
+          --times-scalar 0 --add-scalar 1 my_norm$suffix.hs my_line_integrals$suffix.hs
 
 echo "===  create prompts"
+export suffix # used in the .par file to determine filenames
 correct_projdata uncorrect_projdata_simulation.par > my_create_prompts.log 2>&1
 if [ $? -ne 0 ]; then 
   echo "ERROR running correct_projdata. Check my_create_prompts.log"; exit 1; 
@@ -71,7 +74,7 @@ fi
 
 echo "===  create additive sinogram for reconstruction"
 # need randoms (and scatter) multiplied by ACF and norm (but we don't have a norm here)
-stir_math -s --mult  my_additive_sinogram.hs my_randoms.hs my_acfs.hs
+stir_math -s --mult  my_additive_sinogram$suffix.hs my_randoms$suffix.hs my_acfs$suffix.hs
 if [ $? -ne 0 ]; then 
   echo "ERROR running stir_math"; exit 1; 
 fi

recon_test_pack/uncorrect_projdata_simulation.par

@@ -1,19 +1,19 @@
 correct_projdata Parameters := 
 
-  input file := my_line_integrals.hs
-  output filename := my_prompts.hs
+  input file := my_line_integrals${suffix}.hs
+  output filename := my_prompts${suffix}.hs
   ; use data (1) or set to one (0) := 
   apply (1) or undo (0) correction := 0
-  randoms projdata filename := my_randoms.hs
+  randoms projdata filename := my_randoms${suffix}.hs
 
   Bin Normalisation type := from projdata
     Bin Normalisation From ProjData :=
     ; only attenuation here for this simulation
-    normalisation projdata filename:= my_acfs.hs
+    normalisation projdata filename:= my_acfs${suffix}.hs
     End Bin Normalisation From ProjData:=
 
   ; scatter term to be subtracted AFTER norm+atten correction
   ; defaults to 0
-  ;scatter projdata filename := scatter.hs
+  ;scatter projdata filename := scatter${suffix}.hs
 
 END:= 

scripts/CMakeLists.txt

@@ -30,7 +30,8 @@ set(${dir_scripts}
   zoom_att_image.sh
   get_num_voxels.sh
   create_fdef_from_listmode.sh
-)
+  force_view_offset_to_zero.sh
+  )
 
 
 

scripts/force_view_offset_to_zero.sh

@@ -0,0 +1,33 @@
+#!/bin/sh
+# Takes an interfile projdata header, and creates a new one with view-offset set to 0.
+#
+#
+#  Copyright (C) 2020, University College London
+#  This file is part of STIR.
+#
+#  This file is free software; you can redistribute it and/or modify
+#  it under the terms of the GNU Lesser General Public License as published by
+#  the Free Software Foundation; either version 2.1 of the License, or
+#  (at your option) any later version.
+
+#  This file is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU Lesser General Public License for more details.
+#
+#  See STIR/LICENSE.txt for details
+#
+# Author: Kris Thielemans
+
+if [ $# -ne 2 ]; then
+  echo "Usage: `basename $0` output_filename projdata_interfile_header_filename"
+  echo "sets view-offset to zero"
+  echo "Warning: this script is unsafe. It doesn't do any checks at all."
+  exit 1
+fi
+
+out_filename=$1
+in_filename=$2
+
+set -e
+sed -e "s/View offset (degrees).*/View offset (degrees) := 0/" $in_filename > $out_filename