From b18e889b8a708786b72c9324896e9ef0f5a33062 Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Thu, 15 Oct 2026 13:51:05 +0200 Subject: [PATCH 1/3] fix: handling of non-finite values in background histogram binning Two bugs, both the symptom and the cause. The cause was insufficient NaN comparison in ImageMode::getBackGuess--on debug builds a bare assert(!isnan(mean)) caught mean=NaN, but on release builds this fell through. Because of -ffast-math the NaN was simply allowed to propagate through the statistics. In KappaSigmaBinning::computeBins this NaN propagation also meant that a NaN could poison the whole cell, resulting in ``ndata==0``. This was also asserted against on debug builds, but allowed to propagate on release builds. Uses the fastmath_isnan and fastmath_isinf from SEUtils since these were introduced exactly for this reason, that we still need to sweep for non-finite values even of `-ffast-math` builds which are created by default on release. --- SEImplementation/CMakeLists.txt | 3 + .../Background/SE/KappaSigmaBinning.h | 18 +- .../src/lib/Background/SE/ImageMode.cpp | 15 +- .../tests/src/Background/ImageMode_test.cpp | 161 ++++++++++++++++++ .../src/Background/KappaSigmaBinning_test.cpp | 33 ++++ 5 files changed, 227 insertions(+), 3 deletions(-) create mode 100644 SEImplementation/tests/src/Background/ImageMode_test.cpp diff --git a/SEImplementation/CMakeLists.txt b/SEImplementation/CMakeLists.txt index 02cc4bafc..35f19a8ed 100644 --- a/SEImplementation/CMakeLists.txt +++ b/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) diff --git a/SEImplementation/SEImplementation/Background/SE/KappaSigmaBinning.h b/SEImplementation/SEImplementation/Background/SE/KappaSigmaBinning.h index 59103e46c..6bae0a033 100644 --- a/SEImplementation/SEImplementation/Background/SE/KappaSigmaBinning.h +++ b/SEImplementation/SEImplementation/Background/SE/KappaSigmaBinning.h @@ -19,6 +19,8 @@ #define SOURCEXTRACTORPLUSPLUS_KAPPASIGMABINNING_H #include + +#include "ElementsKernel/Exception.h" #include "Histogram/Histogram.h" namespace SourceXtractor { @@ -67,6 +69,11 @@ class KappaSigmaBinning: public Euclid::Histogram::BinStrategy { */ template 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 +102,16 @@ class KappaSigmaBinning: public Euclid::Histogram::BinStrategy { } std::tie(mean, sigma, ndata) = stats.get(); - assert(ndata > 0); + // Every comparison against a NaN cut is false, so an empty selection here is + // what a non-finite input value looks like by the time it reaches this point. + // Throw in this case (previously there was a bare assert which, when compiled + // out, would just fall down to the code below allowing NaNs to propagate + // silently, or eventually segfault). + if (ndata == 0) { + throw Elements::Exception() + << "Can't compute background histogram bins: no data points within " << m_kappa + << " sigma of the mean; the input may contain NaN or infinite values"; + } // Number of bins m_nbins = computeBinCount(ndata); diff --git a/SEImplementation/src/lib/Background/SE/ImageMode.cpp b/SEImplementation/src/lib/Background/SE/ImageMode.cpp index 65ed307f8..259e6161d 100644 --- a/SEImplementation/src/lib/Background/SE/ImageMode.cpp +++ b/SEImplementation/src/lib/Background/SE/ImageMode.cpp @@ -17,9 +17,12 @@ #include // From Alexandria +#include "ElementsKernel/Exception.h" + #include "SEFramework/Image/ImageChunk.h" #include "SEImplementation/Background/SE/ImageMode.h" #include "SEImplementation/Background/SE/KappaSigmaBinning.h" +#include "SEUtils/IsNan.h" using Euclid::Histogram::Histogram; @@ -96,7 +99,11 @@ std::tuple ImageMode::getBackGuess(const std::vector& 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; + // std::isnan would have been compiled away by -ffast-math in any case + if (fastmath_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); @@ -137,7 +144,11 @@ void ImageMode::processCell(const Image& img, int x, int y, 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) + // 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. fastmath_isnan/isinf rather + // than std::isnan/isfinite because -ffast-math compiles the latter away + if (v != m_invalid && !fastmath_isnan(v) && !fastmath_isinf(v)) filtered.emplace_back(v); } } diff --git a/SEImplementation/tests/src/Background/ImageMode_test.cpp b/SEImplementation/tests/src/Background/ImageMode_test.cpp new file mode 100644 index 000000000..009973a6a --- /dev/null +++ b/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 +#include +#include + +#include + +#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::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 goodCell() { + std::vector values; + values.reserve(CELL * CELL); + + for (int i = 0; i < CELL * CELL; i++) { + values.push_back(100.0f + static_cast((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 modeWithPixel(size_t index, float value) { + auto values = goodCell(); + values[index] = value; + + ImageMode histo(VectorImage::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::quiet_NaN(), + std::numeric_limits::infinity(), + -std::numeric_limits::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 values(CELL * CELL, std::numeric_limits::quiet_NaN()); + + ImageMode histo(VectorImage::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::quiet_NaN(); + } + + ImageMode histo(VectorImage::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::quiet_NaN(); + } + + ImageMode histo(VectorImage::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() diff --git a/SEImplementation/tests/src/Background/KappaSigmaBinning_test.cpp b/SEImplementation/tests/src/Background/KappaSigmaBinning_test.cpp index b71e94dbe..5c382e6ba 100644 --- a/SEImplementation/tests/src/Background/KappaSigmaBinning_test.cpp +++ b/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 empty; + KappaSigmaBinning 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 nans(8, std::numeric_limits::quiet_NaN()); + KappaSigmaBinning 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::quiet_NaN(); + KappaSigmaBinning binning; + BOOST_CHECK_THROW(binning.computeBins(values.begin(), values.end()), Elements::Exception); +} + +//----------------------------------------------------------------------------- + BOOST_AUTO_TEST_SUITE_END() //----------------------------------------------------------------------------- From 84ba73dc5d848a689dce2ffb228f79c51bd8b4c2 Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Thu, 22 Oct 2026 16:23:48 +0200 Subject: [PATCH 2/3] fix: use variance algorithm not as succeptible to catastrophic cancellation This turned out to be a more serious issue affecting this code on nearly flat background cells with close adjacent values, resulting in some caes in catastrophic cancellation. I asked Claude to benchmark different alternative algorithms, namely shifted data and Welford's, against similar data that reproduced the problem, both for speed and accuracy. There was also a comment in the code that suggested it was faster that boost::accumulator. It turns out that actually has two algorithms for calculating the variance: `accumlators::lazy_variance` and `accumulators::variance`. That claim was probably against the "lazy" variant, whereas the non-lazy variant is essentially Weldford's. In any case the claim still seems to hold that this is a bit faster. Additionally I realized that the exception message thrown (replacing the previous bare assert()) was maybe not accurate since from the other fixes we should not be getting non-finite values in this code in the first place. They should already be filtered out, and the exception is almost unreachable code otherwise (with this fix). Leave it in anyways for future catches like this: it's still a more helpful error message than just crashing as before. --- .../Background/SE/KappaSigmaBinning.h | 45 +++++++++++++------ 1 file changed, 31 insertions(+), 14 deletions(-) diff --git a/SEImplementation/SEImplementation/Background/SE/KappaSigmaBinning.h b/SEImplementation/SEImplementation/Background/SE/KappaSigmaBinning.h index 6bae0a033..c34498445 100644 --- a/SEImplementation/SEImplementation/Background/SE/KappaSigmaBinning.h +++ b/SEImplementation/SEImplementation/Background/SE/KappaSigmaBinning.h @@ -19,6 +19,7 @@ #define SOURCEXTRACTORPLUSPLUS_KAPPASIGMABINNING_H #include +#include #include "ElementsKernel/Exception.h" #include "Histogram/Histogram.h" @@ -102,15 +103,14 @@ class KappaSigmaBinning: public Euclid::Histogram::BinStrategy { } std::tie(mean, sigma, ndata) = stats.get(); - // Every comparison against a NaN cut is false, so an empty selection here is - // what a non-finite input value looks like by the time it reaches this point. - // Throw in this case (previously there was a bare assert which, when compiled - // out, would just fall down to the code below allowing NaNs to propagate - // silently, or eventually segfault). + // 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: no data points within " << m_kappa - << " sigma of the mean; the input may contain NaN or infinite values"; + << "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 @@ -149,30 +149,47 @@ class KappaSigmaBinning: public Euclid::Histogram::BinStrategy { /** * 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 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::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; } }; From d87bc01cfc35d73f824c7615b12312efd0aa9760 Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Tue, 27 Oct 2026 17:37:13 +0100 Subject: [PATCH 3/3] fix: disable -ffast-math-only by default on all builds - As discussed in gh-649 this results in buggy, ill-defined behavior - Can still be enabled via -DFINITE_MATH_ONLY=ON if desired, or for comparison - Deletes the IsNan.h workaround which was not a reliable fix anyways, esp. as observed on recent clang versions --- SEFramework/src/lib/Psf/VariablePsf.cpp | 6 +- .../src/lib/Background/SE/ImageMode.cpp | 22 ++--- .../Segmentation/LutzSegmentation_test.cpp | 62 ++++++++++++ SEUtils/CMakeLists.txt | 3 + SEUtils/SEUtils/IsNan.h | 80 --------------- SEUtils/tests/src/NonFiniteMath_test.cpp | 98 +++++++++++++++++++ cmake/LocalBuildFlags.cmake | 23 ++++- 7 files changed, 195 insertions(+), 99 deletions(-) delete mode 100644 SEUtils/SEUtils/IsNan.h create mode 100644 SEUtils/tests/src/NonFiniteMath_test.cpp diff --git a/SEFramework/src/lib/Psf/VariablePsf.cpp b/SEFramework/src/lib/Psf/VariablePsf.cpp index 6e66e98a9..9d8cdd2ad 100644 --- a/SEFramework/src/lib/Psf/VariablePsf.cpp +++ b/SEFramework/src/lib/Psf/VariablePsf.cpp @@ -22,9 +22,9 @@ */ #include +#include #include #include -#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"; } } diff --git a/SEImplementation/src/lib/Background/SE/ImageMode.cpp b/SEImplementation/src/lib/Background/SE/ImageMode.cpp index 259e6161d..095e08bad 100644 --- a/SEImplementation/src/lib/Background/SE/ImageMode.cpp +++ b/SEImplementation/src/lib/Background/SE/ImageMode.cpp @@ -15,6 +15,8 @@ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ +#include + #include // From Alexandria #include "ElementsKernel/Exception.h" @@ -22,7 +24,6 @@ #include "SEFramework/Image/ImageChunk.h" #include "SEImplementation/Background/SE/ImageMode.h" #include "SEImplementation/Background/SE/KappaSigmaBinning.h" -#include "SEUtils/IsNan.h" using Euclid::Histogram::Histogram; @@ -37,7 +38,8 @@ ImageMode::ImageMode(const std::shared_ptr>& 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; @@ -99,9 +101,8 @@ std::tuple ImageMode::getBackGuess(const std::vector& data) const { std::tie(mean, median, sigma) = histo.getStats(); T prev_sigma = sigma * 10; - // processCell() filters non-finite pixels out, so this should be unreachable; - // std::isnan would have been compiled away by -ffast-math in any case - if (fastmath_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"; } @@ -141,14 +142,13 @@ void ImageMode::processCell(const Image& img, int x, int y, std::vector 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); + 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. fastmath_isnan/isinf rather - // than std::isnan/isfinite because -ffast-math compiles the latter away - if (v != m_invalid && !fastmath_isnan(v) && !fastmath_isinf(v)) + // infinity does the same by way of Inf - Inf + if (v != m_invalid && !std::isnan(v) && !std::isinf(v)) filtered.emplace_back(v); } } diff --git a/SEImplementation/tests/src/Segmentation/LutzSegmentation_test.cpp b/SEImplementation/tests/src/Segmentation/LutzSegmentation_test.cpp index 51eb44ada..b521c2685 100644 --- a/SEImplementation/tests/src/Segmentation/LutzSegmentation_test.cpp +++ b/SEImplementation/tests/src/Segmentation/LutzSegmentation_test.cpp @@ -25,6 +25,10 @@ #include +#include +#include +#include + #include "SEFramework/Image/VectorImage.h" #include "SEFramework/Image/ConstantImage.h" #include "SEFramework/Source/SimpleSourceFactory.h" @@ -183,5 +187,63 @@ BOOST_FIXTURE_TEST_CASE( lutz_test, LutzFixture ) { //----------------------------------------------------------------------------- +/// The bit pattern of a 32-bit quiet NaN, read through volatile so that no +/// amount of constant folding can decide the answers below at compile time +volatile uint32_t g_quiet_nan_bits = 0x7FC00000u; + +float runtimeNaN() { + uint32_t bits = g_quiet_nan_bits; + float value; + std::memcpy(&value, &bits, sizeof(value)); + return value; +} + +/// A NaN is not above the detection threshold, so blank pixels must not be labelled. +/// Under -ffinite-math-only `value > 0.0` in Lutz::labelImage is undefined for a NaN and +/// Clang answers true, which on a real mosaic swallows a whole unexposed region +/// represented by NaN into a single component. +BOOST_FIXTURE_TEST_CASE( lutz_nan_test, LutzFixture ) { + // Blank the top four rows the way a mosaic blanks its unexposed area, and put a single + // 3x3 source well clear of it + DetectionImage::PixelType nan_value = runtimeNaN(); + + std::vector pixels(100, 0.0); + for (int idx = 0; idx < 40; idx++) { + pixels[idx] = nan_value; + } + for (int y = 6; y <= 8; y++) { + for (int x = 6; x <= 8; x++) { + pixels[y * 10 + x] = 1.0; + } + } + + auto image = VectorImage::create(10, 10, pixels); + + Segmentation segmentation(nullptr); + segmentation.setLabelling(std::make_shared()); + segmentation.Observable::addObserver(source_observer); + + auto detection_frame = std::make_shared(image); + detection_frame->setBackgroundLevel( + ConstantImage::create(image->getWidth(), image->getHeight(), 0), 0.); + detection_frame->setVarianceMap( + ConstantImage::create(image->getWidth(), image->getHeight(), 0.25)); + segmentation.processFrame(detection_frame); + + // Only the 3x3 source, with none of the blank pixels dragged in with it + BOOST_CHECK_EQUAL(source_observer->m_list.size(), 1); + + for (auto& source : source_observer->m_list) { + auto& coords = source.getCoordinateList(); + BOOST_CHECK_EQUAL(coords.size(), 9u); + for (auto& pixel : coords) { + BOOST_CHECK_GE(pixel.m_y, 6); + BOOST_CHECK_GE(pixel.m_x, 6); + } + } +} + +//----------------------------------------------------------------------------- + BOOST_AUTO_TEST_SUITE_END () diff --git a/SEUtils/CMakeLists.txt b/SEUtils/CMakeLists.txt index 8203f66f7..cd91ea259 100644 --- a/SEUtils/CMakeLists.txt +++ b/SEUtils/CMakeLists.txt @@ -67,6 +67,9 @@ elements_add_unit_test(QuadTree_test tests/src/QuadTree_test.cpp elements_add_unit_test(KdTree_test tests/src/KdTree_test.cpp LINK_LIBRARIES SEUtils TYPE Boost) +elements_add_unit_test(NonFiniteMath_test tests/src/NonFiniteMath_test.cpp + LINK_LIBRARIES SEUtils + TYPE Boost) if(GMOCK_FOUND) elements_add_unit_test(Observable_test tests/src/Observable_test.cpp diff --git a/SEUtils/SEUtils/IsNan.h b/SEUtils/SEUtils/IsNan.h deleted file mode 100644 index c4dea6479..000000000 --- a/SEUtils/SEUtils/IsNan.h +++ /dev/null @@ -1,80 +0,0 @@ -#ifndef _SEUTILS_ISNAN_H -#define _SEUTILS_ISNAN_H - -#include - -namespace SourceXtractor { - -/* - * From: - * https://github.com/searchivarius/BlogCode/tree/2d947d8c42ab0e10f1f28a2ad036eee10389997a/2017/6/1 - * - * IEEE 754 compliant simple functions to test for NANs and INFs. - * This functions are necessary, b/c isnan doesn't work with -Ofast -ffast-math flags - * see http://searchivarius.org/blog/gcc-disables-isnan-and-isinf-when-compiling-ffast-math-flag - * - * See also test files. The file ./regular is compiled without -ffast-math flag and it - * checks for a large number of values that the output of my functions is the same - * as the output of standard functions. For single precision numbers, i.e., for floats - * these checks are exhaustive. That is, I go over the set of all 4B+ possible float values. - * For doubles, this is not possible, so I use tests where the lower 32 bits of mantissa - * are set to zero. - * - * Copyright (c) 2017 Leonid Boytsov - * - * This code is released under the - * Apache License Version 2.0 http://www.apache.org/licenses/. - * - */ - -// A mask to extract an exponent from the single-precision floating point number -// 01111111100000000000000000000000 -const unsigned FLOAT_EXP_MASK = 0x7F800000; -// A mask to extract a mantissa/fractional part from the single-precision floating point number -// 00000000011111111111111111111111 -const unsigned FLOAT_FRAC_PART_MASK = 0x7FFFFF; - -inline bool fastmath_isnan(float x) { - union { - uint32_t u; - float f; - } conv; - conv.f = x; - return ((conv.u & FLOAT_EXP_MASK) == FLOAT_EXP_MASK) && ((conv.u & FLOAT_FRAC_PART_MASK) != 0); -}; - -inline bool fastmath_isinf(float x) { - union { - uint32_t u; - float f; - } conv; - conv.f = x; - return ((conv.u & FLOAT_EXP_MASK) == FLOAT_EXP_MASK) && ((conv.u & FLOAT_FRAC_PART_MASK) == 0); -}; - -// 0111 1111 1111 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 -const uint64_t DOUBLE_EXP_MASK = 0x7FF0000000000000ul; -// 0000 0000 0000 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 -const uint64_t DOUBLE_FRAC_PART_MASK = 0x000FFFFFFFFFFFFFul; - -inline bool fastmath_isnan(double x) { - union { - uint64_t u; - double f; - } conv; - conv.f = x; - return ((conv.u & DOUBLE_EXP_MASK) == DOUBLE_EXP_MASK) && ((conv.u & DOUBLE_FRAC_PART_MASK) != 0); -}; - -inline bool fastmath_isinf(double x) { - union { - uint64_t u; - double f; - } conv; - conv.f = x; - return ((conv.u & DOUBLE_EXP_MASK) == DOUBLE_EXP_MASK) && ((conv.u & DOUBLE_FRAC_PART_MASK) == 0); -}; - -} // namespace SourceXtractor - -#endif diff --git a/SEUtils/tests/src/NonFiniteMath_test.cpp b/SEUtils/tests/src/NonFiniteMath_test.cpp new file mode 100644 index 000000000..8f9b0091d --- /dev/null +++ b/SEUtils/tests/src/NonFiniteMath_test.cpp @@ -0,0 +1,98 @@ +/** Copyright © 2019 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 + */ +/** + * @file tests/src/NonFiniteMath_test.cpp + * + * Detection images legitimately contain NaN: for floating-point FITS the standard + * has no other way to mark an undefined pixel, since BLANK is integer-only, and + * mosaics blank their unexposed area that way. Under -ffinite-math-only a NaN is + * undefined behaviour rather than merely unlikely, so every isnan/isinf test is + * compiled out and the result of any comparison against one becomes arbitrary -- + * differing between compilers and between architectures. + * + * These checks fail loudly at build-verification time rather than leaving it to be + * discovered as a crash or a hang on real data. + */ + +#include + +#include +#include +#include +#include + +namespace { + +/// The bit pattern of a 32-bit quiet NaN, read through volatile so that no +/// amount of constant folding can decide the answers below at compile time +volatile uint32_t g_quiet_nan_bits = 0x7FC00000u; + +float runtimeNaN() { + uint32_t bits = g_quiet_nan_bits; + float value; + std::memcpy(&value, &bits, sizeof(value)); + return value; +} + +} // anonymous namespace + +BOOST_AUTO_TEST_SUITE(NonFiniteMath_test) + +//----------------------------------------------------------------------------- + +/// Both GCC and Clang define this to 1 exactly when -ffinite-math-only is in effect, +/// so unlike the behavioural checks below it cannot accidentally pass +BOOST_AUTO_TEST_CASE(finite_math_only_is_not_enabled_test) { +#ifdef __FINITE_MATH_ONLY__ + BOOST_CHECK_MESSAGE(__FINITE_MATH_ONLY__ == 0, + "built with -ffinite-math-only (implied by -ffast-math): NaN/Inf behavior " + "is ill defined. Configure with -DFINITE_MATH_ONLY=OFF, " + "see cmake/LocalBuildFlags.cmake"); +#endif +} + +//----------------------------------------------------------------------------- + +BOOST_AUTO_TEST_CASE(nan_is_classified_test) { + float nan_value = runtimeNaN(); + + BOOST_CHECK(std::isnan(nan_value)); + BOOST_CHECK(!std::isfinite(nan_value)); + // The canonical idiom, and the one an optimiser is most likely to fold away + BOOST_CHECK(nan_value != nan_value); +} + +//----------------------------------------------------------------------------- + +/// Two comparisons this code base actually relies on: the sentinel test in +/// ImageMode::processCell and the above-threshold test in Lutz::labelImage +BOOST_AUTO_TEST_CASE(nan_comparisons_are_ieee_test) { + float nan_value = runtimeNaN(); + const float sentinel = std::numeric_limits::lowest(); + + BOOST_CHECK(nan_value != sentinel); + BOOST_CHECK(!(nan_value == sentinel)); + + BOOST_CHECK(!(nan_value > 0.0f)); + BOOST_CHECK(!(nan_value < 0.0f)); + BOOST_CHECK(!(nan_value >= 0.0f)); + BOOST_CHECK(!(nan_value <= 0.0f)); +} + +//----------------------------------------------------------------------------- + +BOOST_AUTO_TEST_SUITE_END() diff --git a/cmake/LocalBuildFlags.cmake b/cmake/LocalBuildFlags.cmake index 6983d9854..714afdcf8 100644 --- a/cmake/LocalBuildFlags.cmake +++ b/cmake/LocalBuildFlags.cmake @@ -1,17 +1,30 @@ include_guard(GLOBAL) +# -ffast-math implies -ffinite-math-only, which makes NaN/Inf operations undefined behavior: +# Every isnan/isinf test is then compiled out, and the result of any comparison involving a NaN +# becomes unpredicatable (compiler/arch-dependent). Set this to OFF to keep the rest of +# -ffast-math while restoring IEEE semantics for non-finite values. +option(FINITE_MATH_ONLY "Let -ffast-math assume that no NaN or Inf ever occurs" OFF) + +set(FAST_MATH_FLAGS "-ffast-math") +if(NOT FINITE_MATH_ONLY) + # Must follow -ffast-math, which would otherwise re-enable the assumption + set(FAST_MATH_FLAGS "${FAST_MATH_FLAGS} -fno-finite-math-only") +endif() + + # Build type compilation flags (if different from default or unknown to CMake) -set(CMAKE_CXX_FLAGS_RELEASE "-O3 -ffast-math" +set(CMAKE_CXX_FLAGS_RELEASE "-O3 ${FAST_MATH_FLAGS}" CACHE STRING "Flags used by the compiler during release builds." FORCE) -set(CMAKE_C_FLAGS_RELEASE "-O3 -ffast-math" +set(CMAKE_C_FLAGS_RELEASE "-O3 ${FAST_MATH_FLAGS}" CACHE STRING "Flags used by the compiler during release builds." FORCE) - -set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "-O3 -ffast-math -g${DEBUG_FORMAT}${DEBUG_LEVEL}" + +set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "-O3 ${FAST_MATH_FLAGS} -g${DEBUG_FORMAT}${DEBUG_LEVEL}" CACHE STRING "Flags used by the compiler during Release with Debug Info builds." FORCE) -set(CMAKE_C_FLAGS_RELWITHDEBINFO "-O3 -ffast-math -g${DEBUG_FORMAT}${DEBUG_LEVEL}" +set(CMAKE_C_FLAGS_RELWITHDEBINFO "-O3 ${FAST_MATH_FLAGS} -g${DEBUG_FORMAT}${DEBUG_LEVEL}" CACHE STRING "Flags used by the compiler during Release with Debug Info builds." FORCE)