Time Of Flight reconstruction

CMakeLists.txt

@@ -99,6 +99,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()

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/root_header.hroot

@@ -8,8 +8,11 @@ Average depth of interaction (cm)        := 0.7
 Distance between rings (cm)              := 0.40625
 Default bin size (cm)                    := 0.208626
 Maximum number of non-arc-corrected bins := 344
+Number of TOF time bins := 410
+Size of timing bin (in picoseconds) := 10.00
+Timing resolution (in picoseconds) := 400.0
 
-
+%TOF mashing factor:= 82
 
 GATE scanner type := GATE_Cylindrical_PET
 GATE_Cylindrical_PET Parameters :=

recon_test_pack/run_test_time_of_fligh.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
+

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 from 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

@@ -33,8 +33,6 @@ option(BUILD_EXECUTABLES
 option(BUILD_SHARED_LIBS 
        "Use shared libraries" OFF)
 
-### Settings for external libraries
-
 if (LLN_FOUND)
   set(HAVE_ECAT ON)
   message(STATUS "ECAT support enabled.")
@@ -68,6 +66,16 @@ else()
   message(STATUS "RDF support disabled.")
 endif()
 
+if (HDF5_FOUND)
+  set(HAVE_HDF5 ON)
+  message(STATUS "HDF5 support enabled.")
+  add_definitions(-D HAVE_HDF5)
+  include_directories(${HDF5_INCLUDE_DIRS})
+
+else()
+  message(STATUS "HDF5 support disabled.")
+endif()
+
 
 if (ITK_FOUND) 
   message(STATUS "ITK libraries added.")

src/IO/CMakeLists.txt

@@ -66,6 +66,12 @@ if (CERN_ROOT_FOUND)
  )
 endif()
 
+if (HDF5_FOUND)
+ list(APPEND ${dir_LIB_SOURCES}
+    GEHDF5Data
+ )
+endif()
+
 if(AVW_FOUND)
   list(APPEND ${dir_LIB_SOURCES}
     stir_AVW
@@ -82,6 +88,10 @@ if (CERN_ROOT_FOUND)
   target_link_libraries(IO ${CERN_ROOT_LIBRARIES})
 endif()
 
+if (HDF5_FOUND)
+  target_link_libraries(IO ${HDF5_CXX_LIBRARIES})
+endif()
+
 if (AVW_FOUND)
   target_link_libraries(IO ${AVW_LIBRARIES})
 endif()

src/IO/GEHDF5Data.cxx

@@ -0,0 +1,29 @@
+#include "stir/IO/GEHDF5Data.h"
+
+START_NAMESPACE_STIR
+
+shared_ptr<Scanner>
+GEHDF5Data::get_scanner_sptr() const
+{
+  return this->scanner_sptr;
+}
+
+shared_ptr<ExamInfo>
+GEHDF5Data::get_exam_info_sptr() const
+{
+  return this->exam_info_sptr;
+}
+
+void
+GEHDF5Data::open(const std::string& filename)
+{
+  this->file.openFile( filename, H5F_ACC_RDONLY );
+
+ warning("CListModeDataGESigna: "
+	  "Assuming this is GESigna, but couldn't find scan start time etc");
+  this->scanner_sptr.reset(new Scanner(Scanner::PETMR_Signa));   // for now we are assuming that the scanner is a PET/MR Signa
+// TODO: read scanner type from the dataset: "/HeaderData/ExamData/scannerDesc"
+  this->exam_info_sptr.reset(new ExamInfo); 
+}
+
+END_NAMESPACE_STIR

src/IO/IO_registries.cxx

@@ -59,10 +59,12 @@
 #endif
 #include "stir/IO/ECAT8_32bitListmodeInputFileFormat.h"
 
+#ifdef HAVE_HDF5
+#include "stir/IO/GESignaListmodeInputFileFormat.h"
+#endif
 //! Addition for SAFIR listmode input file format
 #include "stir/IO/SAFIRCListmodeInputFileFormat.h"
 
-
 //! Addition for ROOT support - Nikos Efthimiou
 #ifdef HAVE_CERN_ROOT
 #include "stir/IO/ROOTListmodeInputFileFormat.h"
@@ -139,5 +141,8 @@ static RegisterInputFileFormat<ecat::ecat7::ECAT966ListmodeInputFileFormat> LMdu
 static RegisterInputFileFormat<ecat::ecat7::ECAT962ListmodeInputFileFormat> LMdummyECAT962(5);
 #endif
 static RegisterInputFileFormat<ecat::ECAT8_32bitListmodeInputFileFormat> LMdummyECAT8(6);
+#ifdef HAVE_HDF5
+static RegisterInputFileFormat<GESignaListmodeInputFileFormat> LMdummyGESigna(7);
+#endif
 
 END_NAMESPACE_STIR

src/IO/InputStreamFromROOTFile.cxx

@@ -1,7 +1,7 @@
 /*
  *  Copyright (C) 2015, 2016 University of Leeds
     Copyright (C) 2016, UCL
-    Copyright (C) 2018 University of Hull
+    Copyright (C) 2016, 2017, 2018 University of Hull
     This file is part of STIR.
 
     This file is free software; you can redistribute it and/or modify
@@ -29,6 +29,7 @@ InputStreamFromROOTFile()
 {
     set_defaults();
     reset();
+    least_significant_clock_bit = 1.0e+12; // TODO remove cst or rename
 }