minor clean-up and clarifications

Author: KrisThielemans

.gitignore

@@ -57,6 +57,7 @@ install_manifest.txt
 ~$*
 *~
 *.bak
+\#*#
 
 # Ignore Visual Studio User-specific files
 *.suo

src/include/stir/recon_buildblock/ProjMatrixByBin.inl

@@ -149,7 +149,7 @@ ProjMatrixByBin::get_tof_value(const float d1, const float d2) const
   if ((d1_n >= 4.f && d2_n >= 4.f) || (d1_n <= -4.f && d2_n <= -4.f))
     return 0.F;
   else
-    return 0.5f * (erf_interpolation(d2_n) - erf_interpolation(d1_n));
+    return static_cast<float>(0.5 * (erf_interpolation(d2_n) - erf_interpolation(d1_n)));
 }
 
 END_NAMESPACE_STIR

src/include/stir/recon_buildblock/ProjMatrixElemsForOneBin.h

@@ -178,14 +178,14 @@ class ProjMatrixElemsForOneBin
 
   //******************** projection operations ********************//
 
-  //! back project a single bin
+  //! back project a single bin (accumulates)
   void back_project(DiscretisedDensity<3, float>&, const Bin&) const;
 
-  //! forward project into a single bin
+  //! forward project into a single bin (accumulates)
   void forward_project(Bin&, const DiscretisedDensity<3, float>&) const;
-  //! back project related bins
+  //! back project related bins (accumulates)
   void back_project(DiscretisedDensity<3, float>&, const RelatedBins&) const;
-  //! forward project related bins
+  //! forward project related bins (accumulates)
   void forward_project(RelatedBins&, const DiscretisedDensity<3, float>&) const;
 
 private:

src/recon_buildblock/distributable.cxx

@@ -685,10 +685,10 @@ LM_distributable_computation_template(const shared_ptr<ProjMatrixByBin> PM_sptr,
 
         PM_sptr->get_proj_matrix_elems_for_one_bin(local_row[thread_num], measured_bin);
         call_back(*local_output_image_sptrs[thread_num],
-                  measured_bin,
-                  *input_image_ptr,
                   local_row[thread_num],
-                  has_add ? record.my_corr : 0);
+                  has_add ? record.my_corr : 0.F,
+                  measured_bin,
+                  *input_image_ptr);
       }
   }
 #ifdef STIR_OPENMP
@@ -710,56 +710,54 @@ LM_distributable_computation_template(const shared_ptr<ProjMatrixByBin> PM_sptr,
 
 constexpr float max_quotient = 10000.F;
 
+/* gradient without the sensitivity term
+
+\sum_e A_e^t (y_e/(A_e lambda+ c))
+*/
 inline void
 LM_gradient(DiscretisedDensity<3, float>& output_image,
-            const Bin& measured_bin,
-            const DiscretisedDensity<3, float>& input_image,
             const ProjMatrixElemsForOneBin& row,
-            const float add_term)
+            const float add_term,
+            const Bin& measured_bin,
+            const DiscretisedDensity<3, float>& input_image)
 {
   Bin fwd_bin = measured_bin;
   fwd_bin.set_bin_value(0.0f);
   row.forward_project(fwd_bin, input_image);
+  const auto fwd = fwd_bin.get_bin_value() + add_term;
 
-  if (add_term)
-    fwd_bin.set_bin_value(fwd_bin.get_bin_value() + add_term);
-
-  float measured_div_fwd = 0.F;
-  if (measured_bin.get_bin_value() <= max_quotient * fwd_bin.get_bin_value())
-    measured_div_fwd = measured_bin.get_bin_value() / fwd_bin.get_bin_value();
-  else
-    return;
+  if (measured_bin.get_bin_value() > max_quotient * fwd)
+    return; // cancel singularity
+  const auto measured_div_fwd = measured_bin.get_bin_value() / fwd;
 
   fwd_bin.set_bin_value(measured_div_fwd);
   row.back_project(output_image, fwd_bin);
 }
 
 /* Hessian
 
-\sum_e A_e^t (y_e/(A_e lambda+ c)^2 A_e x)
+\sum_e A_e^t (y_e/(A_e lambda+ c)^2 A_e rhs)
 */
 inline void
 LM_Hessian(DiscretisedDensity<3, float>& output_image,
+           const ProjMatrixElemsForOneBin& row,
+           const float add_term,
            const Bin& measured_bin,
            const DiscretisedDensity<3, float>& input_image,
-           const DiscretisedDensity<3, float>& rhs,
-           const ProjMatrixElemsForOneBin& row,
-           const float add_term)
+           const DiscretisedDensity<3, float>& rhs)
 {
   Bin fwd_bin = measured_bin;
   fwd_bin.set_bin_value(0.0f);
   row.forward_project(fwd_bin, input_image);
+  const auto fwd = fwd_bin.get_bin_value() + add_term;
 
-  if (add_term)
-    fwd_bin.set_bin_value(fwd_bin.get_bin_value() + add_term);
-
-  if (measured_bin.get_bin_value() > max_quotient * fwd_bin.get_bin_value())
-    return;
-  float measured_div_fwd2 = measured_bin.get_bin_value() / square(fwd_bin.get_bin_value());
+  if (measured_bin.get_bin_value() > max_quotient * fwd)
+    return; // cancel singularity
+  const auto measured_div_fwd2 = measured_bin.get_bin_value() / square(fwd);
 
+  // forward project rhs
   fwd_bin.set_bin_value(0.0f);
   row.forward_project(fwd_bin, rhs);
-
   if (fwd_bin.get_bin_value() == 0)
     return;