Fix memory leak from repeated ProjDataInfo::clone() calls (UCL#1673)

* Fix memory leak in ScatterSimulation from repeated ProjDataInfo::clone() calls

* Small edit to add a new method to ScatterSimulation class to set the template projection data info using a ProjDataInfo object directly, instead of only through a shared pointer. This new method creates a shared pointer from the given ProjDataInfo object and then calls the existing method that takes a shared pointer.

* minor clean-up of test_blocks_on_cylindrical.

Co-authored-by: Kris Thielemans <KrisThielemans@users.noreply.github.com>

Author: NikEfth

documentation/release_6.4.htm

@@ -73,6 +73,12 @@ <h3>Bug fixes</h3>
       segments had problems, see <a href=https://github.com/UCL/STIR/issues/1537># 1537</a>.<br>
       <a href=https://github.com/UCL/STIR/pull/1675># 1675</a>. Extra tests were introduced in <a href=https://github.com/UCL/STIR/pull/1674># 1674</a>.
     </li>
+      <li>
+    Fixed a small memory leak in scatter estimation caused by multiple calls to <code>ProjDataInfo::clone()</code> combined with <code>dynamic_cast</code>. 
+    The object is now cloned once and ownership is transferred safely after type checking.<br>
+    <a href=https://github.com/UCL/STIR/pull/1673>PR #1673</a>
+       </li>
+    <br>
   </ul>
 
   <h3>Build system</h3>

src/include/stir/scatter/ScatterSimulation.h

@@ -124,7 +124,9 @@ class ScatterSimulation : public RegisteredObject<ScatterSimulation>
 
   void set_template_proj_data_info(const std::string&);
 
-  void set_template_proj_data_info(const ProjDataInfo&);
+  void set_template_proj_data_info(const ProjDataInfo& arg);
+
+  void set_template_proj_data_info(shared_ptr<const ProjDataInfo> arg);
 
   void set_activity_image_sptr(const shared_ptr<const DiscretisedDensity<3, float>>);
 
@@ -334,7 +336,7 @@ class ScatterSimulation : public RegisteredObject<ScatterSimulation>
 
   std::string template_proj_data_filename;
 
-  shared_ptr<ProjDataInfo> proj_data_info_sptr;
+  shared_ptr<const ProjDataInfo> proj_data_info_sptr;
   //! \details Exam info extracted from the scanner template
   shared_ptr<ExamInfo> template_exam_info_sptr;
 

src/recon_test/test_blocks_on_cylindrical_projectors.cxx

@@ -195,7 +195,7 @@ BlocksTests::run_symmetry_test(ForwardProjectorByBin& forw_projector1, ForwardPr
   shared_ptr<DiscretisedDensity<3, float>> image1_sptr(image.clone());
   write_to_file("image_for", *image1_sptr);
 
-  image = *image.get_empty_copy();
+  image.fill(0);
   for (int i = 15; i < 360; i += 30)
     {
       theta2 = i * _PI / 180;
@@ -307,7 +307,8 @@ BlocksTests::run_plane_symmetry_test(ForwardProjectorByBin& forw_projector1, For
 
   //    rotate by 30 degrees
   phi2 = 30 * _PI / 180;
-  VoxelsOnCartesianGrid<float> image2 = *image.get_empty_copy();
+  VoxelsOnCartesianGrid<float> image2(image);
+  image2.fill(0);
   const Array<2, float> direction2 = make_array(make_1d_array(1.F, 0.F, 0.F),
                                                 make_1d_array(0.F, cos(float(_PI) - phi2), sin(float(_PI) - phi2)),
                                                 make_1d_array(0.F, -sin(float(_PI) - phi2), cos(float(_PI) - phi2)));
@@ -591,7 +592,7 @@ BlocksTests::run_map_orientation_test(ForwardProjectorByBin& forw_projector1, Fo
   shared_ptr<DiscretisedDensity<3, float>> image1_sptr(image.clone());
   write_to_file("image_to_fwp", *image1_sptr);
 
-  image = *image.get_empty_copy();
+  image.fill(0);
 
   shared_ptr<const DetectorCoordinateMap> map_sptr;
   auto scannerBlocks_sptr = std::make_shared<Scanner>(Scanner::SAFIRDualRingPrototype);
@@ -700,7 +701,7 @@ BlocksTests::run_projection_test(ForwardProjectorByBin& forw_projector1, Forward
   shared_ptr<DiscretisedDensity<3, float>> image1_sptr(image.clone());
   write_to_file("image_with_voxel_at_30_0", *image1_sptr);
 
-  image = *image.get_empty_copy();
+  image.fill(0);
   image[(image.get_min_index() + image.get_max_index()) / 2 * grid_spacing.z()][0][-25] = 1;
 
   shared_ptr<DiscretisedDensity<3, float>> image2_sptr(image.clone());

src/scatter_buildblock/ScatterEstimation.cxx

@@ -610,12 +610,12 @@ ScatterEstimation::set_up()
            "discarded)");
       if (run_in_2d_projdata)
         {
-          this->scatter_simulation_sptr->set_template_proj_data_info(*this->input_projdata_2d_sptr->get_proj_data_info_sptr());
+          this->scatter_simulation_sptr->set_template_proj_data_info(this->input_projdata_2d_sptr->get_proj_data_info_sptr());
           this->scatter_simulation_sptr->set_exam_info(this->input_projdata_2d_sptr->get_exam_info());
         }
       else
         {
-          this->scatter_simulation_sptr->set_template_proj_data_info(*this->input_projdata_sptr->get_proj_data_info_sptr());
+          this->scatter_simulation_sptr->set_template_proj_data_info(this->input_projdata_sptr->get_proj_data_info_sptr());
           this->scatter_simulation_sptr->set_exam_info(this->input_projdata_sptr->get_exam_info());
         }
     }

src/scatter_buildblock/ScatterSimulation.cxx

@@ -352,14 +352,14 @@ ScatterSimulation::set_up()
   {
     CartesianCoordinate3D<float> detector_coord_A, detector_coord_B;
     // check above statement
-    if (dynamic_cast<ProjDataInfoCylindricalNoArcCorr*>(proj_data_info_sptr.get()))
+    if (dynamic_cast<const ProjDataInfoCylindricalNoArcCorr*>(proj_data_info_sptr.get()))
       {
-        auto ptr = dynamic_cast<ProjDataInfoCylindricalNoArcCorr*>(proj_data_info_sptr.get());
+        const auto* ptr = dynamic_cast<const ProjDataInfoCylindricalNoArcCorr*>(proj_data_info_sptr.get());
         ptr->find_cartesian_coordinates_of_detection(detector_coord_A, detector_coord_B, Bin(0, 0, 0, 0));
       }
     else
       {
-        auto ptr = dynamic_cast<ProjDataInfoBlocksOnCylindricalNoArcCorr*>(proj_data_info_sptr.get());
+        const auto* ptr = dynamic_cast<const ProjDataInfoBlocksOnCylindricalNoArcCorr*>(proj_data_info_sptr.get());
         ptr->find_cartesian_coordinates_of_detection(detector_coord_A, detector_coord_B, Bin(0, 0, 0, 0));
       }
 
@@ -374,14 +374,14 @@ ScatterSimulation::set_up()
     assert(fabs(m_last + m_first) < m_last * 10E-4);
   }
 #endif
-  if (dynamic_cast<ProjDataInfoCylindricalNoArcCorr*>(proj_data_info_sptr.get()))
+  if (dynamic_cast<const ProjDataInfoCylindricalNoArcCorr*>(proj_data_info_sptr.get()))
     {
       this->shift_detector_coordinates_to_origin
           = CartesianCoordinate3D<float>(this->proj_data_info_sptr->get_m(Bin(0, 0, 0, 0)), 0, 0);
     }
   else
     {
-      if (dynamic_cast<ProjDataInfoBlocksOnCylindricalNoArcCorr*>(proj_data_info_sptr.get()))
+      if (dynamic_cast<const ProjDataInfoBlocksOnCylindricalNoArcCorr*>(proj_data_info_sptr.get()))
         {
           // align BlocksOnCylindrical scanner ring 0 to z=0.
           this->shift_detector_coordinates_to_origin
@@ -736,22 +736,31 @@ ScatterSimulation::set_template_proj_data_info(const std::string& filename)
 
   this->set_exam_info(template_proj_data_sptr->get_exam_info());
 
-  this->set_template_proj_data_info(*template_proj_data_sptr->get_proj_data_info_sptr());
+  this->set_template_proj_data_info(template_proj_data_sptr->get_proj_data_info_sptr());
 }
 
 void
 ScatterSimulation::set_template_proj_data_info(const ProjDataInfo& arg)
 {
-  this->_already_set_up = false;
-  this->proj_data_info_sptr.reset(dynamic_cast<ProjDataInfoBlocksOnCylindricalNoArcCorr*>(arg.clone()));
+  shared_ptr<const ProjDataInfo> sptr(arg.create_shared_clone());
+  this->set_template_proj_data_info(sptr);
+}
 
-  if (is_null_ptr(this->proj_data_info_sptr))
+void
+ScatterSimulation::set_template_proj_data_info(shared_ptr<const ProjDataInfo> arg)
+{
+  this->_already_set_up = false;
+  if (auto p = std::dynamic_pointer_cast<const ProjDataInfoBlocksOnCylindricalNoArcCorr>(arg))
     {
-      this->proj_data_info_sptr.reset(dynamic_cast<ProjDataInfoCylindricalNoArcCorr*>(arg.clone()));
-      if (is_null_ptr(this->proj_data_info_sptr))
-        {
-          error("ScatterSimulation: Can only handle non-arccorrected data");
-        }
+      this->proj_data_info_sptr = p;
+    }
+  else if (auto q = std::dynamic_pointer_cast<const ProjDataInfoCylindricalNoArcCorr>(arg))
+    {
+      this->proj_data_info_sptr = q;
+    }
+  else
+    {
+      error("ScatterSimulation: Can only handle non-arccorrected data");
     }
 
   // find final size of detection_points_vector
@@ -944,7 +953,7 @@ ScatterSimulation::downsample_scanner(int new_num_rings, int new_num_dets)
                                     false));
 
   info(format("ScatterSimulation: down-sampled scanner info:\n{}", templ_proj_data_info_sptr->parameter_info()), 3);
-  this->set_template_proj_data_info(*templ_proj_data_info_sptr);
+  this->set_template_proj_data_info(templ_proj_data_info_sptr);
   this->set_output_proj_data(this->output_proj_data_filename);
 
   return Succeeded::yes;

src/scatter_buildblock/scatter_detection_modelling.cxx

@@ -65,14 +65,14 @@ ScatterSimulation::find_detectors(unsigned& det_num_A, unsigned& det_num_B, cons
     error("ScatterSimulation::find_detectors: need to call set_up() first");
 #endif
   CartesianCoordinate3D<float> detector_coord_A, detector_coord_B;
-  auto ptr = dynamic_cast<ProjDataInfoBlocksOnCylindricalNoArcCorr*>(proj_data_info_sptr.get());
+  const auto* ptr = dynamic_cast<const ProjDataInfoBlocksOnCylindricalNoArcCorr*>(proj_data_info_sptr.get());
   if (ptr)
     {
       ptr->find_cartesian_coordinates_of_detection(detector_coord_A, detector_coord_B, bin);
     }
   else
     {
-      auto ptr = dynamic_cast<ProjDataInfoCylindricalNoArcCorr*>(proj_data_info_sptr.get());
+      const auto* ptr = dynamic_cast<const ProjDataInfoCylindricalNoArcCorr*>(proj_data_info_sptr.get());
       if (ptr)
         {
           ptr->find_cartesian_coordinates_of_detection(detector_coord_A, detector_coord_B, bin);

src/test/test_ScatterSimulation.cxx

@@ -92,7 +92,7 @@ ScatterSimulationTests::test_downsampling_ProjDataInfo()
                                     false)));
 
   unique_ptr<SingleScatterSimulation> sss(new SingleScatterSimulation());
-  sss->set_template_proj_data_info(*original_projdata);
+  sss->set_template_proj_data_info(original_projdata);
   {
     auto sss_projdata = dynamic_cast<const ProjDataInfoCylindricalNoArcCorr*>(sss->get_template_proj_data_info_sptr().get());
     check(*original_projdata == *sss_projdata, "Check the ProjDataInfo has been set correctly.");
@@ -183,7 +183,7 @@ ScatterSimulationTests::test_downsampling_DiscretisedDensity()
 
   unique_ptr<SingleScatterSimulation> sss(new SingleScatterSimulation());
 
-  sss->set_template_proj_data_info(*original_projdata);
+  sss->set_template_proj_data_info(original_projdata);
   sss->set_density_image_sptr(atten_density);
 
   //    int total_scatter_points_orig = sss.get_num_scatter_points();
@@ -391,7 +391,7 @@ ScatterSimulationTests::test_scatter_simulation()
   //// sss settings
   sss->set_randomly_place_scatter_points(false);
 
-  sss->set_template_proj_data_info(*original_projdata_info);
+  sss->set_template_proj_data_info(original_projdata_info);
   sss->downsample_scanner(original_projdata_info->get_scanner_sptr()->get_num_rings(), -1);
 #if 1
   set_tolerance(.02);
@@ -448,7 +448,7 @@ ScatterSimulationTests::test_scatter_simulation()
   }
 
   // a few tests with more downsampled scanner
-  sss->set_template_proj_data_info(*original_projdata_info);
+  sss->set_template_proj_data_info(original_projdata_info);
   sss->downsample_scanner(original_projdata_info->get_scanner_sptr()->get_num_rings() / 2, -1);
   {
     const int new_size_z = 14;

src/test/test_interpolate_projdata.cxx

@@ -850,7 +850,7 @@ InterpolationTests::transaxial_upsampling_interpolation_test_blocks()
 
   // use the code in scatter simulation to downsample the scanner
   auto scatter_simulation = SingleScatterSimulation();
-  scatter_simulation.set_template_proj_data_info(*proj_data_info);
+  scatter_simulation.set_template_proj_data_info(proj_data_info);
   scatter_simulation.set_exam_info(exam_info);
   scatter_simulation.downsample_scanner(-1, 96 / 4); // number of detectors per ring reduced by factor of four
   auto downsampled_proj_data_info = scatter_simulation.get_template_proj_data_info_sptr();