Generate dynamic images

src/Shape_buildblock/GenerateImage.cxx

@@ -3,6 +3,7 @@
 /*
     Copyright (C) 2003-2011, Hammersmith Imanet Ltd
     Copyright (C) 2018-2022, University College London
+    Copyright (C) 2023, Athinoula A. Martinos Center for Biomedical Imaging
     This file is part of STIR.
 
     SPDX-License-Identifier: Apache-2.0
@@ -70,7 +71,7 @@ set_defaults()
   rel_start_time = 0;
   output_filename.resize(0);
   output_file_format_sptr =
-          OutputFileFormat<DiscretisedDensity<3,float> >::default_sptr();
+          OutputFileFormat<DynamicDiscretisedDensity>::default_sptr();
 }
 
 void GenerateImage::set_imaging_modality()
@@ -122,6 +123,7 @@ initialise_keymap()
 
   add_key("image duration (sec)", &image_duration);
   add_key("image relative start time (sec)", &rel_start_time);
+  add_key("time frame definition filename", &frame_definition_filename);
 
   add_parsing_key("shape type", &current_shape_ptr);
   add_key("value", &current_value);
@@ -144,6 +146,32 @@ post_processing()
   exam_info_sptr->patient_position.set_rotation(static_cast<PatientPosition::RotationValue>(patient_rotation_index));
   exam_info_sptr->patient_position.set_orientation(static_cast<PatientPosition::OrientationValue>(patient_orientation_index));
 
+  if (frame_definition_filename.size()!=0)
+    {
+    TimeFrameDefinitions frame_defs(frame_definition_filename);
+    exam_info_sptr->set_time_frame_definitions(frame_defs);
+    }
+  else
+    {
+    if (image_duration>0.0)
+      {
+        std::vector<double> start_times(1, rel_start_time);
+        std::vector<double> durations(1, image_duration);
+        TimeFrameDefinitions frame_defs(start_times, durations);
+        exam_info_sptr->set_time_frame_definitions(frame_defs);
+      }
+    else
+      {
+        warning("image duration not set, so time frame definitions will not be initialised");
+        std::vector<double> start_times(1, rel_start_time);
+        std::vector<double> durations(1, 1);
+        TimeFrameDefinitions frame_defs(start_times, durations);
+        exam_info_sptr->set_time_frame_definitions(frame_defs);
+      }
+    //    std::vector<std::pair<double, double> > frame_times(1, std::pair<double, double>(0, 1));
+    //    frame_defs = TimeFrameDefinitions(frame_times);
+    }
+
   if (!is_null_ptr( current_shape_ptr))
   {
     shape_ptrs.push_back(current_shape_ptr);
@@ -204,6 +232,29 @@ post_processing()
     warning("number of samples to take in x-direction should be strictly positive\n");
     return true;
   }
+
+  tmpl_image.reset( new VoxelsOnCartesianGrid<float>(exam_info_sptr,
+                                                    IndexRange3D(0,output_image_size_z-1,
+                                                                 -(output_image_size_y/2),
+                                                                 -(output_image_size_y/2)+output_image_size_y-1,
+                                                                 -(output_image_size_x/2),
+                                                                 -(output_image_size_x/2)+output_image_size_x-1),
+                                                    CartesianCoordinate3D<float>(0,0,0),
+                                                    CartesianCoordinate3D<float>(output_voxel_size_z,
+                                                                                 output_voxel_size_y,
+                                                                                 output_voxel_size_x)));
+
+  shared_ptr<Scanner> scn(Scanner::get_scanner_from_name((exam_info_sptr->originating_system)));
+  out_density_ptr.reset(new DynamicDiscretisedDensity(exam_info_sptr->get_time_frame_definitions(),
+                                                      rel_start_time,
+                                                      scn,
+                                                      tmpl_image) );
+
+  for (unsigned int frame_num=1; frame_num <= exam_info_sptr->get_time_frame_definitions().get_num_frames(); ++frame_num)
+  {
+    out_density_ptr->get_density_sptr(frame_num)->fill(0.f);
+  }
+
   return false;
 }
 
@@ -252,42 +303,30 @@ compute()
 
 #else
 
-  if (image_duration>0.0)
-  {
-    std::vector<double> start_times(1, rel_start_time);
-    std::vector<double> durations(1, image_duration);
-    TimeFrameDefinitions frame_defs(start_times, durations);
-    exam_info_sptr->set_time_frame_definitions(frame_defs);
-  }
-  else
-  {
-    warning("image duration not set, so time frame definitions will not be initialised");
-  }
-  VoxelsOnCartesianGrid<float>
-          current_image(exam_info_sptr,
-                        IndexRange3D(0,output_image_size_z-1,
-                                     -(output_image_size_y/2),
-                                     -(output_image_size_y/2)+output_image_size_y-1,
-                                     -(output_image_size_x/2),
-                                     -(output_image_size_x/2)+output_image_size_x-1),
-                        CartesianCoordinate3D<float>(0,0,0),
-                        CartesianCoordinate3D<float>(output_voxel_size_z,
-                                                     output_voxel_size_y,
-                                                     output_voxel_size_x));
-  this->out_density_ptr.reset(current_image.clone());
+  for(int iframe = 1; iframe <= exam_info_sptr->get_time_frame_definitions().get_num_time_frames(); ++iframe)
+    {
+    shared_ptr<VoxelsOnCartesianGrid<float> > current_image = std::dynamic_pointer_cast<VoxelsOnCartesianGrid<float> >(out_density_ptr->get_density_sptr(iframe));
+//  this->out_density_ptr.reset(current_image.clone());
+    info(boost::format("Processing time frame %d ...")%iframe, 2);
 #endif
   std::vector<float >::const_iterator value_iter = values.begin();
   for (std::vector<shared_ptr<Shape3D> >::const_iterator iter = shape_ptrs.begin();
-          iter != shape_ptrs.end();
-  ++iter, ++value_iter)
-  {
-    info("Processing next shape...", 2);
-    current_image.fill(0);
-    (**iter).construct_volume(current_image, num_samples);
-    current_image *= *value_iter;
-    *out_density_ptr += current_image;
-  }
-  return Succeeded::yes;
+       iter != shape_ptrs.end();
+       ++iter, ++value_iter)
+    {
+    if( (**iter).is_in_frame(iframe))
+      {
+        info("Processing next shape...", 2);
+        VoxelsOnCartesianGrid<float> tmp_image = *tmpl_image->clone();
+
+        (**iter).construct_volume(tmp_image, num_samples);
+        tmp_image *= *value_iter;
+        *current_image += tmp_image;
+      }
+//    out_density_ptr->set_density(current_image, iframe+1);
+    }
+}
+return Succeeded::yes;
 }
 
 Succeeded
@@ -300,9 +339,15 @@ save_image()
 
 shared_ptr<DiscretisedDensity<3, float>>
 GenerateImage::
-get_output_sptr()
+get_output_sptr(unsigned int frame)
 {
-  return out_density_ptr;
+  return out_density_ptr->get_density_sptr(frame);
 }
 
+shared_ptr<DynamicDiscretisedDensity>
+    GenerateImage::
+    get_all_outputs_sptr()
+{
+  return out_density_ptr;
+}
 END_NAMESPACE_STIR

src/Shape_buildblock/Shape3D.cxx

@@ -16,6 +16,7 @@
 
   \author Kris Thielemans
   \author Sanida Mustafovic
+  \author Nikos Efthimiou
 */
 #include "stir/Shape/Shape3D.h"
 #include "stir/Shape/DiscretisedShape3D.h"
@@ -229,6 +230,15 @@ Shape3D::
 initialise_keymap()
 {
   this->parser.add_key("origin (in mm)", &origin);
+  this->parser.add_key("frames", &frames);
+}
+
+bool
+Shape3D::is_in_frame(const unsigned int this_frame) const
+{
+  if(frames.size() == 0)
+    return true;
+  return std::find(frames.begin(), frames.end(), this_frame-1)!=frames.end() ? true : false;
 }
 
 std::string 

src/buildblock/DynamicDiscretisedDensity.cxx

@@ -7,11 +7,13 @@
   \author Kris Thielemans
   \author Charalampos Tsoumpas
   \author Richard Brown
-
+  \author Nikos Efthimiou
+
 */
 /*
     Copyright (C) 2005- 2011, Hammersmith Imanet Ltd
     Copyright (C) 2018, University College London
+    Copyright (C) 2023, Athinoula A. Martinos Center for Biomedical Imaging
     This file is part of STIR.
 
     SPDX-License-Identifier: Apache-2.0
@@ -103,6 +105,12 @@ DynamicDiscretisedDensity::
 get_density(const unsigned int frame_num)
 {  return *this->_densities.at(frame_num-1) ; }
 
+shared_ptr<DiscretisedDensity<3, float> >
+DynamicDiscretisedDensity::get_density_sptr(const unsigned int frame_num) const
+{
+    return this->_densities.at(frame_num-1);
+}
+
 const float 
 DynamicDiscretisedDensity::
 get_isotope_halflife() const

src/include/stir/DynamicDiscretisedDensity.h

@@ -7,12 +7,14 @@
   \author Kris Thielemans
   \author Charalampos Tsoumpas
   \author Richard Brown
-
+  \author Nikos Efthimiou
+
 */
 /*
     Copyright (C) 2005 - 2011-01-12, Hammersmith Imanet Ltd
     Copyright (C) 2011-07-01 - 2011, Kris Thielemans
     Copyright (C) 2018-2020 University College London
+    Copyright (C) 2023, Athinoula A. Martinos Center for Biomedical Imaging
     This file is part of STIR.
 
     SPDX-License-Identifier: Apache-2.0
@@ -143,6 +145,9 @@ class DynamicDiscretisedDensity: public ExamData
   const singleDiscDensT & 
     get_density(const unsigned int frame_num) const ;
 
+  shared_ptr<DiscretisedDensity<3, float> >
+      get_density_sptr(const unsigned int frame_num) const;
+
   const singleDiscDensT & 
     operator[](const unsigned int frame_num) const 
     { return this->get_density(frame_num); }

src/include/stir/Shape/GenerateImage.h

@@ -3,6 +3,7 @@
 /*
     Copyright (C) 2003-2011, Hammersmith Imanet Ltd
     Copyright (C) 2018-2022, University College London
+    Copyright (C) 2023, Athinoula A. Martinos Center for Biomedical Imaging
     This file is part of STIR.
 
     SPDX-License-Identifier: Apache-2.0
@@ -18,6 +19,7 @@
   \author Kris Thielemans
   \author Sanida Mustafovic
   \author Robert Twyman
+  \author Nikos Efthimiou
 
   \par Example .par file
   \code
@@ -54,6 +56,9 @@
     scale_to_write_data:= 1
   End Interfile Output File Format Parameters:=
 
+  ; Used with simulation of dynamic images
+  ; time frame definition filename :=  frames.fdef
+
   X output image size (in pixels):= 13
   Y output image size (in pixels):= 13
   Z output image size (in pixels):= 15
@@ -77,6 +82,7 @@
      radius-y (in mm):= 2
      length-z (in mm):= 3
      origin (in mm):= {z,y,x}
+     ; frames := {1, 2 , 7}
      END:=
   value := 10
 
@@ -111,6 +117,7 @@
 #include "stir/Succeeded.h"
 #include "stir/IO/OutputFileFormat.h"
 #include "stir/VoxelsOnCartesianGrid.h"
+#include "stir/DynamicDiscretisedDensity.h"
 #include <iostream>
 
 
@@ -134,7 +141,9 @@ class GenerateImage : public KeyParser
     Succeeded save_image();
 
     //! Returns the discretised density with computed shapes.
-    shared_ptr<DiscretisedDensity<3, float>> get_output_sptr();
+    shared_ptr<DiscretisedDensity<3, float>> get_output_sptr(unsigned int frame = 1);
+
+    shared_ptr<DynamicDiscretisedDensity> get_all_outputs_sptr();
 
 private:
 
@@ -149,14 +158,16 @@ class GenerateImage : public KeyParser
     int patient_orientation_index;
     int patient_rotation_index;
 
-    shared_ptr<DiscretisedDensity<3, float> > out_density_ptr;
+    shared_ptr<DynamicDiscretisedDensity> out_density_ptr;
+    shared_ptr<VoxelsOnCartesianGrid<float> >
+        tmpl_image;
 
     std::vector<shared_ptr<Shape3D> > shape_ptrs;
     shared_ptr<Shape3D> current_shape_ptr;
     std::vector<float> values;
     float current_value;
     std::string output_filename;
-    shared_ptr<OutputFileFormat<DiscretisedDensity<3,float> > > output_file_format_sptr;
+    shared_ptr<OutputFileFormat<DynamicDiscretisedDensity> >output_file_format_sptr;
 
     void increment_current_shape_num();
 
@@ -171,6 +182,8 @@ class GenerateImage : public KeyParser
 
     double image_duration;
     double rel_start_time;
+    std::string frame_definition_filename;
+
 };
 
 END_NAMESPACE_STIR

src/include/stir/Shape/Shape3D.h

@@ -2,6 +2,7 @@
 //
 /*
     Copyright (C) 2000- 2007, Hammersmith Imanet Ltd
+    Copyright (C) 2023, Athinoula A. Martinos Center for Biomedical Imaging
     This file is part of STIR.
 
     SPDX-License-Identifier: Apache-2.0
@@ -16,6 +17,7 @@
 
   \author Kris Thielemans
   \author Sanida Mustafovic
+  \author Nikos Efthimiou
 */
 
 #ifndef __stir_Shape_Shape3D_h__
@@ -42,6 +44,9 @@ template <typename elemT> class VoxelsOnCartesianGrid;
   However, this then needs some special treatment for some member 
   functions, and you have to be somewhat careful with that class.
 
+  When dynamic images are generated the frames parameter can be used to
+  tell in which time frames this shape will be drawn.
+
   \todo This could/should be generalised to allow general fuzzy shapes.
   Probably the only thing to change is to let is_inside_shape() return
   a float (between 0 and 1). This would solve some issues with 
@@ -60,6 +65,7 @@ template <typename elemT> class VoxelsOnCartesianGrid;
   \verbatim
   ; specify origin as {z,y,x}
   origin (in mm):= <float> ;defaults to {0,0,0}
+  frames := {0,0,}
   \endverbatim
 */
 class Shape3D :
@@ -117,7 +123,9 @@ class Shape3D :
     \todo replace by floating point return value?
   */
   virtual bool is_inside_shape(const CartesianCoordinate3D<float>& coord) const = 0;
-
+
+  //! Draw if in the correct time dynamic time frame
+  bool is_in_frame(const unsigned int this_frame) const;
   //! translate the whole shape by shifting its origin 
   /*! Uses set_origin().
 
@@ -201,7 +209,8 @@ class Shape3D :
 private:
   //! origin of the shape
   CartesianCoordinate3D<float> origin;
-
+  //! time frames that we will add this shape
+  std::vector<int> frames;
 };
 
 END_NAMESPACE_STIR