change Poisson log-likelihood truncation/thresholding strategy

src/include/stir/ArrayFunction.h

@@ -1,6 +1,7 @@
 /*
     Copyright (C) 2000 PARAPET partners
     Copyright (C) 2000- 2007, Hammersmith Imanet Ltd
+    Copyright (C) 2024, University College London
     This file is part of STIR.
 
     SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license
@@ -22,37 +23,30 @@
    functions which work on all stir::Array objects, and which change every element of the
    array:
    <ul>
-     <li> stir::in_place_log, stir::in_place_exp (these work only well when elements are float or double)
-     <li>stir::in_place_abs (does not work on complex numbers)
-     <li>stir::in_place_apply_function
-     <li>stir::in_place_apply_array_function_on_1st_index
-     <li>stir::in_place_apply_array_functions_on_each_index
+     <li>stir::in_place_log, stir::in_place_exp (these work only well when elements are float or double)</li>
+     <li>stir::in_place_abs (does not work on complex numbers)</li>
+     <li>stir::in_place_apply_function</li>
+     <li>stir::in_place_apply_array_function_on_1st_index</li>
+     <li>stir::in_place_apply_array_functions_on_each_index</li>
    </ul>
    All these functions return a reference to the (modified) array
+   </li>
    <li>
    Analoguous functions that take out_array and in_array
    <ul>
-      <li>stir::apply_array_function_on_1st_index
-      <li>stir::apply_array_functions_on_each_index
+      <li>stir::apply_array_function_on_1st_index</li>
+      <li>stir::apply_array_functions_on_each_index</li>
    </ul>
+   </li>
+   <li>
+   Functions that apply a binary function element-wise on arrays:
+   <ul>
+   <li>std::transform</li>
+   </ul>
+   </li>
    </ul>
-
-   \warning Compilers without partial specialisation of templates are
-   catered for by explicit instantiations. If you need it for any other
-   types, you'd have to add them by hand.
  */
 
-/* History:
-
-  KT 21/05/2001
-  added in_place_apply_array_function_on_1st_index,
-  in_place_apply_array_function_on_each_index
-
-  KT 06/12/2001
-  added apply_array_function_on_1st_index,
-  apply_array_function_on_each_index
-*/
-
 #ifndef __stir_ArrayFunction_H__
 #define __stir_ArrayFunction_H__
 
@@ -253,15 +247,77 @@ inline void apply_array_functions_on_each_index(Array<1, elemT>& out_array,
                                                 ActualFunctionObjectPtrIter start,
                                                 ActualFunctionObjectPtrIter stop);
 
-template <typename elemT, typename FunctionObjectPtrIter>
 //! 1d specialisation for general function objects
 /*! \ingroup Array
  */
+template <typename elemT, typename FunctionObjectPtrIter>
 inline void apply_array_functions_on_each_index(Array<1, elemT>& out_array,
                                                 const Array<1, elemT>& in_array,
                                                 FunctionObjectPtrIter start,
                                                 FunctionObjectPtrIter stop);
 
+//! \name  apply binary function element-wise
+/*!  \ingroup Array
+  arrays need to have the same number of elements, but can have different shapes.
+*/
+//@{
+//! apply binary function element-wise and store in other array
+template <int num_dim, typename elemTout, typename elemT1, typename elemT2, typename BinaryFunctionT>
+inline void apply_binary_func_element_wise(Array<num_dim, elemTout>& out,
+                                           const Array<num_dim, elemT1>& in1,
+                                           const Array<num_dim, elemT2>& in2,
+                                           BinaryFunctionT f);
+
+//! conditionally apply binary function element-wise and store in other array
+/*!
+  arrays need to have the same number of elements, but can have different shapes.
+
+  Element-wise loop, only computing/storing results if <code>predicate(*in1_full_iter, *in2_full_iter)==true</code>.
+  Other elements in \c out are unassigned.
+*/
+template <int num_dim,
+          typename elemTout,
+          typename elemT1,
+          typename elemT2,
+          typename PredicateBinaryFunctionT,
+          typename BinaryFunctionT>
+inline void apply_binary_func_element_wise(Array<num_dim, elemTout>& out,
+                                           const Array<num_dim, elemT1>& in1,
+                                           const Array<num_dim, elemT2>& in2,
+                                           PredicateBinaryFunctionT predicate,
+                                           BinaryFunctionT f);
+
+//! conditionally apply binary function element-wise and store in other array
+/*!
+  arrays need to have the same number of elements, but can have different shapes.
+
+  Element-wise loop, only computing/storing results if <code>bool(*where_full_iter)==true</code>.
+  Other elements in \c out are unassigned.
+*/
+template <int num_dim, typename elemTout, typename elemT1, typename elemT2, typename elemT3, typename BinaryFunctionT>
+inline void apply_binary_func_element_wise(Array<num_dim, elemTout>& out,
+                                           const Array<num_dim, elemT1>& in1,
+                                           const Array<num_dim, elemT2>& in2,
+                                           const Array<num_dim, elemT3>& where,
+                                           BinaryFunctionT f);
+
+template <int num_dim, typename elemTinout, typename elemT2, typename BinaryFunctionT>
+inline void
+in_place_apply_binary_func_element_wise(Array<num_dim, elemTinout>& inout, const Array<num_dim, elemT2>& in2, BinaryFunctionT f);
+
+template <int num_dim, typename elemTinout, typename elemT2, typename PredicateBinaryFunctionT, typename BinaryFunctionT>
+inline void in_place_apply_binary_func_element_wise(Array<num_dim, elemTinout>& inout,
+                                                    const Array<num_dim, elemT2>& in2,
+                                                    PredicateBinaryFunctionT predicate,
+                                                    BinaryFunctionT f);
+
+template <int num_dim, typename elemTinout, typename elemT2, typename elemT3, typename BinaryFunctionT>
+inline void in_place_apply_binary_func_element_wise(Array<num_dim, elemTinout>& inout,
+                                                    const Array<num_dim, elemT2>& in2,
+                                                    const Array<num_dim, elemT3>& where,
+                                                    BinaryFunctionT f);
+//@}
+
 template <int num_dim, typename elemT>
 inline void transform_array_to_periodic_indices(Array<num_dim, elemT>& out_array, const Array<num_dim, elemT>& in_array);
 template <int num_dim, typename elemT>

src/include/stir/ArrayFunction.inl

@@ -1,6 +1,7 @@
 /*
     Copyright (C) 2000 PARAPET partners
     Copyright (C) 2000- 2007, Hammersmith Imanet Ltd
+    Copyright (C) 2024, University College London
     This file is part of STIR.
 
     SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license
@@ -16,24 +17,14 @@
   \author Kris Thielemans (some functions based on some earlier work by Darren Hague)
   \author PARAPET project
 
-
-  \warning Compilers without partial specialisation of templates are
-   catered for by explicit instantiations. If you need it for any other
-   types, you'd have to add them by hand.
  */
 #include "stir/BasicCoordinate.h"
 #include "stir/array_index_functions.h"
 #include "stir/modulo.h"
 
 #include <cmath>
 #include <complex>
-#ifdef BOOST_NO_STDC_NAMESPACE
-namespace std
-{
-using ::log;
-using ::exp;
-} // namespace std
-#endif
+#include <algorithm>
 
 START_NAMESPACE_STIR
 
@@ -110,6 +101,122 @@ in_place_apply_function(T& v, FUNCTION f)
   return v;
 }
 
+template <int num_dim, typename elemTout, typename elemT1, typename elemT2, typename BinaryFunctionT>
+inline void
+apply_binary_func_element_wise(Array<num_dim, elemTout>& out,
+                               const Array<num_dim, elemT1>& in1,
+                               const Array<num_dim, elemT2>& in2,
+                               BinaryFunctionT f)
+{
+  std::transform(in1.begin_all(), in1.end_all(), in2.begin_all(), out.begin_all(), f);
+}
+
+template <int num_dim,
+          typename elemTout,
+          typename elemT1,
+          typename elemT2,
+          typename PredicateBinaryFunctionT,
+          typename BinaryFunctionT>
+inline void
+apply_binary_func_element_wise(Array<num_dim, elemTout>& out,
+                               const Array<num_dim, elemT1>& in1,
+                               const Array<num_dim, elemT2>& in2,
+                               PredicateBinaryFunctionT predicate,
+                               BinaryFunctionT f)
+{
+  auto in1_iter = in1.begin_all();
+  const auto in1_end = in1.end_all();
+  auto in2_iter = in2.begin_all();
+  auto out_iter = out.begin_all();
+  while (in1_iter != in1_end)
+    {
+      if (predicate(*in1_iter, *in2_iter))
+        *out_iter = f(*in1_iter, *in2_iter);
+      ++in1_iter;
+      ++in2_iter;
+      ++out_iter;
+    }
+}
+
+template <int num_dim, typename elemTout, typename elemT1, typename elemT2, typename elemT3, typename BinaryFunctionT>
+inline void
+apply_binary_func_element_wise(Array<num_dim, elemTout>& out,
+                               const Array<num_dim, elemT1>& in1,
+                               const Array<num_dim, elemT2>& in2,
+                               const Array<num_dim, elemT3>& where,
+                               BinaryFunctionT f)
+{
+  auto in1_iter = in1.begin_all();
+  const auto in1_end = in1.end_all();
+  auto in2_iter = in2.begin_all();
+  auto out_iter = out.begin_all();
+  auto where_iter = where.begin_all();
+  while (in1_iter != in1_end)
+    {
+      if (*where_iter)
+        *out_iter = f(*in1_iter, *in2_iter);
+      ++in1_iter;
+      ++in2_iter;
+      ++where_iter;
+      ++out_iter;
+    }
+}
+
+template <int num_dim, typename elemTinout, typename elemT2, typename BinaryFunctionT>
+inline void
+in_place_apply_binary_func_element_wise(Array<num_dim, elemTinout>& inout, const Array<num_dim, elemT2>& in2, BinaryFunctionT f)
+{
+  auto inout_iter = inout.begin_all();
+  const auto inout_end = inout.end_all();
+  auto in2_iter = in2.begin_all();
+  while (inout_iter != inout_end)
+    {
+      *inout_iter = f(*inout_iter, *in2_iter);
+      ++inout_iter;
+      ++in2_iter;
+    }
+}
+
+template <int num_dim, typename elemTinout, typename elemT2, typename PredicateBinaryFunctionT, typename BinaryFunctionT>
+inline void
+in_place_apply_binary_func_element_wise(Array<num_dim, elemTinout>& inout,
+                                        const Array<num_dim, elemT2>& in2,
+                                        PredicateBinaryFunctionT predicate,
+                                        BinaryFunctionT f)
+{
+  auto inout_iter = inout.begin_all();
+  const auto inout_end = inout.end_all();
+  auto in2_iter = in2.begin_all();
+  while (inout_iter != inout_end)
+    {
+      if (predicate(*inout_iter, *in2_iter))
+        *inout_iter = f(*inout_iter, *in2_iter);
+      ++inout_iter;
+      ++in2_iter;
+    }
+}
+
+template <int num_dim, typename elemTinout, typename elemT2, typename elemT3, typename BinaryFunctionT>
+inline void
+in_place_apply_binary_func_element_wise(Array<num_dim, elemTinout>& inout,
+                                        const Array<num_dim, elemT2>& in2,
+                                        const Array<num_dim, elemT3>& where,
+                                        BinaryFunctionT f)
+{
+  auto inout_iter = inout.begin_all();
+  const auto inout_end = inout.end_all();
+  auto in2_iter = in2.begin_all();
+  auto where_iter = where.begin_all();
+  while (inout_iter != inout_end)
+    {
+      if (*where_iter)
+        *inout_iter = f(*inout_iter, *in2_iter);
+      ++inout_iter;
+      ++in2_iter;
+      ++where_iter;
+    }
+}
+
 template <int num_dim, typename elemT, typename FunctionObjectPtr>
 inline void
 in_place_apply_array_function_on_1st_index(Array<num_dim, elemT>& array, FunctionObjectPtr f)

src/include/stir/OSMAPOSL/OSMAPOSLReconstruction.h

@@ -150,12 +150,6 @@ class OSMAPOSLReconstruction
   //! determines whether non-positive values in the initial image will be set to small positive ones
   bool enforce_initial_positivity;
 
-  //! determines wether voxels outside a circular FOV will be set to 0 or not
-  /*! Currently this circular FOV is slightly smaller than the actual image size (1 pixel at each end or so).
-      \deprecated
-  */
-  bool do_rim_truncation;
-
   //! subiteration interval at which to apply inter-update filters
   int inter_update_filter_interval;
 

src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx

@@ -27,6 +27,7 @@
 #include "stir/recon_buildblock/TrivialBinNormalisation.h"
 #include "stir/Succeeded.h"
 #include "stir/RelatedViewgrams.h"
+#include "stir/ArrayFunction.h"
 #include "stir/stream.h"
 #include "stir/info.h"
 #include "stir/warning.h"
@@ -780,49 +781,6 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjData<TargetT>::actual_compute_o
   return accum;
 }
 
-#if 0
-template<typename TargetT>
-float 
-PoissonLogLikelihoodWithLinearModelForMeanAndProjData<TargetT>::
-sum_projection_data() const
-{
-
-  float counts=0.0F;
-
-  for (int segment_num = -max_segment_num_to_process; segment_num <= max_segment_num_to_process; ++segment_num)
-  {
-	  for (int timing_pos_num = -max_timing_pos_num_to_process; timing_pos_num <= max_timing_pos_num_to_process; ++timing_pos_num)
-	  {
-		for (int view_num = proj_data_sptr->get_min_view_num();
-			 view_num <= proj_data_sptr->get_max_view_num();
-			 ++view_num)
-		{
-
-		  Viewgram<float>  viewgram=proj_data_sptr->get_viewgram(view_num,segment_num,false,timing_pos_num);
-
-		  //first adjust data
-
-		  // KT 05/07/2000 made parameters.zero_seg0_end_planes int
-		  if(segment_num==0 && zero_seg0_end_planes!=0)
-		  {
-			viewgram[viewgram.get_min_axial_pos_num()].fill(0);
-			viewgram[viewgram.get_max_axial_pos_num()].fill(0);
-		  }
-
-		  truncate_rim(viewgram,rim_truncation_sino);
-
-		  //now take totals
-		  counts+=viewgram.sum();
-		}
-	  }
-  }
-
-  return counts;
-
-}
-
-#endif
-
 template <typename TargetT>
 void
 PoissonLogLikelihoodWithLinearModelForMeanAndProjData<TargetT>::add_subset_sensitivity(TargetT& sensitivity,
@@ -1183,11 +1141,22 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjData<TargetT>::actual_accumulat
       // final_viewgram starts as measured data
       RelatedViewgrams<float> final_viewgram = this->get_proj_data().get_related_viewgrams(vg_idx_to_process[i], symmetries_sptr);
       {
+#if 0
         // Mult input_viewgram
         final_viewgram *= input_viewgrams_vec[i];
+        // Divide final_viewgram by ybar_sq_viewgram
         int tmp1 = 0, tmp2 = 0; // ignore counters returned by divide_and_truncate
-        // Divide final_viewgeam by ybar_sq_viewgram
         divide_and_truncate(final_viewgram, ybar_sq_viewgram, 0, tmp1, tmp2);
+#else
+        // use division where 0/x = 0, but truncate measured data to be >= 0 (as in the value and gradient)
+        constexpr auto div = [](float a, float b) { return a <= 0.F ? 0.F : a / b; };
+        auto f_iter = final_viewgram.begin();
+        auto yb_sq_iter = ybar_sq_viewgram.begin();
+        for (; f_iter != final_viewgram.end(); ++f_iter, ++yb_sq_iter)
+          in_place_apply_binary_func_element_wise(*f_iter, *yb_sq_iter, div);
+        // Mult input_viewgram
+        final_viewgram *= input_viewgrams_vec[i];
+#endif
       }
 
       // back-project final_viewgram