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/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..c34498445 100644 --- a/SEImplementation/SEImplementation/Background/SE/KappaSigmaBinning.h +++ b/SEImplementation/SEImplementation/Background/SE/KappaSigmaBinning.h @@ -19,6 +19,9 @@ #define SOURCEXTRACTORPLUSPLUS_KAPPASIGMABINNING_H #include +#include + +#include "ElementsKernel/Exception.h" #include "Histogram/Histogram.h" namespace SourceXtractor { @@ -67,6 +70,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 +103,15 @@ class KappaSigmaBinning: public Euclid::Histogram::BinStrategy { } 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 { /** * 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; } }; diff --git a/SEImplementation/src/lib/Background/SE/ImageMode.cpp b/SEImplementation/src/lib/Background/SE/ImageMode.cpp index 65ed307f8..095e08bad 100644 --- a/SEImplementation/src/lib/Background/SE/ImageMode.cpp +++ b/SEImplementation/src/lib/Background/SE/ImageMode.cpp @@ -15,8 +15,12 @@ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ +#include + #include // 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::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; @@ -96,7 +101,10 @@ 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 + 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::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); - 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); } } 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() //----------------------------------------------------------------------------- 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)