Time Of Flight reconstruction

CMakeLists.txt

@@ -116,6 +116,10 @@ if(NOT DISABLE_CERN_ROOT)
   endif()
 endif()
 
+if(NOT DISABLE_HDF5_SUPPORT)
+    find_package(HDF5 COMPONENTS CXX)
+endif()
+
 if(NOT DISABLE_AVW)
   find_package(AVW)
 endif()

examples/PET_simulation/generate_input_data.sh

@@ -41,6 +41,7 @@ template_sino=my_DSTE_3D_rd1_template.hs
 cat > my_input.txt <<EOF
 Discovery STE
 1
+0
 n
 
 0

recon_test_pack/OSMAPOSL_test_lmf.par

@@ -0,0 +1,38 @@
+OSMAPOSLParameters := 
+objective function type:=  PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin
+PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin Parameters:= 
+	list mode filename := PET_ACQ_small.l.hdr.STIR
+	max ring difference num to process :=  2
+	projector pair type := Matrix
+  Projector Pair Using Matrix Parameters :=
+  Matrix type := Ray Tracing
+  Ray tracing matrix parameters :=
+   ; use multiple (almost) parallel LORs for every bin in the sinogram 
+   ; to avoid discretisation artefacts
+   number of rays in tangential direction to trace for each bin:= 3
+   ; you could disable some symmetries if you have enough memory
+   ; this would for instance allow you to increase the number of subsets
+   ; do_symmetry_90degrees_min_phi:=0
+  End Ray tracing matrix parameters :=
+  End Projector Pair Using Matrix Parameters :=
+
+        Bin Normalisation type := From ProjData
+	  Bin Normalisation From ProjData :=
+	   normalisation projdata filename:= my_acfs.hs
+	End Bin Normalisation From ProjData:=
+
+	;num_events_to_store := 100
+	recompute sensitivity :=1
+	use subset sensitivities:= 0
+	sensitivity filename:=  my_sens_t_lm_pr_seg2.hv
+	zoom :=  1
+
+additive sinogram := my_MLrandoms_f1.hs
+
+end PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin Parameters:= 
+enforce initial positivity condition:=  1
+number of subsets:= 1
+number of subiterations:= 1
+save estimates at subiteration intervals:= 1
+output filename prefix := my_output_t_lm_pr_seg2
+END := 

recon_test_pack/Siemens_mMR_seg2.hs

@@ -44,4 +44,6 @@ Number of crystals per singles unit in transaxial direction := 9
 end scanner parameters:=
 effective central bin size (cm) := 0.208815
 number of time frames := 1
+start vertical bed position (mm) := 0
+start horizontal bed position (mm) := 0
 !END OF INTERFILE :=

recon_test_pack/lm_generate_atten_cylinder.par

@@ -1,8 +1,8 @@
 generate_image Parameters :=
 output filename:=my_atten_image
-X output image size (in pixels):=111
-Y output image size (in pixels):=111
-Z output image size (in pixels):=65
+X output image size (in pixels):=100
+Y output image size (in pixels):=100
+Z output image size (in pixels):=127
 X voxel size (in mm):= 3
 Y voxel size (in mm):= 3
 Z voxel size (in mm) := 0.40625
@@ -16,7 +16,7 @@ Ellipsoidal Cylinder Parameters:=
    radius-x (in mm):=100
    radius-y (in mm):=100
    length-z (in mm):=110
-   origin (in mm):={70,10,20}
+   origin (in mm):={70,0,0}
    END:=
 value := 0.096
 

recon_test_pack/root_header.hroot

@@ -11,7 +11,11 @@ Default bin size (cm)                    := 0.208626
 Maximum number of non-arc-corrected bins := 344
 Default number of arc-corrected bins 	 := 344
 View offset (degrees)					 := 0
+Number of TOF time bins := 410
+Size of timing bin (ps) := 10.00
+Timing resolution (ps) := 400.0
 
+%TOF mashing factor:= 82
 
 GATE scanner type := GATE_Cylindrical_PET
 GATE_Cylindrical_PET Parameters :=

recon_test_pack/run_test_listmode_recon.sh

@@ -70,7 +70,7 @@ ThereWereErrors=0
 echo "=== Simulate normalisation data"
 # For normalisation data we are going to use a cylinder in the center,
 # with water attenuation values
-echo "=== Gnerete fake emission image"
+echo "=== Generate fake emission image"
 generate_image  lm_generate_atten_cylinder.par
 echo "=== Calculate ACFs"
 calculate_attenuation_coefficients --ACF my_acfs.hs my_atten_image.hv Siemens_mMR_seg2.hs > my_create_acfs.log 2>&1

recon_test_pack/run_test_simulate_and_recon.sh

@@ -70,6 +70,7 @@ LC_ALL=C
 export LC_ALL
 
 ./simulate_PET_data_for_tests.sh
+
 if [ $? -ne 0 ]; then
   echo "Error running simulation"
   exit 1

recon_test_pack/run_test_simulate_and_recon_with_motion.sh

@@ -116,7 +116,7 @@ fi
 echo "===  create template sinogram (DSTE in 3D with max ring diff 2 to save time)"
 template_sino=my_DSTE_3D_rd2_template.hs
 cat > my_input.txt <<EOF
-Discovery STE
+Discovery STE nonTOF
 1
 n
 

recon_test_pack/run_test_time_of_flight.sh

@@ -0,0 +1,106 @@
+#! /bin/sh
+# A script to check to see if Time Of Flight data are binned and used properly
+#
+#  Copyright (C) 2016, University of Leeds
+#  Copyright (C) 2017, University of Hull
+#  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 Nikos Efthimiou
+# 
+
+# Scripts should exit with error code when a test fails:
+if [ -n "$TRAVIS" ]; then
+    # The code runs inside Travis
+    set -e
+fi
+
+echo This script should work with STIR version '>'3.0. If you have
+echo a later version, you might have to update your test pack.
+echo Please check the web site.
+echo
+
+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
+export LC_ALL
+
+echo "===  create template sinogram. We'll use a test_scanner which is small and 
+has TOF info"
+template_sino=my_test_scanner_template.hs
+cat > my_input.txt <<EOF
+test_scanner
+1
+82
+N
+
+
+
+EOF
+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; 
+fi
+
+export INPUT_ROOT_FILE=test_PET_GATE.root
+export EXCLUDE_RANDOM=1
+export EXCLUDE_SCATTERED=1
+
+INPUT=root_header.hroot TEMPLATE=$template_sino OUT_PROJDATA_FILE=my_tof_sinogram lm_to_projdata  --test_timing_positions lm_to_projdata.par > my_write_TOF_values_${template_sino}.log 2>&1
+
+if [ $? -ne 0 ]; then 
+  echo "ERROR running lm_to_projdata  --test_timing_positions. Check my_write_TOF_values_${template_sino}.log"; exit 1; 
+fi
+
+echo "Comparing values in TOF sinogram ..."
+list_projdata_info --all my_tof_sinogram_f179g1d0b0.hs > my_sino_values_$template_sino.log 2>&1
+if [ $? -ne 0 ]; then 
+  echo "ERROR running list_projdata_info. Check my_sino_values_$template_sino.log"; 
+  exit 1; 
+fi
+
+
+TOF_bins=$(grep 'Total number of timing positions' my_sino_values_$template_sino.log | awk -F ':' '{ print  $2 }')
+echo "Total number of TOF bins:" $TOF_bins
+
+Timming_Locations=$(grep 'Timing location' my_sino_values_$template_sino.log | awk -F ':' '{ print  $2 }')
+echo "Timming_Locations:" $Timming_Locations
+
+Data_mins=$(grep 'Data min' my_sino_values_$template_sino.log | awk -F ':' '{ print  $2 }')
+echo "Data mins:" $Data_mins
+
+Data_maxs=$(grep 'Data max' my_sino_values_$template_sino.log | awk -F ':' '{ print  $2 }')
+echo "Data maxs:" $Data_maxs
+
+for i in $(seq 5)
+do
+  if [ $(( $((i-1 - TOF_bins/2)) - Data_mins[i])) -ne 0 ]; then
+    echo "Wrong values in TOF sinogram. Error. $(( $(($((i-1)) -TOF_bins/2)) - Data_mins[i]))"
+    exit 1
+    fi
+done
+
+
+echo
+echo '--------------- End of Time-Of-Flight tests -------------'
+echo
+echo "Everything seems to be fine !"
+echo 'You could remove all generated files using "rm -f my_* *.log"'
+exit 0
+

recon_test_pack/run_test_zoom_image.sh

@@ -76,7 +76,7 @@ get_ROI_value() {
 # warning: return 0 if they are different (which is non-standard I guess)
 compare_values() {
     if [ $# -ne 3 ]; then
-      echo "Something wrong with call to compare_values"
+      echo "Someting wrong with call to compare_values"
     exit 1
     fi
     error_bigger_than_x=`echo $1 $2 $3 | awk '{ print(($2/$1 - 1)*($2/$1 - 1)> ($3 * $3)) }'`
@@ -148,6 +148,7 @@ template_sino=my_DSTE_3D_rd2_template.hs
 cat > my_input.txt <<EOF
 Discovery STE
 1
+410
 n
 
 0
@@ -167,19 +168,19 @@ if [ $? -ne 0 ]; then
   exit 1
 fi
 
-org_sum=`list_projdata_info --sum my_prompts.hs | awk -F: '{ print $2}'`
+org_sum=`list_projdata_info --sum my_prompts.hs | awk -F: '/sum/ { print $2}'`
 
 ./simulate_data.sh my_zoom_test4.hv my_atten_image.hv ${template_sino} 0
 if [ $? -ne 0 ]; then
   echo "Error running simulation"
   exit 1
 fi
 
-new_sum=`list_projdata_info --sum my_prompts.hs | awk -F: '{ print $2}'`
+new_sum=`list_projdata_info --sum my_prompts.hs | awk -F: '/sum/ { print $2}'`
 
 if compare_values $org_sum $new_sum .01
 then
-  echo "DIFFERENCE IN su  of prompts IS TOO LARGE."
+  echo "DIFFERENCE IN sum of prompts IS TOO LARGE."
   exit 1
 fi
 

recon_test_pack/simulate_PET_data_for_tests.sh

@@ -77,6 +77,7 @@ template_sino=my_DSTE_3D_rd2_template.hs
 cat > my_input.txt <<EOF
 Discovery STE
 1
+0
 n
 
 0

recon_test_pack/template_for_ROOT_scanner.hs

@@ -34,6 +34,11 @@ Default bin size (cm)                    := 0.208626
 View offset (degrees)                    := 0
 Maximum number of non-arc-corrected bins := 344
 Default number of arc-corrected bins     := 344
+Energy resolution := 0
+Reference energy (in keV) := 511
+Number of TOF time bins :=410
+Size of timing bin (ps) := 10
+Timing resolution (ps) := 400
 Number of blocks per bucket in transaxial direction         := 1
 Number of blocks per bucket in axial direction              := 1
 Number of crystals per block in axial direction             := 4

scripts/plot_TOF_bins.m

@@ -0,0 +1,70 @@
+%%
+%Plot TOF bins. 
+% Nikos Efthimiou. 2018/11/01
+% University of Hull 
+
+%This scripts loads LOR files, exported by test_time_of_flight to the disk and plots them. 
+%%
+clc; clear all;
+%Path to TOF files. 
+path_name ='/home/nikos/Desktop/conv_LOR/'
+
+pre_sort_files_in_path = dir(path_name)
+nums = []
+names = []
+
+for i = 1: size(pre_sort_files_in_path)
+    cur_file = pre_sort_files_in_path(i).name   
+   if strfind (cur_file, 'glor')
+       num = sscanf(cur_file,'glor_%d')
+       % The following number can change accordingly.
+       if ((mod(num,1)==0) || num == 500000000)
+        nums{end+1} = int32(num);
+        names{end+1} = cur_file;
+      end
+
+   end   
+end
+
+clear cur_file
+sorted_filenames = cell(numel(nums),2);
+[Sorted_A, Index_A] = sort(cell2mat(nums));
+sorted_filenames(:,2) = names(Index_A);
+
+% hold x values
+x_values = [];
+% hold the tof bins.
+y_tf_values = [];
+% hold the non tof LOR
+y__ntf_values = [];
+
+for i = 1 : size(sorted_filenames,1)
+    cur_file = sorted_filenames{i,2};
+
+    if strfind (cur_file, 'glor')
+
+        if strfind(cur_file, '500000000')
+            cur_full_path = fullfile(path_name, cur_file);
+
+            A = importdata(cur_full_path);
+            y_ntf_values = A(:,2);
+        else
+            cur_full_path = fullfile(path_name, cur_file);
+
+            A = importdata(cur_full_path);
+
+            if size(x_values) == 0
+                x_values = A(:,1);
+            end
+
+            y_tf_values = [y_tf_values A(:,2)];
+
+        end
+    end
+end
+
+sum_of_all_bins = sum(y_tf_values,2);
+x_v = x_values/0.299; 
+
+%% Create Plot
+plot(x_v,y_tf_values(:,:), x_v, sum_of_all_bins, x_v, y_ntf_values)

src/CMakeLists.txt

@@ -67,7 +67,6 @@ else()
   message(STATUS "HDF5 support disabled.")
 endif()
 
-
 if (ITK_FOUND) 
   message(STATUS "ITK libraries added.")
   set(HAVE_ITK ON)

src/IO/GEHDF5Wrapper.cxx

@@ -344,7 +344,6 @@ shared_ptr<Scanner> GEHDF5Wrapper::get_scanner_from_HDF5()
     int num_rings  = num_axial_blocks_per_bucket*num_axial_crystals_per_block*axial_modules_per_system;
     int num_detectors_per_ring = num_transaxial_blocks_per_bucket*num_transaxial_crystals_per_block*radial_modules_per_system;
     float ring_spacing = detector_axial_size/num_rings;
-
     //PW Bin Size, default num of arc corrected bins and inner ring radius not found in RDF header.
     // They will be set from the default STIR values
     shared_ptr<Scanner> scanner_sptr(Scanner::get_scanner_from_name(read_str_scanner));
@@ -389,7 +388,6 @@ shared_ptr<Scanner> GEHDF5Wrapper::get_scanner_from_HDF5()
       }
 
     return scanner_sptr;
-
 }
 
 void GEHDF5Wrapper::initialise_proj_data_info_from_HDF5()

src/IO/IO_registries.cxx

@@ -56,7 +56,6 @@
 #ifdef HAVE_HDF5
 #include "stir/IO/GEHDF5ListmodeInputFileFormat.h"
 #endif
-
 //! Addition for SAFIR listmode input file format
 #include "stir/IO/SAFIRCListmodeInputFileFormat.h"