2nd - nested Patlak and generalized Patlak OSMAPOSL (POSMAPOSL & GPOSMAPOSL) and nested Richardson-Lucy Motion Deconvolution Image Reconstruction (RLMCIR)

recon_test_pack/run_tests_modelling.sh

@@ -126,7 +126,7 @@ generate_image ${INPUTDIR}generate_p5.par
 fi
 assemble_images p0005-p5.${imgext} p0005.hv p5.hv
 # test
-extract_single_images_from_parametric_image p0005-p5_%d.hv p0005-p5.${imgext}
+extract_single_images_from_parametric_image p0005-p5_{}.hv p0005-p5.${imgext}
 if compare_image p0005-p5_1.hv p0005.hv; then
     : # ok
 else
@@ -145,14 +145,14 @@ get_dynamic_images_from_parametric_images dyn_from_p0005-p5.hv p0005-p5.${imgext
 # run and test Patlak estimation
 apply_patlak_to_images indirect_Patlak.hv dyn_from_p0005-p5.hv ${INPUTDIR}PatlakPlot.par
 echo "Test Patlak round-trip"
-extract_single_images_from_parametric_image indirect_Patlak_img_%d.hv indirect_Patlak.hv
+extract_single_images_from_parametric_image indirect_Patlak_img_{}.hv indirect_Patlak.hv
 echo "indirect to original"
 for par in 1 2; do
     compare_image indirect_Patlak_img_${par}.hv p0005-p5_${par}.hv
 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
+extract_single_images_from_dynamic_image dyn_from_p0005-p5_img_f{}g1d0b0.hv dyn_from_p0005-p5.hv
 # if [ ! -r fwd_dyn_from_p0005-p5.S ]; then
 
 for fr in `count 23 28`; do
@@ -184,14 +184,49 @@ for direct in OSMAPOSL OSSPS ; do
     INPUT=fwd_dyn_from_p0005-p5.hs INIT=indirect_Patlak.hv ${MPIRUN} P${direct} P${direct}.par > P${direct}.log 2>&1
 
     echo "Compare the direct parametric images to the original ones"
-    extract_single_images_from_parametric_image p0005-p5_img_f%dg1d0b0.hv p0005-p5.${imgext}
-    extract_single_images_from_parametric_image P${direct}_${ITER}_img_f%dg1d0b0.hv P${direct}_${ITER}.hv
+    extract_single_images_from_parametric_image p0005-p5_img_f{}g1d0b0.hv p0005-p5.${imgext}
+    extract_single_images_from_parametric_image P${direct}_${ITER}_img_f{}g1d0b0.hv P${direct}_${ITER}.hv
     for par in 1 2; do
         compare_image -t .01 P${direct}_${ITER}_img_f${par}g1d0b0.hv p0005-p5_img_f${par}g1d0b0.hv  
     done
     echo "Comparison is OK"
 
 done # POSMAPOSL POSSPS#
+
+for nested in NESTPOSMAPOSL ; do #NESTGPOSMAPOSL ; do 
+    cp ${INPUTDIR}${nested}.par . 
+    echo "Test the nested ${nested} Patlak Plot reconstruction" 
+    # rm -f ${nested}.log
+    INPUT=fwd_dyn_from_p0005-p5.hs INIT=indirect_Patlak.hv ${MPIRUN} ${nested} ${nested}.par > ${nested}.log 2>&1
+
+    case ${nested} in
+        NESTPOSMAPOSL)
+            par_map="1:1 2:2"
+            ;;
+        NESTGPOSMAPOSL)
+            par_map="1:1 3:2"
+            ;;
+        *)
+            echo "Unknown algorithm ${nested}" >&2
+            exit 1
+            ;;
+    esac
+
+    echo "Compare the nested direct parametric images to the original ones"
+    extract_single_images_from_parametric_image p0005-p5_img_f{}g1d0b0.hv p0005-p5.${imgext}
+    extract_single_images_from_parametric_image ${nested}_${ITER}_img_f{}g1d0b0.hv ${nested}_${ITER}.hv
+
+    for pair in ${par_map}; do
+        recon_par=${pair%%:*}
+        truth_par=${pair##*:}
+        compare_image -t .01 \
+            ${nested}_${ITER}_img_f${recon_par}g1d0b0.hv \
+            p0005-p5_img_f${truth_par}g1d0b0.hv
+    done
+
+    echo "Comparison is OK"
+
+done # NESTPOSMAPOSL NESTGPOSMAPOSL
 
 echo "Test the utility: 'mult_model_with_dyn_images'"
 echo "Multiply the  dynamic images with the model matrix to get images in the parametric space."

recon_test_pack/test_modelling_input/NESTGPOSMAPOSL.par

@@ -0,0 +1,100 @@
+OSMAPOSLParameters :=
+
+objective function type:=PoissonNestedLogLikelihoodWithGeneralizedPatlakAndDynamicProjectionData
+PoissonNestedLogLikelihoodWithGeneralizedPatlakAndDynamicProjectionData Parameters :=
+
+input file := ${INPUT}
+
+; if disabled, defaults to maximum segment number in the file
+maximum absolute segment number to process := ${MAXSEG}
+; see User's Guide to see when you need this
+zero end planes of segment 0:= 0
+
+;zoom := .15645
+;xy output image size (in pixels) := 21
+;z output image size (in pixels) := 17
+;z offset (in mm) := 0
+
+projector pair type := Matrix
+  Projector Pair Using Matrix Parameters :=
+  Matrix type := Ray Tracing
+  Ray tracing matrix parameters :=
+  End Ray tracing matrix parameters :=
+  End Projector Pair Using Matrix Parameters :=
+
+Bin Normalisation type := From ProjData
+  Bin Normalisation From ProjData :=
+    normalisation projdata filename := all_ones.hs
+  End Bin Normalisation From ProjData :=
+
+; we need this for backwards compatibility with the testing script
+use subset sensitivities:=0
+sensitivity filename:= sens.img
+; if next is set to 1, sensitivity will be recomputed
+; and written to file (if "sensitivity filename" is set)
+recompute sensitivity := 1
+
+; specify additive projection data to handle randoms or so
+; see User's Guide for more info
+additive sinograms := 0
+
+; patlak related files
+Kinetic Model type := Generalized Patlak Plot
+Generalized Patlak Plot Parameters :=
+time frame definition filename := time.fdef
+starting frame := 23
+calibration factor := 9432.31
+blood data filename := plasma.if 
+Time Shift := 0
+In total counts := 1
+;In correct scale := 0
+kloss lower bound := 0
+kloss upper bound := 0
+; Default: end_of_frame = 0
+;mid_of_frame = 1
+frame reference time: = 0
+number of kloss samples := 3
+end Generalized Patlak Plot Parameters :=
+
+;; Nested subiterations (regular nested iterations using the Generalized Patlak matrix)
+number of nested subiterations := 0
+;; run the nested linear (standard) Patlak EM to convergence first, 
+;;  and use that as the starting point for the nested generalized Patlak EM,
+;; Global initial subiterations for initialization (using standard Patlak model, 
+;;   i.e. ModelMatrix) before switching to Generalized Patlak iterations
+number of global initialization subiterations := 40
+;; Nested subiterations for initialization (using standard Patlak model, 
+;;   i.e. ModelMatrix instead of GenearlizedPatlakMatrix)
+number of nested initialization subiterations := 1
+;; Choose whether both standard and generalized Patlak models 
+;; (1) will be used in an alternating fashion for the initialization or 
+;; (0) only the standard Patlak model will be utilized.
+alternating initialization mode := 0 
+; maximum nested relative change := 
+; minimum nested relative change :=
+
+End PoissonNestedLogLikelihoodWithGeneralizedPatlakAndDynamicProjectionData Parameters:=
+
+; Number of subsets should be a divisor of num_views/4
+number of subsets:=${NUMSUBS}
+; Use for starting the numbering from something else than 1
+start at subiteration number:=1
+; Use if you want to start from another subset than 0 (but why?)
+start at subset:= 0
+number of subiterations:= ${ITER}
+save estimates at subiteration intervals:= ${SAVITER}
+
+initial estimate := ${INIT}
+; enable this when you read an initial estimate with negative data
+enforce initial positivity condition:=1
+
+inter-update filter subiteration interval:= 0
+inter-update filter type := None
+
+inter-iteration filter subiteration interval:= 0
+
+post-filter type := None
+
+output filename prefix:= NESTGPOSMAPOSL
+
+END :=

recon_test_pack/test_modelling_input/NESTPOSMAPOSL.par

@@ -0,0 +1,87 @@
+OSMAPOSLParameters :=
+
+objective function type:=PoissonNestedLogLikelihoodWithLinearKineticModelAndDynamicProjectionData
+PoissonNestedLogLikelihoodWithLinearKineticModelAndDynamicProjectionData Parameters :=
+
+input file := ${INPUT}
+
+; if disabled, defaults to maximum segment number in the file
+maximum absolute segment number to process := ${MAXSEG}
+; see User's Guide to see when you need this
+zero end planes of segment 0:= 0
+
+;zoom := .15645
+;xy output image size (in pixels) := 21
+;z output image size (in pixels) := 17
+;z offset (in mm) := 0
+
+projector pair type := Matrix
+  Projector Pair Using Matrix Parameters :=
+  Matrix type := Ray Tracing
+  Ray tracing matrix parameters :=
+  End Ray tracing matrix parameters :=
+  End Projector Pair Using Matrix Parameters :=
+
+Bin Normalisation type := From ProjData
+  Bin Normalisation From ProjData :=
+    normalisation projdata filename := all_ones.hs
+  End Bin Normalisation From ProjData :=
+
+; we need this for backwards compatibility with the testing script
+use subset sensitivities:=0
+sensitivity filename:= sens.img
+; if next is set to 1, sensitivity will be recomputed
+; and written to file (if "sensitivity filename" is set)
+recompute sensitivity := 1
+
+; specify additive projection data to handle randoms or so
+; see User's Guide for more info
+additive sinograms := 0
+
+; patlak related files
+Kinetic Model type := Patlak Plot
+Patlak Plot Parameters :=
+time frame definition filename := time.fdef
+starting frame := 23
+calibration factor := 9432.31
+blood data filename := plasma.if 
+Time Shift := 0
+In total counts := 1
+;In correct scale := 0
+end Patlak Plot Parameters :=
+
+;; These are unique parameters in Nested methods
+;; 20 subsets as per: Whole-body direct 4D parametric PET imaging employing nested 
+;; generalized Patlak expectation–maximization reconstruction, Phys Med Biol 2016
+number of nested subiterations := 1
+; restrict updates (larger nested relative updates will be thresholded)
+; maximum nested relative change := 
+; restrict updates (smaller nested relative updates will be thresholded)
+; minimum nested relative change := 
+
+
+End PoissonNestedLogLikelihoodWithLinearKineticModelAndDynamicProjectionData Parameters:=
+
+; Number of subsets should be a divisor of num_views/4
+number of subsets:=${NUMSUBS}
+; Use for starting the numbering from something else than 1
+start at subiteration number:=1
+; Use if you want to start from another subset than 0 (but why?)
+start at subset:= 0
+number of subiterations:= ${ITER}
+save estimates at subiteration intervals:= ${SAVITER}
+
+initial estimate := ${INIT}
+; enable this when you read an initial estimate with negative data
+enforce initial positivity condition:=1
+
+inter-update filter subiteration interval:= 0
+inter-update filter type := None
+
+inter-iteration filter subiteration interval:= 0
+
+post-filter type := None
+
+output filename prefix:=NESTPOSMAPOSL
+
+END :=

src/IO/CMakeLists.txt

@@ -23,6 +23,7 @@ set(${dir_LIB_SOURCES}
   InterfileHeader.cxx
   InterfilePDFSHeaderSPECT.cxx
   InputFileFormatRegistry.cxx
+  InterfileDynamicDiscretisedDensityOutputFileFormat.cxx 
 ) 
 
 if (NOT MINI_STIR)

src/IO/ECAT7ParametricDensityOutputFileFormat.cxx

@@ -14,6 +14,7 @@
   \brief Implementation of class stir::ecat::ecat7::ECAT7ParametricDensityOutputFileFormat
 
   \author Kris Thielemans
+  \author Nicolas A Karakatsanis
 
 */
 #include "stir/IO/ECAT7ParametricDensityOutputFileFormat.h"
@@ -155,7 +156,8 @@ ECAT7ParametricDensityOutputFileFormat<DiscretisedDensityT>::actual_write_to_fil
   return Succeeded::yes;
 }
 
-template class ECAT7ParametricDensityOutputFileFormat<ParametricVoxelsOnCartesianGridBaseType>;
+template class ECAT7ParametricDensityOutputFileFormat<Parametric2VoxelsOnCartesianGridBaseType>;
+template class ECAT7ParametricDensityOutputFileFormat<Parametric3VoxelsOnCartesianGridBaseType>;
 
 END_NAMESPACE_ECAT7
 END_NAMESPACE_ECAT

src/IO/IO_registries.cxx

@@ -20,6 +20,8 @@
 
   \author Kris Thielemans
   \author Berta Marti Fuster
+  \author Nicolas A Karakatsanis
+  \author Nikos Efthimiou
 */
 
 #include "stir/IO/InterfileOutputFileFormat.h"
@@ -93,9 +95,12 @@ static RegisterInputFileFormat<InterfileImageInputFileFormat> idummy0(0);
 static ITKOutputFileFormat::RegisterIt dummyITK1;
 #  endif
 static InterfileDynamicDiscretisedDensityOutputFileFormat::RegisterIt dummydynIntfOut;
-static InterfileParametricDiscretisedDensityOutputFileFormat<ParametricVoxelsOnCartesianGridBaseType>::RegisterIt dummyparIntfOut;
+static InterfileParametricDiscretisedDensityOutputFileFormat<Parametric2VoxelsOnCartesianGridBaseType>::RegisterIt
+    dummyparIntfOut;
+static InterfileParametricDiscretisedDensityOutputFileFormat<Parametric3VoxelsOnCartesianGridBaseType>::RegisterIt
+    dummyGenPatIntfIn;
 static MultiDynamicDiscretisedDensityOutputFileFormat::RegisterIt dummydynMultiOut;
-static MultiParametricDiscretisedDensityOutputFileFormat<ParametricVoxelsOnCartesianGridBaseType>::RegisterIt dummyparMultiOut;
+static MultiParametricDiscretisedDensityOutputFileFormat<Parametric2VoxelsOnCartesianGridBaseType>::RegisterIt dummyparMultiOut;
 
 //! Support for SAFIR listmode file format
 static RegisterInputFileFormat<SAFIRCListmodeInputFileFormat<CListEventDataSAFIR>> LMdummySAFIR(4);
@@ -119,7 +124,8 @@ END_NAMESPACE_ECAT6
 START_NAMESPACE_ECAT7
 static ECAT7OutputFileFormat::RegisterIt dummy3;
 static ECAT7DynamicDiscretisedDensityOutputFileFormat::RegisterIt dummydynecat7In;
-static ECAT7ParametricDensityOutputFileFormat<ParametricVoxelsOnCartesianGridBaseType>::RegisterIt dummyparecat7In;
+static ECAT7ParametricDensityOutputFileFormat<Parametric2VoxelsOnCartesianGridBaseType>::RegisterIt dummyparecat7In;
+static ECAT7ParametricDensityOutputFileFormat<Parametric3VoxelsOnCartesianGridBaseType>::RegisterIt dummyGenPatecat7In;
 END_NAMESPACE_ECAT7
 END_NAMESPACE_ECAT
 #  endif
@@ -139,7 +145,8 @@ static RegisterInputFileFormat<ITKImageInputFileFormat<DiscretisedDensity<3, flo
 static RegisterInputFileFormat<ITKImageInputFileFormat<DiscretisedDensity<3, CartesianCoordinate3D<float>>>> idummy7(10000);
 #  endif
 static RegisterInputFileFormat<InterfileDynamicDiscretisedDensityInputFileFormat> dyndummy_intf(1);
-static RegisterInputFileFormat<InterfileParametricDiscretisedDensityInputFileFormat> paradummy_intf(1);
+static RegisterInputFileFormat<InterfileParametricDiscretisedDensityInputFileFormat<2>> paradummy_intf(1);
+static RegisterInputFileFormat<InterfileParametricDiscretisedDensityInputFileFormat<3>> paradummy_int3f(1);
 static RegisterInputFileFormat<MultiDynamicDiscretisedDensityInputFileFormat> dynim_dummy_multi(1);
 static RegisterInputFileFormat<MultiParametricDiscretisedDensityInputFileFormat> parim_dummy_multi(1);
 

src/IO/InputFileFormatRegistry.cxx

@@ -13,6 +13,7 @@
   \brief Instantiations for class stir::InputFileFormatRegistry
 
   \author Kris Thielemans
+  \author Nicolas A Karakatsanis
 
 */
 
@@ -28,8 +29,9 @@ START_NAMESPACE_STIR
 // instantiations
 template class InputFileFormatRegistry<DiscretisedDensity<3, float>>;
 template class InputFileFormatRegistry<ParametricVoxelsOnCartesianGrid>;
+template class InputFileFormatRegistry<Parametric3VoxelsOnCartesianGrid>;
 template class InputFileFormatRegistry<DynamicDiscretisedDensity>;
 template class InputFileFormatRegistry<ListModeData>;
 template class InputFileFormatRegistry<DiscretisedDensity<3, CartesianCoordinate3D<float>>>;
 
-END_NAMESPACE_STIR
+END_NAMESPACE_STIR

src/IO/InterfileHeader.cxx

@@ -389,6 +389,9 @@ InterfileHeader::post_processing()
 
   exam_info_sptr->time_frame_definitions = TimeFrameDefinitions(image_relative_start_times, image_durations);
 
+  // Added for old implementations relying on this->time_frame_definitions variable
+  this->time_frame_definitions = exam_info_sptr->time_frame_definitions;
+
   return false;
 }
 

src/IO/InterfileParametricDiscretisedDensityOutputFileFormat.cxx

@@ -103,6 +103,7 @@ InterfileParamDiscDensity::actual_write_to_file(std::string& filename,
 #undef ParamDiscDensity
 #undef TEMPLATE
 
-template class InterfileParametricDiscretisedDensityOutputFileFormat<ParametricVoxelsOnCartesianGridBaseType>;
+template class InterfileParametricDiscretisedDensityOutputFileFormat<Parametric2VoxelsOnCartesianGridBaseType>;
+template class InterfileParametricDiscretisedDensityOutputFileFormat<Parametric3VoxelsOnCartesianGridBaseType>;
 
 END_NAMESPACE_STIR

src/IO/MultiParametricDiscretisedDensityOutputFileFormat.cxx

@@ -126,6 +126,6 @@ ParamDiscDensityOutputFileFormat::actual_write_to_file(std::string& filename,
 #undef ParamDiscDensity
 #undef TEMPLATE
 
-template class MultiParametricDiscretisedDensityOutputFileFormat<ParametricVoxelsOnCartesianGridBaseType>;
+template class MultiParametricDiscretisedDensityOutputFileFormat<Parametric2VoxelsOnCartesianGridBaseType>;
 
 END_NAMESPACE_STIR