fix: handling of non-finite values in background histogram binning

SEFramework/src/lib/Psf/VariablePsf.cpp

@@ -22,9 +22,9 @@
  */
 
 #include <algorithm>
+#include <cmath>
 #include <ElementsKernel/Logging.h>
 #include <ElementsKernel/Exception.h>
-#include "SEUtils/IsNan.h"
 
 #include "SEFramework/Psf/VariablePsf.h"
 
@@ -136,10 +136,10 @@ void VariablePsf::selfTest() {
     }
     for (auto x = 0; x < psf_width; ++x){
       for (auto y = 0; y < psf_height; ++y) {
-        if (fastmath_isnan(coeff->at(x, y))) {
+        if (std::isnan(coeff->at(x, y))) {
           throw Elements::Exception() << "Malformed variable PSF, coefficient matrices contains NANs";
         }
-        else if (fastmath_isinf(coeff->at(x, y))){
+        else if (std::isinf(coeff->at(x, y))){
           throw Elements::Exception() << "Malformed variable PSF, coefficient matrices contains INFs";
         }
       }

SEImplementation/CMakeLists.txt

@@ -129,6 +129,9 @@ elements_add_python_module(SourceXtractorPy
 #                       INCLUDE_DIRS ElementsExamples
 #                       LINK_LIBRARIES ElementsExamples TYPE Boost)
 #===============================================================================
+elements_add_unit_test(ImageMode_test tests/src/Background/ImageMode_test.cpp
+		             LINK_LIBRARIES SEImplementation
+		             TYPE Boost)
 elements_add_unit_test(KappaSigmaBinning_test tests/src/Background/KappaSigmaBinning_test.cpp
 		             LINK_LIBRARIES SEImplementation
 		             TYPE Boost)

SEImplementation/SEImplementation/Background/SE/KappaSigmaBinning.h

@@ -19,6 +19,9 @@
 #define SOURCEXTRACTORPLUSPLUS_KAPPASIGMABINNING_H
 
 #include <iostream>
+#include <iterator>
+
+#include "ElementsKernel/Exception.h"
 #include "Histogram/Histogram.h"
 
 namespace SourceXtractor {
@@ -67,6 +70,11 @@ class KappaSigmaBinning: public Euclid::Histogram::BinStrategy<VarType> {
    */
   template<typename Iterator>
   void computeBins(Iterator begin, Iterator end) {
+    // *begin is dereferenced below, and there are no bins to compute anyway
+    if (begin == end) {
+      throw Elements::Exception() << "Can't compute background histogram bins: no data points";
+    }
+
     // Compute mean and standard deviation of the original data set
     double mean, sigma;
     size_t ndata;
@@ -95,7 +103,15 @@ class KappaSigmaBinning: public Euclid::Histogram::BinStrategy<VarType> {
     }
     std::tie(mean, sigma, ndata) = stats.get();
 
-    assert(ndata > 0);
+    // This case shouldn't be reachable anymore with fixes in #649, but still
+    // throw here just in case (previously was a bare assert); this way at least
+    // it gives a more informative error message.
+    if (ndata == 0) {
+      throw Elements::Exception()
+          << "Can't compute background histogram bins: none of the " << std::distance(begin, end)
+          << " input values lie within " << m_kappa << " sigma of their mean (" << mean
+          << "), which should not be possible for finite data";
+    }
 
     // Number of bins
     m_nbins = computeBinCount(ndata);
@@ -133,30 +149,47 @@ class KappaSigmaBinning: public Euclid::Histogram::BinStrategy<VarType> {
 
   /**
    * Online mean and standard deviation computation
-   * Similar to boost::accumulators, but this turned out to be faster
    */
   struct Stats {
-    double mean = 0, sigma = 0;
+    /// Everything is accumulated as an offset from the first value seen, so that the
+    /// sums stay on the scale of the spread rather than that of the data
+    double shift = 0, sum = 0, sum_sq = 0;
     size_t ndata = 0;
 
+    // Shifting by a value from the data largely avoids the catastrophic cancellation
+    // that E[X^2] - E[X]^2 suffers when the spread is small next to the magnitude.
+    // See https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance, under
+    // "Computing shifted data"; Welford's algorithm there fixes it too, but costs
+    // about 4x as much, as it performs a division per sample.
     void operator() (VarType v) {
-      mean += v;
-      sigma += v * v;
+      if (ndata == 0) {
+        shift = v;
+      }
+
+      double offset = v - shift;
+      sum += offset;
+      sum_sq += offset * offset;
       ++ndata;
     }
 
     std::tuple<double, double, size_t> get() {
-      mean /= ndata;
-      sigma = sigma / ndata - mean * mean;
+      if (ndata == 0) {
+        return std::make_tuple(0., 0., ndata);
+      }
+
+      double mean = shift + sum / ndata;
+      double sigma = (sum_sq - sum * sum / ndata) / ndata;
+
       if (sigma > std::numeric_limits<double>::epsilon())
         sigma = std::sqrt(sigma);
       else
         sigma = 0.;
+
       return std::make_tuple(mean, sigma, ndata);
     }
 
     void reset() {
-      mean = sigma = 0;
+      shift = sum = sum_sq = 0;
       ndata = 0;
     }
   };

SEImplementation/src/lib/Background/SE/ImageMode.cpp

@@ -15,8 +15,12 @@
  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  */
 
+#include <cmath>
+
 #include <Histogram/Histogram.h> // From Alexandria
 
+#include "ElementsKernel/Exception.h"
+
 #include "SEFramework/Image/ImageChunk.h"
 #include "SEImplementation/Background/SE/ImageMode.h"
 #include "SEImplementation/Background/SE/KappaSigmaBinning.h"
@@ -34,7 +38,8 @@ ImageMode<T>::ImageMode(const std::shared_ptr<Image<T>>& image, const std::share
                                                             m_cell_w(cell_w), m_cell_h(cell_h),
                                                             m_invalid(invalid_value),
                                                             m_kappa1(kappa1), m_kappa2(kappa2), m_kappa3(kappa3),
-                                                            m_rtol(rtol), m_max_iter(max_iter), m_pearson_factor(pearson_factor) {
+                                                            m_rtol(rtol), m_pearson_factor(pearson_factor),
+                                                            m_max_iter(max_iter) {
   auto hist_width = std::div(image->getWidth(), m_cell_w);
   if (hist_width.rem)
     ++hist_width.quot;
@@ -96,7 +101,10 @@ std::tuple<T, T> ImageMode<T>::getBackGuess(const std::vector<T>& data) const {
   std::tie(mean, median, sigma) = histo.getStats();
   T prev_sigma = sigma * 10;
 
-  assert(!std::isnan(mean));
+  // processCell() filters non-finite pixels out, so this should be unreachable
+  if (std::isnan(mean)) {
+    throw Elements::Exception() << "Background histogram produced a NaN mean";
+  }
 
   for (size_t iter = 0; iter < m_max_iter && sigma > atol && std::abs(sigma / prev_sigma - 1.0) > m_rtol; ++iter) {
     histo.clip(median - sigma * m_kappa3, median + sigma * m_kappa3);
@@ -134,10 +142,13 @@ void ImageMode<T>::processCell(const Image<T>& img, int x, int y,
   std::vector<T> filtered;
   filtered.reserve(w * h);
 
-  for (int y = 0; y < h; ++y) {
-    for (int x = 0; x < w; ++x) {
-      auto v = img_chunk.getValue(x, y);
-      if (v != m_invalid)
+  for (int py = 0; py < h; ++py) {
+    for (int px = 0; px < w; ++px) {
+      auto v = img_chunk.getValue(px, py);
+      // A NaN compares unequal to everything, so the sentinel test alone lets it
+      // through, and one is enough to poison the whole cell's statistics.  An
+      // infinity does the same by way of Inf - Inf
+      if (v != m_invalid && !std::isnan(v) && !std::isinf(v))
         filtered.emplace_back(v);
     }
   }

SEImplementation/tests/src/Background/ImageMode_test.cpp

@@ -0,0 +1,161 @@
+/** Copyright © 2019-2026 Université de Genève, LMU Munich - Faculty of Physics, IAP-CNRS/Sorbonne Université
+ *
+ * This library is free software; you can redistribute it and/or modify it under
+ * the terms of the GNU Lesser General Public License as published by the Free
+ * Software Foundation; either version 3.0 of the License, or (at your option)
+ * any later version.
+ *
+ * This library is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this library; if not, write to the Free Software Foundation, Inc.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#include <limits>
+#include <tuple>
+#include <vector>
+
+#include <boost/test/unit_test.hpp>
+
+#include "SEFramework/Image/VectorImage.h"
+#include "SEImplementation/Background/SE/ImageMode.h"
+
+using namespace SourceXtractor;
+
+namespace {
+
+constexpr int CELL = 10;
+const float INVALID = std::numeric_limits<float>::lowest();
+
+/**
+ * One cell's worth of pixels: a repeating ramp over [100.0, 101.6]
+ *
+ * The values have to be spread rather than constant, since a constant cell gives
+ * sigma == 0, which takes a different branch in both computeBins() and
+ * getBackGuess() than the one these tests are about.
+ */
+std::vector<float> goodCell() {
+  std::vector<float> values;
+  values.reserve(CELL * CELL);
+
+  for (int i = 0; i < CELL * CELL; i++) {
+    values.push_back(100.0f + static_cast<float>((i * 3) % 17) * 0.1f);
+  }
+
+  return values;
+}
+
+/// The mode and sigma of a single cell whose pixel at `index` is set to `value`
+std::tuple<float, float> modeWithPixel(size_t index, float value) {
+  auto values = goodCell();
+  values[index] = value;
+
+  ImageMode<float> histo(VectorImage<float>::create(CELL, CELL, values), nullptr, CELL, CELL,
+                         INVALID);
+
+  return std::make_tuple(histo.getModeImage()->getValue(0, 0),
+                         histo.getSigmaImage()->getValue(0, 0));
+}
+
+}  // anonymous namespace
+
+//-----------------------------------------------------------------------------
+
+BOOST_AUTO_TEST_SUITE(ImageMode_test)
+
+//-----------------------------------------------------------------------------
+
+/**
+ * A non-finite pixel must be discarded exactly like the invalid sentinel
+ *
+ * A NaN compares unequal to everything, so testing only against the sentinel let
+ * it into the cell's statistics, where a single one was enough to make the mean
+ * NaN, every histogram cut NaN, and the clipped selection empty.  An infinity
+ * does the same by way of Inf - Inf.
+ *
+ * Comparing against the sentinel is the sharp form of the check: all of them
+ * discard the same pixel, so all of them must agree exactly.
+ */
+BOOST_AUTO_TEST_CASE(non_finite_pixel_discarded_like_invalid) {
+  const size_t pixel = 42;
+
+  float expected_mode, expected_sigma;
+  std::tie(expected_mode, expected_sigma) = modeWithPixel(pixel, INVALID);
+
+  // a sanity check on the fixture itself: the good cell must produce a usable
+  // background, or the comparisons below would hold vacuously
+  BOOST_CHECK_NE(expected_mode, INVALID);
+
+  for (float bad : {std::numeric_limits<float>::quiet_NaN(),
+                    std::numeric_limits<float>::infinity(),
+                    -std::numeric_limits<float>::infinity()}) {
+    float mode, sigma;
+    std::tie(mode, sigma) = modeWithPixel(pixel, bad);
+    BOOST_CHECK_EQUAL(expected_mode, mode);
+    BOOST_CHECK_EQUAL(expected_sigma, sigma);
+  }
+}
+
+//-----------------------------------------------------------------------------
+
+/**
+ * A cell with no usable pixels is marked invalid rather than failing
+ *
+ * This is the path SEBackgroundLevelAnalyzer's ReplaceUndef interpolates over,
+ * and it is what a mostly-empty mosaic relies on: a resampled survey skycell,
+ * say, which is NaN everywhere outside its footprint (similar to a real world
+ * case from Roman which tripped this).
+ */
+BOOST_AUTO_TEST_CASE(all_nan_cell_is_invalid) {
+  std::vector<float> values(CELL * CELL, std::numeric_limits<float>::quiet_NaN());
+
+  ImageMode<float> histo(VectorImage<float>::create(CELL, CELL, values), nullptr, CELL, CELL,
+                         INVALID);
+
+  BOOST_CHECK_EQUAL(INVALID, histo.getModeImage()->getValue(0, 0));
+  BOOST_CHECK_EQUAL(INVALID, histo.getSigmaImage()->getValue(0, 0));
+}
+
+//-----------------------------------------------------------------------------
+
+/// A mostly-NaN cell is invalid too, by the pre-existing "at least half usable" rule
+BOOST_AUTO_TEST_CASE(mostly_nan_cell_is_invalid) {
+  auto values = goodCell();
+
+  for (size_t i = 0; i < 60; i++) {
+    values[i] = std::numeric_limits<float>::quiet_NaN();
+  }
+
+  ImageMode<float> histo(VectorImage<float>::create(CELL, CELL, values), nullptr, CELL, CELL,
+                         INVALID);
+
+  BOOST_CHECK_EQUAL(INVALID, histo.getModeImage()->getValue(0, 0));
+}
+
+//-----------------------------------------------------------------------------
+
+/// A cell that is NaN only in part still yields a background, from what is left
+BOOST_AUTO_TEST_CASE(partially_nan_cell_is_usable) {
+  auto values = goodCell();
+
+  for (size_t i = 0; i < 30; i++) {
+    values[i] = std::numeric_limits<float>::quiet_NaN();
+  }
+
+  ImageMode<float> histo(VectorImage<float>::create(CELL, CELL, values), nullptr, CELL, CELL,
+                         INVALID);
+
+  auto mode = histo.getModeImage()->getValue(0, 0);
+  BOOST_CHECK_NE(mode, INVALID);
+  // the surviving pixels all sit within [100, 101.6]
+  BOOST_CHECK_GT(mode, 99.0f);
+  BOOST_CHECK_LT(mode, 102.0f);
+}
+
+//-----------------------------------------------------------------------------
+
+BOOST_AUTO_TEST_SUITE_END()

SEImplementation/tests/src/Background/KappaSigmaBinning_test.cpp

@@ -218,6 +218,39 @@ BOOST_AUTO_TEST_CASE(compat4Test) {
 
 //-----------------------------------------------------------------------------
 
+/// There are no bins to compute for an empty range, and *begin used to be
+/// dereferenced unconditionally
+BOOST_AUTO_TEST_CASE(emptyRange) {
+  std::vector<float> empty;
+  KappaSigmaBinning<float> binning;
+  BOOST_CHECK_THROW(binning.computeBins(empty.begin(), empty.end()), Elements::Exception);
+}
+
+//-----------------------------------------------------------------------------
+
+/// Every comparison against a NaN cut is false, so nothing is selected and no
+/// bins can be computed.  This has to throw an exception, not a bare assert:
+/// otherwise release builds with -ffast-math allowed NaN to propagate silently
+/// into the background model
+BOOST_AUTO_TEST_CASE(allNaN) {
+  std::vector<float> nans(8, std::numeric_limits<float>::quiet_NaN());
+  KappaSigmaBinning<float> binning;
+  BOOST_CHECK_THROW(binning.computeBins(nans.begin(), nans.end()), Elements::Exception);
+}
+
+//-----------------------------------------------------------------------------
+
+/// Likewise for a single NaN among otherwise good values: it poisons the mean
+/// and sigma, so the cuts become NaN and select nothing
+BOOST_FIXTURE_TEST_CASE(singleNaN, KappaSigmaBinningFixture) {
+  auto values = data;
+  values[5] = std::numeric_limits<float>::quiet_NaN();
+  KappaSigmaBinning<float> binning;
+  BOOST_CHECK_THROW(binning.computeBins(values.begin(), values.end()), Elements::Exception);
+}
+
+//-----------------------------------------------------------------------------
+
 BOOST_AUTO_TEST_SUITE_END()
 
 //-----------------------------------------------------------------------------