From 5dd2da6debe24458520abea8a79a8c1a0436fc8e Mon Sep 17 00:00:00 2001 From: Marc Schefer <31961290+marcschefer@users.noreply.github.com> Date: Thu, 22 May 2025 17:17:14 +0200 Subject: [PATCH 01/40] Update CMakeLists.txt Bump develop version --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index e658ef8f..265b89ed 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -6,4 +6,4 @@ find_package(ElementsProject) #--------------------------------------------------------------- # Declare project name and version -elements_project(SourceXtractorPlusPlus 1.0.3 USE Alexandria 2.32.0 DESCRIPTION "SourceXtractor++, the next generation SExtractor") +elements_project(SourceXtractorPlusPlus 1.1.0 USE Alexandria 2.32.0 DESCRIPTION "SourceXtractor++, the next generation SExtractor") From 519ad8f8d6f4553d7a4e06e1dcc5cd055ddd8e1a Mon Sep 17 00:00:00 2001 From: Marc Schefer <31961290+marcschefer@users.noreply.github.com> Date: Thu, 22 May 2025 17:35:18 +0200 Subject: [PATCH 02/40] Update README.md --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index b04dbdd0..3caad664 100644 --- a/README.md +++ b/README.md @@ -17,13 +17,13 @@ Documentation: https://astrorama.github.io/SourceXtractorPlusPlus/ SourceXtractor++ is available on [Anaconda Cloud for Linux and MacOSX](https://anaconda.org/astrorama/sourcextractor) ```bash -conda install -c conda-forge -c astrorama sourcextractor==1.0.2 +conda install -c conda-forge -c astrorama sourcextractor==1.0.3 ``` We would recommend, however, to install it into its own environment. ```bash -conda create -n sourcex -c astrorama -c conda-forge sourcextractor==1.0.2 +conda create -n sourcex -c astrorama -c conda-forge sourcextractor==1.0.3 conda activate sourcex ``` From 9bf1525ab767bba6ea465c889eddf16b3fd70cf9 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Fri, 27 Jun 2025 17:24:43 +0200 Subject: [PATCH 03/40] Optimize VariablePsfStack using KdTree --- SEBenchmarks/CMakeLists.txt | 2 + .../src/program/BenchVariablePsfStack.cpp | 156 +++++++++++ .../SEFramework/Psf/VariablePsfStack.h | 39 ++- SEFramework/src/lib/Psf/VariablePsfStack.cpp | 88 +++--- SEUtils/CMakeLists.txt | 3 + SEUtils/SEUtils/KdTree.h | 2 + SEUtils/SEUtils/_impl/KdTree.icpp | 48 ++++ SEUtils/tests/src/KdTree_test.cpp | 258 ++++++++++++++++++ 8 files changed, 553 insertions(+), 43 deletions(-) create mode 100644 SEBenchmarks/src/program/BenchVariablePsfStack.cpp create mode 100644 SEUtils/tests/src/KdTree_test.cpp diff --git a/SEBenchmarks/CMakeLists.txt b/SEBenchmarks/CMakeLists.txt index 8c45f204..0b5b7dd3 100644 --- a/SEBenchmarks/CMakeLists.txt +++ b/SEBenchmarks/CMakeLists.txt @@ -56,6 +56,8 @@ elements_add_executable(BenchRendering src/program/BenchRendering.cpp LINK_LIBRARIES SEFramework SEImplementation ${Boost_LIBRARIES}) elements_add_executable(BenchBackgroundModel src/program/BenchBackgroundModel.cpp LINK_LIBRARIES SEFramework SEImplementation ${Boost_LIBRARIES}) +elements_add_executable(BenchVariablePsfStack src/program/BenchVariablePsfStack.cpp + LINK_LIBRARIES SEFramework ${Boost_LIBRARIES}) #=============================================================================== # Declare the Boost tests here diff --git a/SEBenchmarks/src/program/BenchVariablePsfStack.cpp b/SEBenchmarks/src/program/BenchVariablePsfStack.cpp new file mode 100644 index 00000000..2e0c823e --- /dev/null +++ b/SEBenchmarks/src/program/BenchVariablePsfStack.cpp @@ -0,0 +1,156 @@ +/** Copyright © 2019-2025 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 src/program/BenchVariablePsfStack.cpp + * @date 06/27/25 + * @author marc schefer + */ + +#include +#include +#include + +#include +#include +#include "ElementsKernel/ProgramHeaders.h" +#include "ElementsKernel/Real.h" +#include "SEFramework/Psf/VariablePsfStack.h" + +namespace po = boost::program_options; +namespace timer = boost::timer; +using namespace SourceXtractor; + +static Elements::Logging logger = Elements::Logging::getLogger("BenchVariablePsfStack"); + +class BenchVariablePsfStack : public Elements::Program { +private: + std::default_random_engine random_generator; + std::uniform_real_distribution random_dist{0.0, 1000.0}; + +public: + + po::options_description defineSpecificProgramOptions() override { + po::options_description options{}; + options.add_options() + ("iterations", po::value()->default_value(100000), "Number of getPsf calls to benchmark") + ("measures", po::value()->default_value(10), "Number of timing measurements to take") + ("fits-file", po::value()->default_value(""), "FITS file containing PSF stack (optional, will use nullptr if not provided)"); + return options; + } + + Elements::ExitCode mainMethod(std::map &args) override { + + auto iterations = args["iterations"].as(); + auto measures = args["measures"].as(); + auto fits_file = args["fits-file"].as(); + + logger.info() << "Benchmarking VariablePsfStack::getPsf() with " << iterations << " iterations"; + logger.info() << "Taking " << measures << " timing measurements"; + + // Initialize VariablePsfStack with FITS file if provided, otherwise nullptr + std::shared_ptr fitsPtr = nullptr; + if (!fits_file.empty()) { + try { + fitsPtr = std::make_shared(fits_file); + logger.info() << "Using FITS file: " << fits_file; + } catch (const std::exception& e) { + logger.error() << "Failed to load FITS file '" << fits_file << "': " << e.what(); + logger.info() << "Continuing with nullptr - this will likely cause exceptions during getPsf() calls"; + } + } else { + logger.info() << "No FITS file provided - using nullptr (will likely cause exceptions)"; + } + + try { + VariablePsfStack psfStack(fitsPtr); + + logger.info() << "VariablePsfStack loaded successfully with " << psfStack.getNumberOfPsfs() << " PSFs"; + logger.info() << "PSF size: " << psfStack.getWidth() << "x" << psfStack.getHeight(); + logger.info() << "Pixel sampling: " << psfStack.getPixelSampling(); + + std::cout << "Iterations,Measurement,Time_nanoseconds" << std::endl; + + for (int m = 0; m < measures; ++m) { + logger.info() << "Measurement " << (m + 1) << "/" << measures; + + timer::cpu_timer timer; + timer.stop(); + + // Prepare test values for getPsf calls + std::vector> testValues; + testValues.reserve(iterations); + + for (int i = 0; i < iterations; ++i) { + testValues.push_back({random_dist(random_generator), random_dist(random_generator)}); + } + + // Start timing + timer.start(); + + for (int i = 0; i < iterations; ++i) { + try { + auto psf = psfStack.getPsf(testValues[i]); + // Prevent compiler optimization by using the result + volatile auto width = psf->getWidth(); + (void)width; // Suppress unused variable warning + } catch (const std::exception& e) { + // Expected to fail with nullptr, but we still measure the timing + // until the exception is thrown + } + } + + timer.stop(); + + auto elapsed_ns = timer.elapsed().wall; + std::cout << iterations << "," << (m + 1) << "," << elapsed_ns << std::endl; + + logger.info() << "Time for " << iterations << " calls: " << (elapsed_ns / 1e9) << " seconds"; + logger.info() << "Average time per call: " << (elapsed_ns / iterations) << " nanoseconds"; + } + + } catch (const std::exception& e) { + logger.error() << "Error initializing VariablePsfStack: " << e.what(); + logger.info() << "This is expected with nullptr parameter - fix the FITS file parameter later"; + + // Still run a basic timing test to measure overhead + std::cout << "Running basic timing test without actual PSF operations..." << std::endl; + std::cout << "Iterations,Measurement,Time_nanoseconds" << std::endl; + + for (int m = 0; m < measures; ++m) { + timer::cpu_timer timer; + timer.stop(); + + timer.start(); + for (int i = 0; i < iterations; ++i) { + // Just measure the overhead of generating random values and vector creation + std::vector values = {random_dist(random_generator), random_dist(random_generator)}; + volatile auto size = values.size(); + (void)size; + } + timer.stop(); + + auto elapsed_ns = timer.elapsed().wall; + std::cout << iterations << "," << (m + 1) << "," << elapsed_ns << std::endl; + } + } + + return Elements::ExitCode::OK; + } +}; + +MAIN_FOR(BenchVariablePsfStack) diff --git a/SEFramework/SEFramework/Psf/VariablePsfStack.h b/SEFramework/SEFramework/Psf/VariablePsfStack.h index 836b9cf1..c71fa045 100644 --- a/SEFramework/SEFramework/Psf/VariablePsfStack.h +++ b/SEFramework/SEFramework/Psf/VariablePsfStack.h @@ -25,8 +25,10 @@ #define _SEIMPLEMENTATION_PSF_VARIABLEPSFSTACK_H_ #include +#include #include #include +#include namespace SourceXtractor { @@ -42,6 +44,18 @@ namespace SourceXtractor { */ class VariablePsfStack final : public Psf { public: + /** + * @brief Structure to hold PSF position data + */ + struct PsfPosition { + SeFloat ra; + SeFloat dec; + SeFloat x; + SeFloat y; + double gridx; + double gridy; + }; + /** * Constructor */ @@ -83,6 +97,13 @@ class VariablePsfStack final : public Psf { return m_components; }; + /** + * @return The number of PSFs loaded in the stack + */ + long getNumberOfPsfs() const { + return m_nrows; + }; + /** * */ @@ -99,12 +120,8 @@ class VariablePsfStack final : public Psf { long m_nrows; - std::vector m_ra_values; - std::vector m_dec_values; - std::vector m_x_values; - std::vector m_y_values; - std::vector m_gridx_values; - std::vector m_gridy_values; + std::vector m_positions; + std::unique_ptr> m_kdtree; std::vector m_components = {"X_IMAGE", "Y_IMAGE"}; @@ -119,6 +136,16 @@ class VariablePsfStack final : public Psf { void selfTest(); }; +/** + * @brief KdTree traits specialization for PsfPosition + */ +template <> +struct KdTreeTraits { + static double getCoord(const VariablePsfStack::PsfPosition& pos, size_t index) { + return (index == 0) ? static_cast(pos.x) : static_cast(pos.y); + } +}; + } // namespace SourceXtractor #endif //_SEIMPLEMENTATION_PSF_VARIABLEPSFSTACK_H_ diff --git a/SEFramework/src/lib/Psf/VariablePsfStack.cpp b/SEFramework/src/lib/Psf/VariablePsfStack.cpp index 2276e55d..9b65da75 100644 --- a/SEFramework/src/lib/Psf/VariablePsfStack.cpp +++ b/SEFramework/src/lib/Psf/VariablePsfStack.cpp @@ -21,6 +21,7 @@ * Author: Martin Kuemmel */ #include +#include #include #include #include "SEFramework/Psf/VariablePsfStack.h" @@ -73,18 +74,33 @@ void VariablePsfStack::setup(std::shared_ptr pFits) { // read the nrows value m_nrows = position_data.rows(); + // Temporary vectors for reading data + std::vector ra_values, dec_values, x_values, y_values; + std::vector gridx_values, gridy_values; + try { // read in all the EXT specific columns - position_data.column("GRIDX", false).read(m_gridx_values, 0, m_nrows); - position_data.column("GRIDY", false).read(m_gridy_values, 0, m_nrows); + position_data.column("GRIDX", false).read(gridx_values, 0, m_nrows); + position_data.column("GRIDY", false).read(gridy_values, 0, m_nrows); } catch (CCfits::Table::NoSuchColumn) { - position_data.column("X_CENTER", false).read(m_gridx_values, 0, m_nrows); - position_data.column("Y_CENTER", false).read(m_gridy_values, 0, m_nrows); + position_data.column("X_CENTER", false).read(gridx_values, 0, m_nrows); + position_data.column("Y_CENTER", false).read(gridy_values, 0, m_nrows); + } + position_data.column("RA", false).read(ra_values, 0, m_nrows); + position_data.column("DEC", false).read(dec_values, 0, m_nrows); + position_data.column("X", false).read(x_values, 0, m_nrows); + position_data.column("Y", false).read(y_values, 0, m_nrows); + + // Populate the positions vector + m_positions.reserve(m_nrows); + for (long i = 0; i < m_nrows; ++i) { + m_positions.push_back({ra_values[i], dec_values[i], x_values[i], y_values[i], gridx_values[i], gridy_values[i]}); + } + + // Build KdTree for fast nearest neighbor searches + if (!m_positions.empty()) { + m_kdtree = std::make_unique>(m_positions); } - position_data.column("RA", false).read(m_ra_values, 0, m_nrows); - position_data.column("DEC", false).read(m_dec_values, 0, m_nrows); - position_data.column("X", false).read(m_x_values, 0, m_nrows); - position_data.column("Y", false).read(m_y_values, 0, m_nrows); } catch (CCfits::FitsException& e) { throw Elements::Exception() << "Error loading stacked PSF file: " << e.message(); @@ -95,54 +111,54 @@ void VariablePsfStack::selfTest() { int naxis1, naxis2; // read in the min/max grid values in x/y - const auto x_grid_minmax = std::minmax_element(begin(m_gridx_values), end(m_gridx_values)); - const auto y_grid_minmax = std::minmax_element(begin(m_gridy_values), end(m_gridy_values)); + auto x_grid_minmax = std::minmax_element(m_positions.begin(), m_positions.end(), + [](const PsfPosition& a, const PsfPosition& b) { return a.gridx < b.gridx; }); + auto y_grid_minmax = std::minmax_element(m_positions.begin(), m_positions.end(), + [](const PsfPosition& a, const PsfPosition& b) { return a.gridy < b.gridy; }); // read the image size m_pFits->extension(1).readKey("NAXIS1", naxis1); m_pFits->extension(1).readKey("NAXIS2", naxis2); // make sure all PSF in the grid are there - if (*x_grid_minmax.first - m_grid_offset < 1) - throw Elements::Exception() << "The PSF at the smallest x-grid starts at: " << *x_grid_minmax.first - m_grid_offset; - if (*y_grid_minmax.first - m_grid_offset < 1) - throw Elements::Exception() << "The PSF at the smallest y-grid starts at: " << *y_grid_minmax.first - m_grid_offset; - if (*x_grid_minmax.second + m_grid_offset > naxis1) - throw Elements::Exception() << "The PSF at the largest x-grid is too large: " << *x_grid_minmax.second + m_grid_offset + if (x_grid_minmax.first->gridx - m_grid_offset < 1) + throw Elements::Exception() << "The PSF at the smallest x-grid starts at: " << x_grid_minmax.first->gridx - m_grid_offset; + if (y_grid_minmax.first->gridy - m_grid_offset < 1) + throw Elements::Exception() << "The PSF at the smallest y-grid starts at: " << y_grid_minmax.first->gridy - m_grid_offset; + if (x_grid_minmax.second->gridx + m_grid_offset > naxis1) + throw Elements::Exception() << "The PSF at the largest x-grid is too large: " << x_grid_minmax.second->gridx + m_grid_offset << " NAXIS1: " << naxis1; - if (*y_grid_minmax.second + m_grid_offset > naxis2) - throw Elements::Exception() << "The PSF at the largest y-grid is too large: " << *y_grid_minmax.second + m_grid_offset - << " NAXIS2: " << naxis1; + if (y_grid_minmax.second->gridy + m_grid_offset > naxis2) + throw Elements::Exception() << "The PSF at the largest y-grid is too large: " << y_grid_minmax.second->gridy + m_grid_offset + << " NAXIS2: " << naxis2; } std::shared_ptr> VariablePsfStack::getPsf(const std::vector& values) const { - long index_min_distance = 0; - double min_distance = 1.0e+32; - // make sure there are only two positions if (values.size() != 2) throw Elements::Exception() << "There can be only two positional value for the stacked PSF!"; - // find the position of minimal distance - for (int act_index = 0; act_index < m_nrows; act_index++) { - double act_distance = (values[0] - m_x_values[act_index]) * (values[0] - m_x_values[act_index]) + - (values[1] - m_y_values[act_index]) * (values[1] - m_y_values[act_index]); - if (act_distance < min_distance) { - index_min_distance = act_index; - min_distance = act_distance; - } - } + // Use KdTree to find the nearest PSF position + KdTree::Coord coord; + coord.coord[0] = values[0]; // x coordinate + coord.coord[1] = values[1]; // y coordinate + + PsfPosition nearest_position = m_kdtree->findNearest(coord); + + // Calculate distance for logging + double min_distance = (values[0] - nearest_position.x) * (values[0] - nearest_position.x) + + (values[1] - nearest_position.y) * (values[1] - nearest_position.y); + // give some feedback stack_logger.debug() << "Distance: " << sqrt(min_distance) << " (" << values[0] << "," << values[1] << ")<-->(" - << m_x_values[index_min_distance] << "," << m_y_values[index_min_distance] - << ") index: " << index_min_distance; + << nearest_position.x << "," << nearest_position.y << ")"; // get the first and last pixels for the PSF to be extracted // NOTE: CCfits has 1-based indices, also the last index is *included* in the reading // NOTE: the +0.5 forces a correct cast/ceiling - std::vector first_vertex{long(m_gridx_values[index_min_distance]+.5) - long(m_grid_offset), long(m_gridy_values[index_min_distance]+.5) - long(m_grid_offset)}; + std::vector first_vertex{long(nearest_position.gridx+.5) - long(m_grid_offset), long(nearest_position.gridy+.5) - long(m_grid_offset)}; stack_logger.debug() << "First vertex: ( " << first_vertex[0] << ", " << first_vertex[1] << ") First vertex alternative: " << - m_gridx_values[index_min_distance]-m_grid_offset << " " << m_gridy_values[index_min_distance]-m_grid_offset << + nearest_position.gridx-m_grid_offset << " " << nearest_position.gridy-m_grid_offset << " grid offset:" << m_grid_offset; std::vector last_vertex{first_vertex[0] + long(m_psf_size) - 1, first_vertex[1] +long( m_psf_size) - 1}; @@ -155,8 +171,6 @@ std::shared_ptr> VariablePsfStack::getPsf(const std::vector m_pFits->extension(1).read(stamp_data, first_vertex, last_vertex, stride); } - //stack_logger.info() << "DDD ( " << first_vertex[0] << ", " << first_vertex[1] << ") --> ( " << last_vertex[0] << ", " << last_vertex[1] << "): " << stamp_data.size(); - // create and return the psf image return VectorImage::create(m_psf_size, m_psf_size, std::begin(stamp_data), std::end(stamp_data)); } diff --git a/SEUtils/CMakeLists.txt b/SEUtils/CMakeLists.txt index 9f3c453a..64bcc685 100644 --- a/SEUtils/CMakeLists.txt +++ b/SEUtils/CMakeLists.txt @@ -66,6 +66,9 @@ elements_add_unit_test(Misc_test tests/src/Misc_test.cpp elements_add_unit_test(QuadTree_test tests/src/QuadTree_test.cpp LINK_LIBRARIES SEUtils TYPE Boost) +elements_add_unit_test(KdTree_test tests/src/KdTree_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/KdTree.h b/SEUtils/SEUtils/KdTree.h index 5c4fc354..3e739c2e 100644 --- a/SEUtils/SEUtils/KdTree.h +++ b/SEUtils/SEUtils/KdTree.h @@ -21,6 +21,7 @@ #include #include #include +#include namespace SourceXtractor { @@ -49,6 +50,7 @@ class KdTree { explicit KdTree(const std::vector& data); std::vector findPointsWithinRadius(Coord coord, double radius) const; + T findNearest(Coord coord) const; private: class Node; diff --git a/SEUtils/SEUtils/_impl/KdTree.icpp b/SEUtils/SEUtils/_impl/KdTree.icpp index 606b939c..ce45e4d2 100644 --- a/SEUtils/SEUtils/_impl/KdTree.icpp +++ b/SEUtils/SEUtils/_impl/KdTree.icpp @@ -21,6 +21,7 @@ template class KdTree::Node { public: virtual std::vector findPointsWithinRadius(Coord coord, double radius) const = 0; + virtual T findNearest(Coord coord, double& best_dist_sq) const = 0; virtual ~Node() = default; }; @@ -45,6 +46,24 @@ public: return selection; } + virtual T findNearest(Coord coord, double& best_dist_sq) const { + T best_point = m_data[0]; + best_dist_sq = std::numeric_limits::max(); + + for (const auto& entry : m_data) { + double square_dist = 0.0; + for (size_t i = 0; i < N; i++) { + double delta = Traits::getCoord(entry, i) - coord.coord[i]; + square_dist += delta * delta; + } + if (square_dist < best_dist_sq) { + best_dist_sq = square_dist; + best_point = entry; + } + } + return best_point; + } + private: const std::vector m_data; }; @@ -101,6 +120,29 @@ public: } } + virtual T findNearest(Coord coord, double& best_dist_sq) const { + // Determine which side to search first + bool search_left_first = coord.coord[m_axis] < m_split_value; + auto first_child = search_left_first ? m_left_child : m_right_child; + auto second_child = search_left_first ? m_right_child : m_left_child; + + // Search the first side + T best_point = first_child->findNearest(coord, best_dist_sq); + + // Check if we need to search the other side + double axis_dist = coord.coord[m_axis] - m_split_value; + if (axis_dist * axis_dist < best_dist_sq) { + double second_best_dist_sq = best_dist_sq; + T second_best = second_child->findNearest(coord, second_best_dist_sq); + if (second_best_dist_sq < best_dist_sq) { + best_dist_sq = second_best_dist_sq; + best_point = second_best; + } + } + + return best_point; + } + private: size_t m_axis; double m_split_value; @@ -124,4 +166,10 @@ std::vector KdTree::findPointsWithinRadius(Coord coord, double radiu return m_root->findPointsWithinRadius(coord, radius); } +template +T KdTree::findNearest(Coord coord) const { + double best_dist_sq; + return m_root->findNearest(coord, best_dist_sq); +} + } diff --git a/SEUtils/tests/src/KdTree_test.cpp b/SEUtils/tests/src/KdTree_test.cpp new file mode 100644 index 00000000..5c5cf9e1 --- /dev/null +++ b/SEUtils/tests/src/KdTree_test.cpp @@ -0,0 +1,258 @@ +/** Copyright © 2021 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 "SEUtils/KdTree.h" + +namespace SourceXtractor { + +struct TestPoint { + double m_x, m_y; + int m_id; // For identification in tests + + TestPoint(double x, double y, int id = 0) : m_x(x), m_y(y), m_id(id) {} +}; + +bool operator==(const TestPoint& a, const TestPoint& b) { + return a.m_x == b.m_x && a.m_y == b.m_y && a.m_id == b.m_id; +} + +template <> +struct KdTreeTraits { + static double getCoord(const TestPoint& t, size_t index) { + if (index == 0) { + return t.m_x; + } else { + return t.m_y; + } + } +}; + +// Test point for 3D tests +struct TestPoint3D { + double m_x, m_y, m_z; + int m_id; + + TestPoint3D(double x, double y, double z, int id = 0) : m_x(x), m_y(y), m_z(z), m_id(id) {} +}; + +bool operator==(const TestPoint3D& a, const TestPoint3D& b) { + return a.m_x == b.m_x && a.m_y == b.m_y && a.m_z == b.m_z && a.m_id == b.m_id; +} + +template <> +struct KdTreeTraits { + static double getCoord(const TestPoint3D& t, size_t index) { + if (index == 0) { + return t.m_x; + } else if (index == 1) { + return t.m_y; + } else { + return t.m_z; + } + } +}; + +BOOST_AUTO_TEST_SUITE (KdTree_test) + +BOOST_AUTO_TEST_CASE( smoke_test ) { + std::vector data = {{1.0, 2.0, 1}, {3.0, 4.0, 2}, {5.0, 6.0, 3}}; + KdTree tree(data); +} + +BOOST_AUTO_TEST_CASE( empty_tree_test ) { + std::vector data = {}; + // This should not crash - testing edge case + BOOST_CHECK_NO_THROW(KdTree tree(data)); +} + +BOOST_AUTO_TEST_CASE( single_point_test ) { + std::vector data = {{1.0, 2.0, 1}}; + KdTree tree(data); + + KdTree::Coord coord; + coord.coord[0] = 1.0; + coord.coord[1] = 2.0; + + auto result = tree.findPointsWithinRadius(coord, 0.1); + BOOST_CHECK_EQUAL(result.size(), 1); + BOOST_CHECK(result[0] == TestPoint(1.0, 2.0, 1)); + + auto nearest = tree.findNearest(coord); + BOOST_CHECK(nearest == TestPoint(1.0, 2.0, 1)); +} + +BOOST_AUTO_TEST_CASE( simple_radius_search_test ) { + std::vector data = { + {1.0, 1.0, 1}, {2.0, 2.0, 2}, {3.0, 3.0, 3}, {10.0, 10.0, 4} + }; + KdTree tree(data); + + KdTree::Coord coord; + coord.coord[0] = 2.0; + coord.coord[1] = 2.0; + + // Search with radius 2.0 - should find points at (1,1), (2,2), (3,3) + auto result = tree.findPointsWithinRadius(coord, 2.0); + BOOST_CHECK_EQUAL(result.size(), 3); + + // Search with radius 1.0 - should find only point at (2,2) + result = tree.findPointsWithinRadius(coord, 1.0); + BOOST_CHECK_EQUAL(result.size(), 1); + BOOST_CHECK(result[0] == TestPoint(2.0, 2.0, 2)); + + // Search with large radius - should find all points + result = tree.findPointsWithinRadius(coord, 15.0); + BOOST_CHECK_EQUAL(result.size(), 4); +} + +BOOST_AUTO_TEST_CASE( nearest_neighbor_test ) { + std::vector data = { + {1.0, 1.0, 1}, {2.0, 2.0, 2}, {3.0, 3.0, 3}, {10.0, 10.0, 4} + }; + KdTree tree(data); + + KdTree::Coord coord; + coord.coord[0] = 2.1; + coord.coord[1] = 1.9; + + auto nearest = tree.findNearest(coord); + BOOST_CHECK(nearest == TestPoint(2.0, 2.0, 2)); + + // Test corner case - equidistant points + coord.coord[0] = 1.5; + coord.coord[1] = 1.5; + nearest = tree.findNearest(coord); + // Should return one of the two closest points (1,1) or (2,2) + BOOST_CHECK(nearest == TestPoint(1.0, 1.0, 1) || nearest == TestPoint(2.0, 2.0, 2)); +} + +BOOST_AUTO_TEST_CASE( large_dataset_test ) { + std::vector data; + + // Create a grid of 1000 points + for (int i = 0; i < 1000; i++) { + double x = 5.0 + (i % 100) / 10.0; + double y = 2.0 + (i / 100) / 5.0; + data.emplace_back(x, y, i); + } + + KdTree tree(data); + + KdTree::Coord coord; + coord.coord[0] = 7.0; + coord.coord[1] = 3.0; + + // Test radius search on large dataset + auto result = tree.findPointsWithinRadius(coord, 1.0); + BOOST_CHECK_GT(result.size(), 0); + + // Test nearest neighbor on large dataset + auto nearest = tree.findNearest(coord); + BOOST_CHECK_GT(nearest.m_id, -1); // Should find a valid point +} + +BOOST_AUTO_TEST_CASE( three_dimensional_test ) { + std::vector data = { + {1.0, 1.0, 1.0, 1}, {2.0, 2.0, 2.0, 2}, {3.0, 3.0, 3.0, 3}, {10.0, 10.0, 10.0, 4} + }; + KdTree tree(data); + + KdTree::Coord coord; + coord.coord[0] = 2.0; + coord.coord[1] = 2.0; + coord.coord[2] = 2.0; + + // Test radius search in 3D + auto result = tree.findPointsWithinRadius(coord, 2.0); + BOOST_CHECK_EQUAL(result.size(), 3); // Should find first 3 points + + // Test nearest neighbor in 3D + coord.coord[0] = 2.1; + coord.coord[1] = 1.9; + coord.coord[2] = 2.1; + + auto nearest = tree.findNearest(coord); + BOOST_CHECK(nearest == TestPoint3D(2.0, 2.0, 2.0, 2)); +} + +BOOST_AUTO_TEST_CASE( small_leaf_size_test ) { + std::vector data = { + {1.0, 1.0, 1}, {2.0, 2.0, 2}, {3.0, 3.0, 3}, {4.0, 4.0, 4}, {5.0, 5.0, 5} + }; + + // Use small leaf size to force tree structure + KdTree tree(data); + + KdTree::Coord coord; + coord.coord[0] = 3.0; + coord.coord[1] = 3.0; + + auto result = tree.findPointsWithinRadius(coord, 1.5); + BOOST_CHECK_GE(result.size(), 1); + + auto nearest = tree.findNearest(coord); + BOOST_CHECK(nearest == TestPoint(3.0, 3.0, 3)); +} + +BOOST_AUTO_TEST_CASE( edge_case_identical_coordinates ) { + std::vector data = { + {1.0, 1.0, 1}, {1.0, 1.0, 2}, {1.0, 1.0, 3}, {2.0, 2.0, 4} + }; + KdTree tree(data); + + KdTree::Coord coord; + coord.coord[0] = 1.0; + coord.coord[1] = 1.0; + + auto result = tree.findPointsWithinRadius(coord, 0.1); + BOOST_CHECK_EQUAL(result.size(), 3); // Should find all three identical points + + auto nearest = tree.findNearest(coord); + // Should return one of the three identical points + BOOST_CHECK(nearest.m_x == 1.0 && nearest.m_y == 1.0); +} + +BOOST_AUTO_TEST_CASE( stress_test_performance ) { + std::vector data; + + // Create 10000 random-ish points + for (int i = 0; i < 10000; i++) { + double x = (i * 17) % 1000 / 10.0; // Pseudo-random distribution + double y = (i * 23) % 1000 / 10.0; + data.emplace_back(x, y, i); + } + + KdTree tree(data); + + KdTree::Coord coord; + coord.coord[0] = 50.0; + coord.coord[1] = 50.0; + + // This should complete reasonably quickly due to tree structure + auto result = tree.findPointsWithinRadius(coord, 10.0); + BOOST_CHECK_GE(result.size(), 0); + + auto nearest = tree.findNearest(coord); + BOOST_CHECK_GE(nearest.m_id, 0); +} + +BOOST_AUTO_TEST_SUITE_END () + +} From da6f69c6b7ddc99f5bfc15ad7038a74640f6a247 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Mon, 30 Jun 2025 11:42:47 +0200 Subject: [PATCH 04/40] benchmark code cleanup --- .../src/program/BenchVariablePsfStack.cpp | 32 ++++--------------- 1 file changed, 6 insertions(+), 26 deletions(-) diff --git a/SEBenchmarks/src/program/BenchVariablePsfStack.cpp b/SEBenchmarks/src/program/BenchVariablePsfStack.cpp index 2e0c823e..e3de2eb8 100644 --- a/SEBenchmarks/src/program/BenchVariablePsfStack.cpp +++ b/SEBenchmarks/src/program/BenchVariablePsfStack.cpp @@ -48,8 +48,8 @@ class BenchVariablePsfStack : public Elements::Program { po::options_description options{}; options.add_options() ("iterations", po::value()->default_value(100000), "Number of getPsf calls to benchmark") - ("measures", po::value()->default_value(10), "Number of timing measurements to take") - ("fits-file", po::value()->default_value(""), "FITS file containing PSF stack (optional, will use nullptr if not provided)"); + ("measures", po::value()->default_value(3), "Number of timing measurements to take") + ("fits-file", po::value()->default_value(""), "FITS file containing PSF stack"); return options; } @@ -70,10 +70,11 @@ class BenchVariablePsfStack : public Elements::Program { logger.info() << "Using FITS file: " << fits_file; } catch (const std::exception& e) { logger.error() << "Failed to load FITS file '" << fits_file << "': " << e.what(); - logger.info() << "Continuing with nullptr - this will likely cause exceptions during getPsf() calls"; + return Elements::ExitCode::DATAERR; } } else { - logger.info() << "No FITS file provided - using nullptr (will likely cause exceptions)"; + logger.error() << "No FITS file provided"; + return Elements::ExitCode::USAGE; } try { @@ -125,28 +126,7 @@ class BenchVariablePsfStack : public Elements::Program { } catch (const std::exception& e) { logger.error() << "Error initializing VariablePsfStack: " << e.what(); - logger.info() << "This is expected with nullptr parameter - fix the FITS file parameter later"; - - // Still run a basic timing test to measure overhead - std::cout << "Running basic timing test without actual PSF operations..." << std::endl; - std::cout << "Iterations,Measurement,Time_nanoseconds" << std::endl; - - for (int m = 0; m < measures; ++m) { - timer::cpu_timer timer; - timer.stop(); - - timer.start(); - for (int i = 0; i < iterations; ++i) { - // Just measure the overhead of generating random values and vector creation - std::vector values = {random_dist(random_generator), random_dist(random_generator)}; - volatile auto size = values.size(); - (void)size; - } - timer.stop(); - - auto elapsed_ns = timer.elapsed().wall; - std::cout << iterations << "," << (m + 1) << "," << elapsed_ns << std::endl; - } + return Elements::ExitCode::DATAERR; } return Elements::ExitCode::OK; From 8f1e8a4b1859d2388de86427cb82934a2f4b1783 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Mon, 30 Jun 2025 11:44:36 +0200 Subject: [PATCH 05/40] update .gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 8b9266d5..7b4781cc 100644 --- a/.gitignore +++ b/.gitignore @@ -11,3 +11,4 @@ __pycache__ /packages /doc/build/ /build/ +.vscode From 658644b44ed5bc405d2eb10db84f585524233185 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Fri, 8 Aug 2025 17:19:00 +0200 Subject: [PATCH 06/40] changed fitting window calculation to use doubles --- .../CoordinateSystem/CoordinateSystem.h | 24 ++++++++++ .../MeasurementFrameRectangle.h | 21 ++++++--- .../FlexibleModelFittingIterativeTask.cpp | 16 ++++--- .../MeasurementFrameRectangleTask.cpp | 21 ++++----- .../MeasurementFrameRectangleTaskNoDetect.cpp | 44 +++++++++---------- 5 files changed, 78 insertions(+), 48 deletions(-) diff --git a/SEFramework/SEFramework/CoordinateSystem/CoordinateSystem.h b/SEFramework/SEFramework/CoordinateSystem/CoordinateSystem.h index d43e26a4..30605da8 100644 --- a/SEFramework/SEFramework/CoordinateSystem/CoordinateSystem.h +++ b/SEFramework/SEFramework/CoordinateSystem/CoordinateSystem.h @@ -45,6 +45,30 @@ struct ImageCoordinate { ImageCoordinate() : m_x(0), m_y(0) {} ImageCoordinate(double x, double y) : m_x(x), m_y(y) {} + + ImageCoordinate operator*(double scalar) const { + return ImageCoordinate(m_x * scalar, m_y * scalar); + } + + ImageCoordinate operator+(const ImageCoordinate& other) const { + return ImageCoordinate(m_x + other.m_x, m_y + other.m_y); + } + + ImageCoordinate& operator+=(const ImageCoordinate& other) { + m_x += other.m_x; + m_y += other.m_y; + return *this; + } + + ImageCoordinate operator-(const ImageCoordinate& other) const { + return ImageCoordinate(m_x - other.m_x, m_y - other.m_y); + } + + ImageCoordinate& operator-=(const ImageCoordinate& other) { + m_x -= other.m_x; + m_y -= other.m_y; + return *this; + } }; class CoordinateSystem { diff --git a/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h b/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h index fb7ae567..db5e37b5 100644 --- a/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h +++ b/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h @@ -24,7 +24,10 @@ #ifndef _SEIMPLEMENTATION_PLUGIN_MEASUREMENTFRAMERECTANGLE_MEASUREMENTFRAMERECTANGLE_H_ #define _SEIMPLEMENTATION_PLUGIN_MEASUREMENTFRAMERECTANGLE_MEASUREMENTFRAMERECTANGLE_H_ +#include + #include "SEUtils/PixelRectangle.h" +#include "SEFramework/CoordinateSystem/CoordinateSystem.h" #include "SEFramework/Property/Property.h" #include "SEFramework/Image/Image.h" @@ -36,11 +39,14 @@ class MeasurementFrameRectangle: public Property { virtual ~MeasurementFrameRectangle() = default; explicit MeasurementFrameRectangle(bool bad_projection): - m_min_coord{-1, -1}, m_max_coord{-1, -1}, m_bad_projection{bad_projection}{} - - MeasurementFrameRectangle(PixelCoordinate min_coord, PixelCoordinate max_coord): - m_min_coord{min_coord}, m_max_coord{max_coord}, m_bad_projection{false} { - assert(min_coord.m_x <= max_coord.m_x && min_coord.m_y <= max_coord.m_y); + m_min_coord{-1, -1}, m_max_coord{-1, -1}, m_min_coord_image{-1, -1}, + m_max_coord_image{-1, -1}, m_bad_projection{bad_projection} {} + + MeasurementFrameRectangle(ImageCoordinate min_coord_image, ImageCoordinate max_coord_image): + m_min_coord_image(min_coord_image), m_max_coord_image(max_coord_image), m_bad_projection{false} { + assert(min_coord_image.m_x <= max_coord_image.m_x && min_coord_image.m_y <= max_coord_image.m_y); + m_min_coord = PixelCoordinate(static_cast(min_coord_image.m_x + 0.5), static_cast(min_coord_image.m_y + 0.5)); + m_max_coord = PixelCoordinate(static_cast(max_coord_image.m_x + 0.5), static_cast(max_coord_image.m_y + 0.5)); } PixelCoordinate getTopLeft() const { @@ -65,8 +71,8 @@ class MeasurementFrameRectangle: public Property { return m_max_coord.m_y - m_min_coord.m_y + 1; } - PixelRectangle getRect() const { - return PixelRectangle(m_min_coord, m_max_coord); + std::tuple getImageRect() const { + return std::make_tuple(m_min_coord_image, m_max_coord_image); } bool badProjection() const { @@ -75,6 +81,7 @@ class MeasurementFrameRectangle: public Property { private: PixelCoordinate m_min_coord, m_max_coord; + ImageCoordinate m_min_coord_image, m_max_coord_image; bool m_bad_projection; }; diff --git a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp index 1e99d9ea..0c1425d5 100644 --- a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp +++ b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp @@ -74,7 +74,8 @@ FlexibleModelFittingIterativeTask::~FlexibleModelFittingIterativeTask() { } PixelRectangle FlexibleModelFittingIterativeTask::getUnclippedFittingRect(SourceInterface& source, int frame_index) const { - auto fitting_rect = source.getProperty(frame_index).getRect(); + ImageCoordinate min_coord, max_coord; + std::tie(min_coord, max_coord) = source.getProperty(frame_index).getImageRect(); if (m_window_type == WindowType::ROTATED_ELLIPSE) { auto ellipse = getFittingEllipse(source, frame_index); @@ -84,14 +85,14 @@ PixelRectangle FlexibleModelFittingIterativeTask::getUnclippedFittingRect(Source return getEllipseRect(ellipse); } - if (fitting_rect.getWidth() <= 0 || fitting_rect.getHeight() <= 0) { + if (max_coord.m_x - min_coord.m_x <= 0 || max_coord.m_y - min_coord.m_y <= 0) { return PixelRectangle(); } else { - auto min = fitting_rect.getTopLeft(); - auto max = fitting_rect.getBottomRight(); + auto min = min_coord; + auto max = max_coord; // FIXME temporary, for now just enlarge the area by a fixed amount of pixels - PixelCoordinate border = (max - min) * .8 + PixelCoordinate(2, 2); + ImageCoordinate border = (max - min) * .8 + ImageCoordinate(2.0, 2.0); min -= border; max += border; @@ -125,7 +126,10 @@ PixelRectangle FlexibleModelFittingIterativeTask::getUnclippedFittingRect(Source max.m_y = min.m_y + size; } - return PixelRectangle(min, max); + auto min_pc = PixelCoordinate(static_cast(min.m_x + 0.5), static_cast(min.m_y + 0.5)); + auto max_pc = PixelCoordinate(static_cast(max.m_x + 0.5), static_cast(max.m_y + 0.5)); + + return PixelRectangle(min_pc, max_pc); } } diff --git a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp index f57d56fb..341ff2f2 100644 --- a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp +++ b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp @@ -67,25 +67,20 @@ void MeasurementFrameRectangleTask::computeProperties(SourceInterface& source) c auto max_x = std::max(coord1.m_x, std::max(coord2.m_x, std::max(coord3.m_x, coord4.m_x))); auto max_y = std::max(coord1.m_y, std::max(coord2.m_y, std::max(coord3.m_y, coord4.m_y))); - PixelCoordinate min_coord, max_coord; - min_coord.m_x = int(min_x); - min_coord.m_y = int(min_y); - max_coord.m_x = int(max_x) + 1; - max_coord.m_y = int(max_y) + 1; - // The full boundaries may lie outside of the frame - if (bad_coordinates || max_coord.m_x < 0 || max_coord.m_y < 0 || - min_coord.m_x >= measurement_frame_info.getWidth() || min_coord.m_y >= measurement_frame_info.getHeight()) { + if (bad_coordinates || max_x < 0.0 || max_y < 0.0 || + min_x >= measurement_frame_info.getWidth() - 0.5 || min_y >= measurement_frame_info.getHeight() - 0.5) { source.setIndexedProperty(m_instance, bad_coordinates); } // Clip the coordinates to fit the available image else { - min_coord.m_x = std::max(0, min_coord.m_x); - min_coord.m_y = std::max(0, min_coord.m_y); - max_coord.m_x = std::min(measurement_frame_info.getWidth() - 1, max_coord.m_x); - max_coord.m_y = std::min(measurement_frame_info.getHeight() - 1, max_coord.m_y); + min_x = std::max(0.0, min_x); + min_y = std::max(0.0, min_y); + max_x = std::min(double(measurement_frame_info.getWidth()) - 0.5, max_x); + max_y = std::min(double(measurement_frame_info.getHeight()) - 0.5, max_y); - source.setIndexedProperty(m_instance, min_coord, max_coord); + source.setIndexedProperty( + m_instance, ImageCoordinate(min_x, min_y), ImageCoordinate(max_x, max_y)); } } diff --git a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp index 68caedfc..c994d17b 100644 --- a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp +++ b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp @@ -23,6 +23,7 @@ #include "SEImplementation/Plugin/DetectionFrameCoordinates/DetectionFrameCoordinates.h" #include #include +#include #include #include @@ -31,23 +32,27 @@ namespace SourceXtractor { void MeasurementFrameRectangleTaskNoDetect::computeProperties(SourceInterface& source) const { auto measurement_frame_coordinates = source.getProperty(m_instance).getCoordinateSystem(); + auto reference_frame_coordinates = source.getProperty().getCoordinateSystem(); + const auto& measurement_frame_info = source.getProperty(m_instance); const auto& world_centroid = source.getProperty(); const auto& assoc_mode = source.getProperty(); - auto coord = world_centroid.getCentroid(); - bool bad_coordinates = false; ImageCoordinate coord1, coord2, coord3, coord4; try { - int w = assoc_mode.getRefFramePixelWidth(); - int h = assoc_mode.getRefFramePixelHeight(); + auto w = assoc_mode.getRefFramePixelWidth(); + auto h = assoc_mode.getRefFramePixelHeight(); - auto c = measurement_frame_coordinates->worldToImage(coord); - coord1 = ImageCoordinate(c.m_x - w, c.m_y - h); - coord2 = ImageCoordinate(c.m_x + w, c.m_y - h); - coord3 = ImageCoordinate(c.m_x - w, c.m_y + h); - coord4 = ImageCoordinate(c.m_x + w, c.m_y + h); + auto c = reference_frame_coordinates->worldToImage(world_centroid.getCentroid()); + coord1 = measurement_frame_coordinates->worldToImage( + reference_frame_coordinates->imageToWorld(ImageCoordinate(c.m_x - w, c.m_y - h))); + coord2 = measurement_frame_coordinates->worldToImage( + reference_frame_coordinates->imageToWorld(ImageCoordinate(c.m_x + w, c.m_y - h))); + coord3 = measurement_frame_coordinates->worldToImage( + reference_frame_coordinates->imageToWorld(ImageCoordinate(c.m_x - w, c.m_y + h))); + coord4 = measurement_frame_coordinates->worldToImage( + reference_frame_coordinates->imageToWorld(ImageCoordinate(c.m_x + w, c.m_y + h))); } catch (const InvalidCoordinatesException&) { bad_coordinates = true; @@ -59,25 +64,20 @@ void MeasurementFrameRectangleTaskNoDetect::computeProperties(SourceInterface& s auto max_x = std::max(coord1.m_x, std::max(coord2.m_x, std::max(coord3.m_x, coord4.m_x))); auto max_y = std::max(coord1.m_y, std::max(coord2.m_y, std::max(coord3.m_y, coord4.m_y))); - PixelCoordinate min_coord, max_coord; - min_coord.m_x = int(min_x); - min_coord.m_y = int(min_y); - max_coord.m_x = int(max_x) + 1; - max_coord.m_y = int(max_y) + 1; - // The full boundaries may lie outside of the frame - if (bad_coordinates || max_coord.m_x < 0 || max_coord.m_y < 0 || - min_coord.m_x >= measurement_frame_info.getWidth() || min_coord.m_y >= measurement_frame_info.getHeight()) { + if (bad_coordinates || max_x < 0.0 || max_y < 0.0 || + min_x >= measurement_frame_info.getWidth() - 0.5 || min_y >= measurement_frame_info.getHeight() - 0.5) { source.setIndexedProperty(m_instance, bad_coordinates); } // Clip the coordinates to fit the available image else { - min_coord.m_x = std::max(0, min_coord.m_x); - min_coord.m_y = std::max(0, min_coord.m_y); - max_coord.m_x = std::min(measurement_frame_info.getWidth() - 1, max_coord.m_x); - max_coord.m_y = std::min(measurement_frame_info.getHeight() - 1, max_coord.m_y); + min_x = std::max(0.0, min_x); + min_y = std::max(0.0, min_y); + max_x = std::min(double(measurement_frame_info.getWidth()) - 0.5, max_x); + max_y = std::min(double(measurement_frame_info.getHeight()) - 0.5, max_y); - source.setIndexedProperty(m_instance, min_coord, max_coord); + source.setIndexedProperty( + m_instance, ImageCoordinate(min_x, min_y), ImageCoordinate(max_x, max_y)); } } From 0bbbf7145c1a914a49abdbb39d49aa15cc5553ef Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Tue, 12 Aug 2025 13:59:47 +0200 Subject: [PATCH 07/40] fix regression --- .../MeasurementFrameRectangle.h | 7 ++++--- .../MeasurementFrameRectangleTask.cpp | 20 ++++++++++++------- .../MeasurementFrameRectangleTaskNoDetect.cpp | 20 ++++++++++++------- 3 files changed, 30 insertions(+), 17 deletions(-) diff --git a/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h b/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h index db5e37b5..d163e632 100644 --- a/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h +++ b/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h @@ -42,11 +42,12 @@ class MeasurementFrameRectangle: public Property { m_min_coord{-1, -1}, m_max_coord{-1, -1}, m_min_coord_image{-1, -1}, m_max_coord_image{-1, -1}, m_bad_projection{bad_projection} {} - MeasurementFrameRectangle(ImageCoordinate min_coord_image, ImageCoordinate max_coord_image): + MeasurementFrameRectangle(PixelCoordinate min_coord, PixelCoordinate max_coord, + ImageCoordinate min_coord_image, ImageCoordinate max_coord_image): + m_min_coord(min_coord), m_max_coord(max_coord), m_min_coord_image(min_coord_image), m_max_coord_image(max_coord_image), m_bad_projection{false} { + assert(min_coord.m_x <= max_coord.m_x && min_coord.m_y <= max_coord.m_y); assert(min_coord_image.m_x <= max_coord_image.m_x && min_coord_image.m_y <= max_coord_image.m_y); - m_min_coord = PixelCoordinate(static_cast(min_coord_image.m_x + 0.5), static_cast(min_coord_image.m_y + 0.5)); - m_max_coord = PixelCoordinate(static_cast(max_coord_image.m_x + 0.5), static_cast(max_coord_image.m_y + 0.5)); } PixelCoordinate getTopLeft() const { diff --git a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp index 341ff2f2..da88a8dc 100644 --- a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp +++ b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp @@ -67,20 +67,26 @@ void MeasurementFrameRectangleTask::computeProperties(SourceInterface& source) c auto max_x = std::max(coord1.m_x, std::max(coord2.m_x, std::max(coord3.m_x, coord4.m_x))); auto max_y = std::max(coord1.m_y, std::max(coord2.m_y, std::max(coord3.m_y, coord4.m_y))); + PixelCoordinate min_coord, max_coord; + min_coord.m_x = int(min_x); + min_coord.m_y = int(min_y); + max_coord.m_x = int(max_x) + 1; + max_coord.m_y = int(max_y) + 1; + // The full boundaries may lie outside of the frame - if (bad_coordinates || max_x < 0.0 || max_y < 0.0 || - min_x >= measurement_frame_info.getWidth() - 0.5 || min_y >= measurement_frame_info.getHeight() - 0.5) { + if (bad_coordinates || max_coord.m_x < 0 || max_coord.m_y < 0 || + min_coord.m_x >= measurement_frame_info.getWidth() || min_coord.m_y >= measurement_frame_info.getHeight()) { source.setIndexedProperty(m_instance, bad_coordinates); } // Clip the coordinates to fit the available image else { - min_x = std::max(0.0, min_x); - min_y = std::max(0.0, min_y); - max_x = std::min(double(measurement_frame_info.getWidth()) - 0.5, max_x); - max_y = std::min(double(measurement_frame_info.getHeight()) - 0.5, max_y); + min_coord.m_x = std::max(0, min_coord.m_x); + min_coord.m_y = std::max(0, min_coord.m_y); + max_coord.m_x = std::min(measurement_frame_info.getWidth() - 1, max_coord.m_x); + max_coord.m_y = std::min(measurement_frame_info.getHeight() - 1, max_coord.m_y); source.setIndexedProperty( - m_instance, ImageCoordinate(min_x, min_y), ImageCoordinate(max_x, max_y)); + m_instance, min_coord, max_coord, ImageCoordinate(min_x, min_y), ImageCoordinate(max_x, max_y)); } } diff --git a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp index c994d17b..2cf63a2c 100644 --- a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp +++ b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp @@ -64,20 +64,26 @@ void MeasurementFrameRectangleTaskNoDetect::computeProperties(SourceInterface& s auto max_x = std::max(coord1.m_x, std::max(coord2.m_x, std::max(coord3.m_x, coord4.m_x))); auto max_y = std::max(coord1.m_y, std::max(coord2.m_y, std::max(coord3.m_y, coord4.m_y))); + PixelCoordinate min_coord, max_coord; + min_coord.m_x = int(min_x); + min_coord.m_y = int(min_y); + max_coord.m_x = int(max_x) + 1; + max_coord.m_y = int(max_y) + 1; + // The full boundaries may lie outside of the frame - if (bad_coordinates || max_x < 0.0 || max_y < 0.0 || - min_x >= measurement_frame_info.getWidth() - 0.5 || min_y >= measurement_frame_info.getHeight() - 0.5) { + if (bad_coordinates || max_coord.m_x < 0 || max_coord.m_y < 0 || + min_coord.m_x >= measurement_frame_info.getWidth() || min_coord.m_y >= measurement_frame_info.getHeight()) { source.setIndexedProperty(m_instance, bad_coordinates); } // Clip the coordinates to fit the available image else { - min_x = std::max(0.0, min_x); - min_y = std::max(0.0, min_y); - max_x = std::min(double(measurement_frame_info.getWidth()) - 0.5, max_x); - max_y = std::min(double(measurement_frame_info.getHeight()) - 0.5, max_y); + min_coord.m_x = std::max(0, min_coord.m_x); + min_coord.m_y = std::max(0, min_coord.m_y); + max_coord.m_x = std::min(measurement_frame_info.getWidth() - 1, max_coord.m_x); + max_coord.m_y = std::min(measurement_frame_info.getHeight() - 1, max_coord.m_y); source.setIndexedProperty( - m_instance, ImageCoordinate(min_x, min_y), ImageCoordinate(max_x, max_y)); + m_instance, min_coord, max_coord, ImageCoordinate(min_x, min_y), ImageCoordinate(max_x, max_y)); } } From 1990c86e70415338363fe404583b08a7c3b29d10 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Fri, 29 Aug 2025 18:02:57 +0200 Subject: [PATCH 08/40] fix helgrind problems --- .../SEFramework/Output/OutputRegistry.h | 2 +- SEFramework/src/lib/Output/OutputRegistry.cpp | 19 ++++++---- .../Measurement/MultithreadedMeasurement.h | 8 ++--- .../Measurement/MultithreadedMeasurement.cpp | 36 +++++++++---------- 4 files changed, 34 insertions(+), 31 deletions(-) diff --git a/SEFramework/SEFramework/Output/OutputRegistry.h b/SEFramework/SEFramework/Output/OutputRegistry.h index 83c5b0ae..1a6a0d88 100644 --- a/SEFramework/SEFramework/Output/OutputRegistry.h +++ b/SEFramework/SEFramework/Output/OutputRegistry.h @@ -169,7 +169,7 @@ class OutputRegistry { return converter(source.getProperty(i)); }; } - Euclid::Table::Row::cell_type operator()(const SourceInterface& source) { + Euclid::Table::Row::cell_type operator()(const SourceInterface& source) const { return m_convert_func(source, index); } std::size_t index = 0; diff --git a/SEFramework/src/lib/Output/OutputRegistry.cpp b/SEFramework/src/lib/Output/OutputRegistry.cpp index 6cef389b..cf284cdf 100644 --- a/SEFramework/src/lib/Output/OutputRegistry.cpp +++ b/SEFramework/src/lib/Output/OutputRegistry.cpp @@ -49,21 +49,26 @@ auto OutputRegistry::getSourceToRowConverter(const std::vector& ena return [this, out_prop_list](const SourceInterface& source) { std::vector info_list {}; std::vector cell_values {}; + + const auto& property_to_names_map = m_property_to_names_map; + const auto& name_to_col_info_map = m_name_to_col_info_map; + const auto& name_to_converter_map = m_name_to_converter_map; + for (const auto& property : out_prop_list) { - if (m_property_to_names_map.count(property) == 0) { + if (property_to_names_map.count(property) == 0) { throw Elements::Exception() << "Missing column generator for " << property.name(); } - for (const auto& name : m_property_to_names_map.at(property)) { - auto& col_info = m_name_to_col_info_map.at(name); - info_list.emplace_back(name, m_name_to_converter_map.at(name).first, - col_info.unit, col_info.description); - cell_values.emplace_back(m_name_to_converter_map.at(name).second(source)); + for (const auto& name : property_to_names_map.at(property)) { + auto& col_info = name_to_col_info_map.at(name); + info_list.push_back(ColumnInfo::info_type {name, name_to_converter_map.at(name).first, + col_info.unit, col_info.description}); + cell_values.push_back(name_to_converter_map.at(name).second(source)); } } if (info_list.empty()) { throw Elements::Exception() << "The given configuration would not generate any output"; } - return Row {std::move(cell_values), std::make_shared(move(info_list))}; + return Row {std::move(cell_values), std::make_shared(std::move(info_list))}; }; } diff --git a/SEImplementation/SEImplementation/Measurement/MultithreadedMeasurement.h b/SEImplementation/SEImplementation/Measurement/MultithreadedMeasurement.h index 7e46a139..3524303d 100644 --- a/SEImplementation/SEImplementation/Measurement/MultithreadedMeasurement.h +++ b/SEImplementation/SEImplementation/Measurement/MultithreadedMeasurement.h @@ -40,11 +40,11 @@ class MultithreadedMeasurement : public Measurement { using SourceToRowConverter = std::function; MultithreadedMeasurement(SourceToRowConverter source_to_row, const std::shared_ptr& thread_pool, - unsigned max_queue_size) + unsigned /*max_queue_size*/) : m_source_to_row(source_to_row), m_thread_pool(thread_pool), m_group_counter(0), - m_input_done(false), m_abort_raised(false), m_semaphore(max_queue_size) {} + m_input_done(false), m_abort_raised(false) {} ~MultithreadedMeasurement() override; @@ -66,10 +66,8 @@ class MultithreadedMeasurement : public Measurement { int m_group_counter; std::atomic_bool m_input_done, m_abort_raised; - std::condition_variable m_new_output; std::list>> m_output_queue; - std::mutex m_output_queue_mutex; - Euclid::Semaphore m_semaphore; + mutable std::mutex m_output_queue_mutex; }; } diff --git a/SEImplementation/src/lib/Measurement/MultithreadedMeasurement.cpp b/SEImplementation/src/lib/Measurement/MultithreadedMeasurement.cpp index b29fdecf..65ee15fc 100644 --- a/SEImplementation/src/lib/Measurement/MultithreadedMeasurement.cpp +++ b/SEImplementation/src/lib/Measurement/MultithreadedMeasurement.cpp @@ -34,7 +34,7 @@ static Elements::Logging logger = Elements::Logging::getLogger("Multithreading") MultithreadedMeasurement::~MultithreadedMeasurement() { - if (m_output_thread->joinable()) { + if (m_output_thread && m_output_thread->joinable()) { m_output_thread->join(); } } @@ -46,7 +46,9 @@ void MultithreadedMeasurement::startThreads() { void MultithreadedMeasurement::stopThreads() { m_input_done = true; m_thread_pool->block(); - m_output_thread->join(); + if (m_output_thread && m_output_thread->joinable()) { + m_output_thread->join(); + } logger.debug() << "All worker threads done!"; } @@ -81,6 +83,7 @@ void MultithreadedMeasurement::receiveSource(std::unique_ptr output_lock(m_output_queue_mutex); + std::lock_guard output_lock(m_output_queue_mutex); m_output_queue.emplace_back(order_number, std::move(source_group)); } - m_new_output.notify_one(); }; auto lambda_copyable = [lambda = std::make_shared(std::move(lambda))](){ (*lambda)(); }; m_thread_pool->submit(lambda_copyable); + ++m_group_counter; } @@ -117,24 +120,21 @@ void MultithreadedMeasurement::outputThreadStatic(MultithreadedMeasurement *meas } void MultithreadedMeasurement::outputThreadLoop() { - while (m_thread_pool->activeThreads() > 0) { - std::unique_lock output_lock(m_output_queue_mutex); - - // Wait for something in the output queue - if (m_output_queue.empty()) { - m_new_output.wait_for(output_lock, std::chrono::milliseconds(100)); - } + while (true) { + { + std::lock_guard output_lock(m_output_queue_mutex); - // Process the output queue - while (!m_output_queue.empty()) { - sendSource(std::move(m_output_queue.front().second)); - m_output_queue.pop_front(); - } + while (!m_output_queue.empty()) { + sendSource(std::move(m_output_queue.front().second)); + m_output_queue.pop_front(); + } - if (m_input_done && m_thread_pool->running() + m_thread_pool->queued() == 0 && + if (m_input_done && m_thread_pool->running() + m_thread_pool->queued() == 0 && m_output_queue.empty()) { - break; + break; + } } + std::this_thread::sleep_for(std::chrono::milliseconds(5)); } } From 661b05cf9a1f175a5aac3fd7d0d27ed14be3cd63 Mon Sep 17 00:00:00 2001 From: mkuemmel Date: Fri, 5 Sep 2025 10:57:35 +0200 Subject: [PATCH 09/40] Includes a nan and inf test that works with fast-math Uses it for checking PSFEx files. --- .../SEFramework/Psf/SEPP_isnan_isinf.h | 76 +++++++++++++++++++ SEFramework/SEFramework/Psf/VariablePsf.h | 2 + SEFramework/src/lib/Psf/VariablePsf.cpp | 15 ++++ 3 files changed, 93 insertions(+) create mode 100644 SEFramework/SEFramework/Psf/SEPP_isnan_isinf.h diff --git a/SEFramework/SEFramework/Psf/SEPP_isnan_isinf.h b/SEFramework/SEFramework/Psf/SEPP_isnan_isinf.h new file mode 100644 index 00000000..248f84a6 --- /dev/null +++ b/SEFramework/SEFramework/Psf/SEPP_isnan_isinf.h @@ -0,0 +1,76 @@ +#ifndef _SEIMPLEMENTATION_PSF_SEPP_ISNAN_ISINF_H +#define _SEIMPLEMENTATION_PSF_SEPP_ISNAN_ISINF_H +//#include + +/* + * 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 sepp_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 sepp_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 sepp_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 sepp_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); +}; + +#endif diff --git a/SEFramework/SEFramework/Psf/VariablePsf.h b/SEFramework/SEFramework/Psf/VariablePsf.h index f50bda38..57a6f4fa 100644 --- a/SEFramework/SEFramework/Psf/VariablePsf.h +++ b/SEFramework/SEFramework/Psf/VariablePsf.h @@ -24,6 +24,8 @@ #ifndef _SEIMPLEMENTATION_PSF_VARIABLEPSF_H_ #define _SEIMPLEMENTATION_PSF_VARIABLEPSF_H_ +#include + #include "SEFramework/Image/VectorImage.h" #include "SEFramework/Psf/Psf.h" diff --git a/SEFramework/src/lib/Psf/VariablePsf.cpp b/SEFramework/src/lib/Psf/VariablePsf.cpp index c8367c35..62bfc7cd 100644 --- a/SEFramework/src/lib/Psf/VariablePsf.cpp +++ b/SEFramework/src/lib/Psf/VariablePsf.cpp @@ -21,10 +21,13 @@ * Author: Alejandro Álvarez Ayllón */ +#include #include #include #include "SEFramework/Psf/VariablePsf.h" +#include "SEFramework/Psf/SEPP_isnan_isinf.h" +static auto stack_logger = Elements::Logging::getLogger("PSFExPsf"); namespace SourceXtractor { @@ -94,6 +97,8 @@ void VariablePsf::selfTest() { if (m_coefficients.size() == 0) { throw Elements::Exception() << "A variable PSF needs at least one set of coefficients"; } + // give some feedback + stack_logger.info() << "In VariablePsf::selfTest()"; // Pre-condition: There is a degree value per unique group std::vector n_component_per_group(m_group_degrees.size()); @@ -130,6 +135,16 @@ void VariablePsf::selfTest() { if (coeff->getWidth() != psf_width || coeff->getHeight() != psf_height) { throw Elements::Exception() << "Malformed variable PSF, coefficient matrices do not have the same dimensions"; } + for (auto x = 0; x < psf_width; ++x) { + for (auto y = 0; y < psf_height; ++y) { + if (sepp_isnan(coeff->at(x, y))) { + throw Elements::Exception() << "Malformed variable PSF, coefficient matrices contains NANs"; + } + else if (sepp_isnan(coeff->at(x, y))) { + throw Elements::Exception() << "Malformed variable PSF, coefficient matrices contains INFs"; + } + } + } } } From cf0f9c80895f7e3423cd624966ad3ba467d6a1ba Mon Sep 17 00:00:00 2001 From: mkuemmel Date: Sat, 20 Sep 2025 23:08:41 +0200 Subject: [PATCH 10/40] Pushed some feedback to debug --- SEFramework/src/lib/Psf/VariablePsf.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SEFramework/src/lib/Psf/VariablePsf.cpp b/SEFramework/src/lib/Psf/VariablePsf.cpp index 62bfc7cd..51d9c243 100644 --- a/SEFramework/src/lib/Psf/VariablePsf.cpp +++ b/SEFramework/src/lib/Psf/VariablePsf.cpp @@ -98,7 +98,7 @@ void VariablePsf::selfTest() { throw Elements::Exception() << "A variable PSF needs at least one set of coefficients"; } // give some feedback - stack_logger.info() << "In VariablePsf::selfTest()"; + stack_logger.debug() << "In VariablePsf::selfTest()"; // Pre-condition: There is a degree value per unique group std::vector n_component_per_group(m_group_degrees.size()); From e9c75415da0ab7d13ee6b602c198463573cf5f1f Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Fri, 26 Sep 2025 15:14:49 +0200 Subject: [PATCH 11/40] fix regression with exceptions in threads --- .../src/lib/Measurement/MultithreadedMeasurement.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SEImplementation/src/lib/Measurement/MultithreadedMeasurement.cpp b/SEImplementation/src/lib/Measurement/MultithreadedMeasurement.cpp index 65ee15fc..b0d7046a 100644 --- a/SEImplementation/src/lib/Measurement/MultithreadedMeasurement.cpp +++ b/SEImplementation/src/lib/Measurement/MultithreadedMeasurement.cpp @@ -120,7 +120,7 @@ void MultithreadedMeasurement::outputThreadStatic(MultithreadedMeasurement *meas } void MultithreadedMeasurement::outputThreadLoop() { - while (true) { + while (m_thread_pool->activeThreads() > 0) { { std::lock_guard output_lock(m_output_queue_mutex); From a092f362fab6482817a273109bb8cfb764a95e1d Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Fri, 26 Sep 2025 16:46:52 +0200 Subject: [PATCH 12/40] moved isnan function to SEUtils, small fixes --- SEFramework/SEFramework/Psf/VariablePsf.h | 2 -- SEFramework/src/lib/Psf/VariablePsf.cpp | 25 +++++++++---------- .../SEUtils/IsNan.h | 20 +++++++++------ 3 files changed, 24 insertions(+), 23 deletions(-) rename SEFramework/SEFramework/Psf/SEPP_isnan_isinf.h => SEUtils/SEUtils/IsNan.h (89%) diff --git a/SEFramework/SEFramework/Psf/VariablePsf.h b/SEFramework/SEFramework/Psf/VariablePsf.h index 57a6f4fa..f50bda38 100644 --- a/SEFramework/SEFramework/Psf/VariablePsf.h +++ b/SEFramework/SEFramework/Psf/VariablePsf.h @@ -24,8 +24,6 @@ #ifndef _SEIMPLEMENTATION_PSF_VARIABLEPSF_H_ #define _SEIMPLEMENTATION_PSF_VARIABLEPSF_H_ -#include - #include "SEFramework/Image/VectorImage.h" #include "SEFramework/Psf/Psf.h" diff --git a/SEFramework/src/lib/Psf/VariablePsf.cpp b/SEFramework/src/lib/Psf/VariablePsf.cpp index 51d9c243..6e66e98a 100644 --- a/SEFramework/src/lib/Psf/VariablePsf.cpp +++ b/SEFramework/src/lib/Psf/VariablePsf.cpp @@ -21,11 +21,12 @@ * Author: Alejandro Álvarez Ayllón */ +#include #include #include -#include +#include "SEUtils/IsNan.h" + #include "SEFramework/Psf/VariablePsf.h" -#include "SEFramework/Psf/SEPP_isnan_isinf.h" static auto stack_logger = Elements::Logging::getLogger("PSFExPsf"); @@ -97,8 +98,6 @@ void VariablePsf::selfTest() { if (m_coefficients.size() == 0) { throw Elements::Exception() << "A variable PSF needs at least one set of coefficients"; } - // give some feedback - stack_logger.debug() << "In VariablePsf::selfTest()"; // Pre-condition: There is a degree value per unique group std::vector n_component_per_group(m_group_degrees.size()); @@ -135,15 +134,15 @@ void VariablePsf::selfTest() { if (coeff->getWidth() != psf_width || coeff->getHeight() != psf_height) { throw Elements::Exception() << "Malformed variable PSF, coefficient matrices do not have the same dimensions"; } - for (auto x = 0; x < psf_width; ++x) { - for (auto y = 0; y < psf_height; ++y) { - if (sepp_isnan(coeff->at(x, y))) { - throw Elements::Exception() << "Malformed variable PSF, coefficient matrices contains NANs"; - } - else if (sepp_isnan(coeff->at(x, y))) { - throw Elements::Exception() << "Malformed variable PSF, coefficient matrices contains INFs"; - } - } + for (auto x = 0; x < psf_width; ++x){ + for (auto y = 0; y < psf_height; ++y) { + if (fastmath_isnan(coeff->at(x, y))) { + throw Elements::Exception() << "Malformed variable PSF, coefficient matrices contains NANs"; + } + else if (fastmath_isinf(coeff->at(x, y))){ + throw Elements::Exception() << "Malformed variable PSF, coefficient matrices contains INFs"; + } + } } } } diff --git a/SEFramework/SEFramework/Psf/SEPP_isnan_isinf.h b/SEUtils/SEUtils/IsNan.h similarity index 89% rename from SEFramework/SEFramework/Psf/SEPP_isnan_isinf.h rename to SEUtils/SEUtils/IsNan.h index 248f84a6..c4dea647 100644 --- a/SEFramework/SEFramework/Psf/SEPP_isnan_isinf.h +++ b/SEUtils/SEUtils/IsNan.h @@ -1,6 +1,9 @@ -#ifndef _SEIMPLEMENTATION_PSF_SEPP_ISNAN_ISINF_H -#define _SEIMPLEMENTATION_PSF_SEPP_ISNAN_ISINF_H -//#include +#ifndef _SEUTILS_ISNAN_H +#define _SEUTILS_ISNAN_H + +#include + +namespace SourceXtractor { /* * From: @@ -31,7 +34,7 @@ const unsigned FLOAT_EXP_MASK = 0x7F800000; // 00000000011111111111111111111111 const unsigned FLOAT_FRAC_PART_MASK = 0x7FFFFF; -inline bool sepp_isnan(float x) { +inline bool fastmath_isnan(float x) { union { uint32_t u; float f; @@ -40,7 +43,7 @@ inline bool sepp_isnan(float x) { return ((conv.u & FLOAT_EXP_MASK) == FLOAT_EXP_MASK) && ((conv.u & FLOAT_FRAC_PART_MASK) != 0); }; -inline bool sepp_isinf(float x) { +inline bool fastmath_isinf(float x) { union { uint32_t u; float f; @@ -54,7 +57,7 @@ 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 sepp_isnan(double x) { +inline bool fastmath_isnan(double x) { union { uint64_t u; double f; @@ -63,8 +66,7 @@ inline bool sepp_isnan(double x) { return ((conv.u & DOUBLE_EXP_MASK) == DOUBLE_EXP_MASK) && ((conv.u & DOUBLE_FRAC_PART_MASK) != 0); }; - -inline bool sepp_isinf(double x) { +inline bool fastmath_isinf(double x) { union { uint64_t u; double f; @@ -73,4 +75,6 @@ inline bool sepp_isinf(double x) { return ((conv.u & DOUBLE_EXP_MASK) == DOUBLE_EXP_MASK) && ((conv.u & DOUBLE_FRAC_PART_MASK) == 0); }; +} // namespace SourceXtractor + #endif From 6aea35fabb0ff6201b51baa8ca81ac6bd27e7e49 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Fri, 26 Sep 2025 17:03:04 +0200 Subject: [PATCH 13/40] print filename of PSF that doesn't load --- SEImplementation/src/lib/Plugin/Psf/PsfPluginConfig.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/SEImplementation/src/lib/Plugin/Psf/PsfPluginConfig.cpp b/SEImplementation/src/lib/Plugin/Psf/PsfPluginConfig.cpp index 462835f9..640cfa65 100644 --- a/SEImplementation/src/lib/Plugin/Psf/PsfPluginConfig.cpp +++ b/SEImplementation/src/lib/Plugin/Psf/PsfPluginConfig.cpp @@ -124,10 +124,11 @@ static std::shared_ptr readPsfEx(std::unique_ptr &pFi } return std::make_shared(pixel_sampling, components, group_degrees, coefficients); - } catch (CCfits::FITS::NoSuchHDU&) { - throw; } catch (CCfits::FitsException &e) { throw Elements::Exception() << "Error loading PSFEx file: " << e.message(); + } catch (...) { + logger.error() << "Failed loading a psf file: " << pFits->name(); + throw; } } From 13905a58a51228cdada6406de7cf167188b66b84 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Mon, 29 Sep 2025 14:56:35 +0200 Subject: [PATCH 14/40] fix regression --- SEImplementation/src/lib/Plugin/Psf/PsfPluginConfig.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/SEImplementation/src/lib/Plugin/Psf/PsfPluginConfig.cpp b/SEImplementation/src/lib/Plugin/Psf/PsfPluginConfig.cpp index 640cfa65..e0c61a9d 100644 --- a/SEImplementation/src/lib/Plugin/Psf/PsfPluginConfig.cpp +++ b/SEImplementation/src/lib/Plugin/Psf/PsfPluginConfig.cpp @@ -124,6 +124,8 @@ static std::shared_ptr readPsfEx(std::unique_ptr &pFi } return std::make_shared(pixel_sampling, components, group_degrees, coefficients); + } catch (CCfits::FITS::NoSuchHDU&) { // Make sure we propagate this specific exception + throw; } catch (CCfits::FitsException &e) { throw Elements::Exception() << "Error loading PSFEx file: " << e.message(); } catch (...) { From 8f86e9cfe0ea9cd05d9d33d9039d75412b9c4fce Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Fri, 31 Oct 2025 17:28:03 +0100 Subject: [PATCH 15/40] added option to specify max nb of simultaneous files --- .../SEFramework/FITS/FitsImageSource.h | 7 ++- SEFramework/src/lib/FITS/FitsImageSource.cpp | 19 +++++- .../Configuration/FileManagerConfig.h | 44 ++++++++++++++ .../lib/Configuration/FileManagerConfig.cpp | 59 +++++++++++++++++++ SEMain/src/program/SourceXtractor.cpp | 8 ++- 5 files changed, 131 insertions(+), 6 deletions(-) create mode 100644 SEImplementation/SEImplementation/Configuration/FileManagerConfig.h create mode 100644 SEImplementation/src/lib/Configuration/FileManagerConfig.cpp diff --git a/SEFramework/SEFramework/FITS/FitsImageSource.h b/SEFramework/SEFramework/FITS/FitsImageSource.h index 5686df98..a66107d6 100644 --- a/SEFramework/SEFramework/FITS/FitsImageSource.h +++ b/SEFramework/SEFramework/FITS/FitsImageSource.h @@ -46,6 +46,7 @@ using Euclid::FilePool::FileHandler; class FitsImageSource : public ImageSource, public std::enable_shared_from_this { public: + static std::shared_ptr getFileManager(unsigned int max_open_files = 500); /** * Constructor @@ -57,14 +58,14 @@ class FitsImageSource : public ImageSource, public std::enable_shared_from_this< */ explicit FitsImageSource(const std::string& filename, int hdu_number = 0, ImageTile::ImageType image_type = ImageTile::AutoType, - std::shared_ptr manager = FileManager::getDefault()); + std::shared_ptr manager = getFileManager()); FitsImageSource(const std::string& filename, int width, int height, ImageTile::ImageType image_type, const std::shared_ptr coord_system = nullptr, bool append = false, bool empty_primary = false, - std::shared_ptr manager = FileManager::getDefault()); + std::shared_ptr manager = getFileManager()); virtual ~FitsImageSource() = default; @@ -137,6 +138,8 @@ class FitsImageSource : public ImageSource, public std::enable_shared_from_this< ImageTile::ImageType m_image_type; int m_current_layer; + + static std::shared_ptr s_file_manager; }; } diff --git a/SEFramework/src/lib/FITS/FitsImageSource.cpp b/SEFramework/src/lib/FITS/FitsImageSource.cpp index f68b0629..776ac420 100644 --- a/SEFramework/src/lib/FITS/FitsImageSource.cpp +++ b/SEFramework/src/lib/FITS/FitsImageSource.cpp @@ -23,10 +23,13 @@ */ #include "SEFramework/FITS/FitsImageSource.h" + #include "SEFramework/FITS/FitsFile.h" #include "SEUtils/VariantCast.h" #include +#include #include +#include "ElementsKernel/Logging.h" #include #include #include @@ -39,6 +42,9 @@ namespace SourceXtractor { +static Elements::Logging logger = Elements::Logging::getLogger("FitsFile"); + + using Euclid::make_unique; namespace { @@ -386,7 +392,18 @@ int FitsImageSource::getImageType() const { return LONGLONG_IMG; } } - //FIXME add missing types + +std::shared_ptr FitsImageSource::getFileManager(unsigned int max_open_files) { + if (!s_file_manager) { + logger.info() << "Creating shared LRUFileManager with limit of " << max_open_files << " open files."; + s_file_manager = std::make_shared(max_open_files); + } + return s_file_manager; } + +std::shared_ptr FitsImageSource::s_file_manager = nullptr; + +} // namespace SourceXtractor + diff --git a/SEImplementation/SEImplementation/Configuration/FileManagerConfig.h b/SEImplementation/SEImplementation/Configuration/FileManagerConfig.h new file mode 100644 index 00000000..08f1e7d8 --- /dev/null +++ b/SEImplementation/SEImplementation/Configuration/FileManagerConfig.h @@ -0,0 +1,44 @@ +/** Copyright © 2019-2025 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 + */ + +#ifndef SEIMPLEMENTATION_CONFIGURATION_FILEMANAGERCONFIG_H +#define SEIMPLEMENTATION_CONFIGURATION_FILEMANAGERCONFIG_H + +#include "Configuration/Configuration.h" + +namespace SourceXtractor { + +class FileManagerConfig : public Euclid::Configuration::Configuration { +public: + explicit FileManagerConfig(long manager_id); + virtual ~FileManagerConfig() = default; + + std::map getProgramOptions() override; + void preInitialize(const UserValues& args) override; + void initialize(const UserValues& args) override; + + int getMaxSimultaneousFiles() const { + return m_max_simultaneous_files; + } + +private: + int m_max_simultaneous_files; +}; + +} // namespace SourceXtractor + +#endif // SEIMPLEMENTATION_CONFIGURATION_FILEMANAGERCONFIG_H \ No newline at end of file diff --git a/SEImplementation/src/lib/Configuration/FileManagerConfig.cpp b/SEImplementation/src/lib/Configuration/FileManagerConfig.cpp new file mode 100644 index 00000000..4117b40e --- /dev/null +++ b/SEImplementation/src/lib/Configuration/FileManagerConfig.cpp @@ -0,0 +1,59 @@ +/** Copyright © 2019-2025 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 "SEImplementation/Configuration/FileManagerConfig.h" + +#include + +#include "ElementsKernel/Exception.h" +#include "Configuration/ConfigManager.h" + +#include "SEFramework/FITS/FitsImageSource.h" + +#include + +using namespace Euclid::Configuration; +namespace po = boost::program_options; + +namespace SourceXtractor { + +static const std::string MAX_SIMULTANEOUS_FILES {"max-simultaneous-files" }; + + + // Constructor +FileManagerConfig::FileManagerConfig(long manager_id) + : Configuration(manager_id), m_max_simultaneous_files(500) { +} + +std::map FileManagerConfig::getProgramOptions() { + std::cout << "Getting FileManagerConfig options" << std::endl; + return { {"File Manager", { + {MAX_SIMULTANEOUS_FILES.c_str(), po::value()->default_value(500), + "Maximum number of simultaneous open files."} + }} }; +} + +void FileManagerConfig::preInitialize(const UserValues& args) { + // We have to configure the FileManager before any other configuration tries to use it + m_max_simultaneous_files = args.find(MAX_SIMULTANEOUS_FILES)->second.as(); + FitsImageSource::getFileManager(m_max_simultaneous_files); +} + +void FileManagerConfig::initialize(const UserValues& args) { +} + +} // namespace SourceXtractor \ No newline at end of file diff --git a/SEMain/src/program/SourceXtractor.cpp b/SEMain/src/program/SourceXtractor.cpp index 49b2d493..e79ae077 100644 --- a/SEMain/src/program/SourceXtractor.cpp +++ b/SEMain/src/program/SourceXtractor.cpp @@ -74,6 +74,7 @@ #include "SEImplementation/Configuration/MemoryConfig.h" #include "SEImplementation/Configuration/OutputConfig.h" #include "SEImplementation/Configuration/SamplingConfig.h" +#include "SEImplementation/Configuration/FileManagerConfig.h" #include "SEImplementation/CheckImages/CheckImages.h" #include "SEImplementation/Prefetcher/Prefetcher.h" @@ -191,6 +192,7 @@ class SEMain : public Elements::Program { config_manager.registerConfiguration(); config_manager.registerConfiguration(); config_manager.registerConfiguration(); + config_manager.registerConfiguration(); CheckImages::getInstance().reportConfigDependencies(config_manager); @@ -319,15 +321,15 @@ class SEMain : public Elements::Program { throw Elements::Exception() << "The configuration file '" << cfg_file << "' does not exist"; } } - + // Create the progress listener and printer ASAP progress_printer_factory.configure(args); auto progress_mediator = progress_printer_factory.createProgressMediator(); - + // Initialize the rest of the components auto& config_manager = ConfigManager::getInstance(config_manager_id); config_manager.initialize(args); - + // Configure TileManager auto memory_config = config_manager.getConfiguration(); TileManager::getInstance()->setOptions(memory_config.getTileSize(), From 03cd7765195c545e0d44ff42f46e6771f3966a37 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Wed, 3 Dec 2025 16:48:25 +0100 Subject: [PATCH 16/40] remove leftover debug print --- SEImplementation/src/lib/Configuration/FileManagerConfig.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/SEImplementation/src/lib/Configuration/FileManagerConfig.cpp b/SEImplementation/src/lib/Configuration/FileManagerConfig.cpp index 4117b40e..e37c683c 100644 --- a/SEImplementation/src/lib/Configuration/FileManagerConfig.cpp +++ b/SEImplementation/src/lib/Configuration/FileManagerConfig.cpp @@ -40,7 +40,6 @@ FileManagerConfig::FileManagerConfig(long manager_id) } std::map FileManagerConfig::getProgramOptions() { - std::cout << "Getting FileManagerConfig options" << std::endl; return { {"File Manager", { {MAX_SIMULTANEOUS_FILES.c_str(), po::value()->default_value(500), "Maximum number of simultaneous open files."} From 22da3e6293b1c03baeeb5b6f8bfbad3b64de2800 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Fri, 20 Feb 2026 14:19:03 +0100 Subject: [PATCH 17/40] fix Dependent parameters not removing observers in destructor --- .../ModelFitting/Parameters/DependentParameter.h | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/ModelFitting/ModelFitting/Parameters/DependentParameter.h b/ModelFitting/ModelFitting/Parameters/DependentParameter.h index c759a110..8dc9e8c3 100644 --- a/ModelFitting/ModelFitting/Parameters/DependentParameter.h +++ b/ModelFitting/ModelFitting/Parameters/DependentParameter.h @@ -65,7 +65,11 @@ class DependentParameter: public BasicParameter { //m_get_value_hook = std::bind(&DependentParameter::getValueHook, this); } - virtual ~DependentParameter() = default; + virtual ~DependentParameter() { + for (auto& parameter_observer : m_parameter_observers) { + std::get<0>(parameter_observer)->removeObserver(std::get<1>(parameter_observer)); + } + } double getValue() const override { if (!this->isObserved()) { @@ -118,13 +122,17 @@ class DependentParameter: public BasicParameter { template void addParameterObserver(int, Param& param) { - param->addObserver([this](double){ + auto id = param->addObserver([this](double){ // Do not bother updating live if there are no observers if (this->isObserved()) { this->update((*m_params)[0]->getValue()); } }); + + m_parameter_observers.emplace_back(param, id); } + + std::vector, size_t>> m_parameter_observers; }; template @@ -133,6 +141,7 @@ std::shared_ptr> createDependentParameter( return std::make_shared>(value_calculator, parameters...); } + } #endif /* MODELFITTING_DEPENDENTPARAMETER_H */ From 0e2d8dda13e4cea9ff4dfe324b2262bdbea403c1 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Fri, 20 Feb 2026 14:50:07 +0100 Subject: [PATCH 18/40] avoid recreating the model fitting parameter object multiple times --- .../FlexibleModelFittingIterativeTask.h | 8 +- .../FlexibleModelFittingIterativeTask.cpp | 96 +++++++------------ 2 files changed, 41 insertions(+), 63 deletions(-) diff --git a/SEImplementation/SEImplementation/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.h b/SEImplementation/SEImplementation/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.h index 9c0ecc3c..839c2ef5 100644 --- a/SEImplementation/SEImplementation/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.h +++ b/SEImplementation/SEImplementation/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.h @@ -21,6 +21,7 @@ #include "ModelFitting/Models/FrameModel.h" #include "ModelFitting/Engine/ResidualEstimator.h" #include "ModelFitting/Engine/LeastSquareEngineManager.h" +#include "ModelFitting/Engine/EngineParameterManager.h" #include "SEUtils/PixelRectangle.h" @@ -31,6 +32,7 @@ #include "SEImplementation/Plugin/FlexibleModelFitting/FlexibleModelFittingParameter.h" #include "SEImplementation/Plugin/FlexibleModelFitting/FlexibleModelFittingFrame.h" #include "SEImplementation/Plugin/FlexibleModelFitting/FlexibleModelFittingPrior.h" +#include "SEImplementation/Plugin/FlexibleModelFitting/FlexibleModelFittingParameterManager.h" #include "SEImplementation/Image/DownSampledImagePsf.h" @@ -85,6 +87,10 @@ class FlexibleModelFittingIterativeTask : public GroupTask { std::vector iterations_per_meta; std::vector fitting_areas_x; std::vector fitting_areas_y; + + FlexibleModelFittingParameterManager parameter_manager; + ModelFitting::EngineParameterManager engine_parameter_manager {}; + int n_free_parameters = 0; }; struct FittingState { @@ -115,7 +121,7 @@ class FlexibleModelFittingIterativeTask : public GroupTask { std::shared_ptr> model, std::shared_ptr> weights, int& data_points) const; int fitSourcePrepareParameters(FlexibleModelFittingParameterManager& parameter_manager, ModelFitting::EngineParameterManager& engine_parameter_manager, - SourceInterface& source, int index, FittingState& state) const; + SourceInterface& source, SourceState& state) const; int fitSourcePrepareModels(FlexibleModelFittingParameterManager& parameter_manager, ModelFitting::ResidualEstimator& res_estimator, int& good_pixels, SourceGroupInterface& group, SourceInterface& source, int index, FittingState& state, double downscaling) const; diff --git a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp index 0c1425d5..5a669b5b 100644 --- a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp +++ b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp @@ -417,6 +417,8 @@ void FlexibleModelFittingIterativeTask::computeProperties(SourceGroupInterface& } } + initial_state.n_free_parameters = fitSourcePrepareParameters(initial_state.parameter_manager, initial_state.engine_parameter_manager, source, initial_state); + fitting_state.source_states.emplace_back(std::move(initial_state)); } @@ -492,45 +494,27 @@ std::shared_ptr> FlexibleModelFittingIterativeTask::createD auto rect = getFittingRect(source, frame_index); double pixel_scale = 1.0; - FlexibleModelFittingParameterManager parameter_manager; - ModelFitting::EngineParameterManager engine_parameter_manager {}; - int n_free_parameters = 0; + auto deblend_image = VectorImage::create(rect.getWidth(), rect.getHeight()); int index = 0; for (auto& src : group) { if (index != source_index) { + // reset parameters to final values after fitting for (auto parameter : m_parameters) { - auto free_parameter = std::dynamic_pointer_cast(parameter); - - if (free_parameter != nullptr) { - ++n_free_parameters; - - // Initial with the values from the current iteration run - parameter_manager.addParameter(src, parameter, - free_parameter->create(parameter_manager, engine_parameter_manager, src, - state.source_states[index].parameters_initial_values.at(free_parameter->getId()), - state.source_states[index].parameters_values.at(free_parameter->getId()))); - } else { - parameter_manager.addParameter(src, parameter, - parameter->create(parameter_manager, engine_parameter_manager, src)); + auto engine_parameter = std::dynamic_pointer_cast(state.source_states[index].parameter_manager.getParameter(src, parameter)); + if (engine_parameter != nullptr) { + engine_parameter->setValue(state.source_states[index].parameters_values.at(parameter->getId())); } } - } - index++; - } - auto deblend_image = VectorImage::create(rect.getWidth(), rect.getHeight()); - index = 0; - for (auto& src : group) { - if (index != source_index) { - auto frame_model = createFrameModel(src, pixel_scale, parameter_manager, frame, rect); - auto final_stamp = frame_model.getImage(); + auto frame_model = createFrameModel(src, pixel_scale, state.source_states[index].parameter_manager, frame, rect); + auto final_stamp = frame_model.getImage(); - for (int y = 0; y < final_stamp->getHeight(); ++y) { - for (int x = 0; x < final_stamp->getWidth(); ++x) { - deblend_image->at(x, y) += final_stamp->at(x, y); - } + for (int y = 0; y < final_stamp->getHeight(); ++y) { + for (int x = 0; x < final_stamp->getWidth(); ++x) { + deblend_image->at(x, y) += final_stamp->at(x, y); } + } } index++; } @@ -541,7 +525,7 @@ std::shared_ptr> FlexibleModelFittingIterativeTask::createD int FlexibleModelFittingIterativeTask::fitSourcePrepareParameters( FlexibleModelFittingParameterManager& parameter_manager, ModelFitting::EngineParameterManager& engine_parameter_manager, - SourceInterface& source, int index, FittingState& state) const { + SourceInterface& source, SourceState& state) const { int free_parameters_nb = 0; for (auto parameter : m_parameters) { auto free_parameter = std::dynamic_pointer_cast(parameter); @@ -552,8 +536,8 @@ int FlexibleModelFittingIterativeTask::fitSourcePrepareParameters( // Initial with the values from the current iteration run parameter_manager.addParameter(source, parameter, free_parameter->create(parameter_manager, engine_parameter_manager, source, - state.source_states[index].parameters_initial_values.at(free_parameter->getId()), - state.source_states[index].parameters_values.at(free_parameter->getId()))); + state.parameters_initial_values.at(free_parameter->getId()), + state.parameters_values.at(free_parameter->getId()))); } else { parameter_manager.addParameter(source, parameter, parameter->create(parameter_manager, engine_parameter_manager, source)); @@ -708,17 +692,19 @@ void FlexibleModelFittingIterativeTask::fitSource(SourceGroupInterface& group, S ////////////////////////////////////////////// // Prepare parameters - FlexibleModelFittingParameterManager parameter_manager; - ModelFitting::EngineParameterManager engine_parameter_manager{}; - int n_free_parameters = fitSourcePrepareParameters( - parameter_manager, engine_parameter_manager, source, index, state); + for (auto parameter : m_parameters) { + auto engine_parameter = std::dynamic_pointer_cast(state.source_states[index].parameter_manager.getParameter(source, parameter)); + if (engine_parameter != nullptr) { + engine_parameter->setValue(state.source_states[index].parameters_values.at(parameter->getId())); + } + } /////////////////////////////////////////////////////////////////////////////////// // Add models for all frames ResidualEstimator res_estimator {}; int n_good_pixels = 0; int valid_frames = fitSourcePrepareModels( - parameter_manager, res_estimator, n_good_pixels, group, source, index, state, down_scaling); + state.source_states[index].parameter_manager, res_estimator, n_good_pixels, group, source, index, state, down_scaling); /////////////////////////////////////////////////////////////////////////////// // Check that we had enough data for the fit @@ -728,7 +714,7 @@ void FlexibleModelFittingIterativeTask::fitSource(SourceGroupInterface& group, S if (valid_frames == 0) { flags = Flags::OUTSIDE; } - else if (n_good_pixels < n_free_parameters) { + else if (n_good_pixels < state.source_states[index].n_free_parameters) { flags = Flags::INSUFFICIENT_DATA; } @@ -745,14 +731,14 @@ void FlexibleModelFittingIterativeTask::fitSource(SourceGroupInterface& group, S //////////////////////////////////////////////////////////////////////////////// // Add priors for (auto prior : m_priors) { - prior->setupPrior(parameter_manager, source, res_estimator); + prior->setupPrior(state.source_states[index].parameter_manager, source, res_estimator); } ///////////////////////////////////////////////////////////////////////////////// // Model fitting auto engine = LeastSquareEngineManager::create(m_least_squares_engine, m_max_iterations); - auto solution = engine->solveProblem(engine_parameter_manager, res_estimator); + auto solution = engine->solveProblem(state.source_states[index].engine_parameter_manager, res_estimator); auto iterations = solution.iteration_no; auto stop_reason = solution.engine_stop_reason; @@ -764,49 +750,34 @@ void FlexibleModelFittingIterativeTask::fitSource(SourceGroupInterface& group, S //////////////////////////////////////////////////////////////////////////////////// // compute chi squared - SeFloat avg_reduced_chi_squared = fitSourceComputeChiSquared(parameter_manager, group, source, index, state); + SeFloat avg_reduced_chi_squared = fitSourceComputeChiSquared(state.source_states[index].parameter_manager, group, source, index, state); //////////////////////////////////////////////////////////////////////////////////// // update state with results - fitSourceUpdateState(parameter_manager, source, avg_reduced_chi_squared, duration, iterations, stop_reason, flags, solution, + fitSourceUpdateState(state.source_states[index].parameter_manager, source, avg_reduced_chi_squared, duration, iterations, stop_reason, flags, solution, index, state); } void FlexibleModelFittingIterativeTask::updateCheckImages(SourceGroupInterface& group, double pixel_scale, FittingState& state) const { - // recreate parameters - - FlexibleModelFittingParameterManager parameter_manager; - ModelFitting::EngineParameterManager engine_parameter_manager {}; - int index = 0; for (auto& src : group) { for (auto parameter : m_parameters) { - auto free_parameter = std::dynamic_pointer_cast(parameter); - - if (free_parameter != nullptr) { - // Initialize with the values from the current iteration run - parameter_manager.addParameter(src, parameter, - free_parameter->create(parameter_manager, engine_parameter_manager, src, - state.source_states[index].parameters_initial_values.at(free_parameter->getId()), - state.source_states[index].parameters_values.at(free_parameter->getId()))); - } else { - parameter_manager.addParameter(src, parameter, - parameter->create(parameter_manager, engine_parameter_manager, src)); + // reset parameters to final values after fitting + auto engine_parameter = std::dynamic_pointer_cast(state.source_states[index].parameter_manager.getParameter(src, parameter)); + if (engine_parameter != nullptr) { + engine_parameter->setValue(state.source_states[index].parameters_values.at(parameter->getId())); } } - index++; - } - for (auto& src : group) { for (auto frame : m_frames) { int frame_index = frame->getFrameNb(); if (isFrameValid(src, frame_index)) { auto stamp_rect = getFittingRect(src, frame_index); - auto frame_model = createFrameModel(src, pixel_scale, parameter_manager, frame, stamp_rect); + auto frame_model = createFrameModel(src, pixel_scale, state.source_states[index].parameter_manager, frame, stamp_rect); auto final_stamp = frame_model.getImage(); auto weight_image = createWeightImage(src, frame_index); @@ -846,6 +817,7 @@ void FlexibleModelFittingIterativeTask::updateCheckImages(SourceGroupInterface& } } + index++; } } From bcb72281967f6a4acc60bf359d7291718aaf92d7 Mon Sep 17 00:00:00 2001 From: Marc Schefer <31961290+marcschefer@users.noreply.github.com> Date: Fri, 6 Mar 2026 16:17:52 +0100 Subject: [PATCH 19/40] Revert "Changed fitting window calculation to use doubles" --- .../CoordinateSystem/CoordinateSystem.h | 24 ----------------- .../MeasurementFrameRectangle.h | 18 ++++--------- .../FlexibleModelFittingIterativeTask.cpp | 16 +++++------- .../MeasurementFrameRectangleTask.cpp | 3 +-- .../MeasurementFrameRectangleTaskNoDetect.cpp | 26 +++++++------------ 5 files changed, 22 insertions(+), 65 deletions(-) diff --git a/SEFramework/SEFramework/CoordinateSystem/CoordinateSystem.h b/SEFramework/SEFramework/CoordinateSystem/CoordinateSystem.h index 30605da8..d43e26a4 100644 --- a/SEFramework/SEFramework/CoordinateSystem/CoordinateSystem.h +++ b/SEFramework/SEFramework/CoordinateSystem/CoordinateSystem.h @@ -45,30 +45,6 @@ struct ImageCoordinate { ImageCoordinate() : m_x(0), m_y(0) {} ImageCoordinate(double x, double y) : m_x(x), m_y(y) {} - - ImageCoordinate operator*(double scalar) const { - return ImageCoordinate(m_x * scalar, m_y * scalar); - } - - ImageCoordinate operator+(const ImageCoordinate& other) const { - return ImageCoordinate(m_x + other.m_x, m_y + other.m_y); - } - - ImageCoordinate& operator+=(const ImageCoordinate& other) { - m_x += other.m_x; - m_y += other.m_y; - return *this; - } - - ImageCoordinate operator-(const ImageCoordinate& other) const { - return ImageCoordinate(m_x - other.m_x, m_y - other.m_y); - } - - ImageCoordinate& operator-=(const ImageCoordinate& other) { - m_x -= other.m_x; - m_y -= other.m_y; - return *this; - } }; class CoordinateSystem { diff --git a/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h b/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h index d163e632..fb7ae567 100644 --- a/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h +++ b/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h @@ -24,10 +24,7 @@ #ifndef _SEIMPLEMENTATION_PLUGIN_MEASUREMENTFRAMERECTANGLE_MEASUREMENTFRAMERECTANGLE_H_ #define _SEIMPLEMENTATION_PLUGIN_MEASUREMENTFRAMERECTANGLE_MEASUREMENTFRAMERECTANGLE_H_ -#include - #include "SEUtils/PixelRectangle.h" -#include "SEFramework/CoordinateSystem/CoordinateSystem.h" #include "SEFramework/Property/Property.h" #include "SEFramework/Image/Image.h" @@ -39,15 +36,11 @@ class MeasurementFrameRectangle: public Property { virtual ~MeasurementFrameRectangle() = default; explicit MeasurementFrameRectangle(bool bad_projection): - m_min_coord{-1, -1}, m_max_coord{-1, -1}, m_min_coord_image{-1, -1}, - m_max_coord_image{-1, -1}, m_bad_projection{bad_projection} {} + m_min_coord{-1, -1}, m_max_coord{-1, -1}, m_bad_projection{bad_projection}{} - MeasurementFrameRectangle(PixelCoordinate min_coord, PixelCoordinate max_coord, - ImageCoordinate min_coord_image, ImageCoordinate max_coord_image): - m_min_coord(min_coord), m_max_coord(max_coord), - m_min_coord_image(min_coord_image), m_max_coord_image(max_coord_image), m_bad_projection{false} { + MeasurementFrameRectangle(PixelCoordinate min_coord, PixelCoordinate max_coord): + m_min_coord{min_coord}, m_max_coord{max_coord}, m_bad_projection{false} { assert(min_coord.m_x <= max_coord.m_x && min_coord.m_y <= max_coord.m_y); - assert(min_coord_image.m_x <= max_coord_image.m_x && min_coord_image.m_y <= max_coord_image.m_y); } PixelCoordinate getTopLeft() const { @@ -72,8 +65,8 @@ class MeasurementFrameRectangle: public Property { return m_max_coord.m_y - m_min_coord.m_y + 1; } - std::tuple getImageRect() const { - return std::make_tuple(m_min_coord_image, m_max_coord_image); + PixelRectangle getRect() const { + return PixelRectangle(m_min_coord, m_max_coord); } bool badProjection() const { @@ -82,7 +75,6 @@ class MeasurementFrameRectangle: public Property { private: PixelCoordinate m_min_coord, m_max_coord; - ImageCoordinate m_min_coord_image, m_max_coord_image; bool m_bad_projection; }; diff --git a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp index 0c1425d5..1e99d9ea 100644 --- a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp +++ b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp @@ -74,8 +74,7 @@ FlexibleModelFittingIterativeTask::~FlexibleModelFittingIterativeTask() { } PixelRectangle FlexibleModelFittingIterativeTask::getUnclippedFittingRect(SourceInterface& source, int frame_index) const { - ImageCoordinate min_coord, max_coord; - std::tie(min_coord, max_coord) = source.getProperty(frame_index).getImageRect(); + auto fitting_rect = source.getProperty(frame_index).getRect(); if (m_window_type == WindowType::ROTATED_ELLIPSE) { auto ellipse = getFittingEllipse(source, frame_index); @@ -85,14 +84,14 @@ PixelRectangle FlexibleModelFittingIterativeTask::getUnclippedFittingRect(Source return getEllipseRect(ellipse); } - if (max_coord.m_x - min_coord.m_x <= 0 || max_coord.m_y - min_coord.m_y <= 0) { + if (fitting_rect.getWidth() <= 0 || fitting_rect.getHeight() <= 0) { return PixelRectangle(); } else { - auto min = min_coord; - auto max = max_coord; + auto min = fitting_rect.getTopLeft(); + auto max = fitting_rect.getBottomRight(); // FIXME temporary, for now just enlarge the area by a fixed amount of pixels - ImageCoordinate border = (max - min) * .8 + ImageCoordinate(2.0, 2.0); + PixelCoordinate border = (max - min) * .8 + PixelCoordinate(2, 2); min -= border; max += border; @@ -126,10 +125,7 @@ PixelRectangle FlexibleModelFittingIterativeTask::getUnclippedFittingRect(Source max.m_y = min.m_y + size; } - auto min_pc = PixelCoordinate(static_cast(min.m_x + 0.5), static_cast(min.m_y + 0.5)); - auto max_pc = PixelCoordinate(static_cast(max.m_x + 0.5), static_cast(max.m_y + 0.5)); - - return PixelRectangle(min_pc, max_pc); + return PixelRectangle(min, max); } } diff --git a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp index da88a8dc..f57d56fb 100644 --- a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp +++ b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp @@ -85,8 +85,7 @@ void MeasurementFrameRectangleTask::computeProperties(SourceInterface& source) c max_coord.m_x = std::min(measurement_frame_info.getWidth() - 1, max_coord.m_x); max_coord.m_y = std::min(measurement_frame_info.getHeight() - 1, max_coord.m_y); - source.setIndexedProperty( - m_instance, min_coord, max_coord, ImageCoordinate(min_x, min_y), ImageCoordinate(max_x, max_y)); + source.setIndexedProperty(m_instance, min_coord, max_coord); } } diff --git a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp index 2cf63a2c..68caedfc 100644 --- a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp +++ b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp @@ -23,7 +23,6 @@ #include "SEImplementation/Plugin/DetectionFrameCoordinates/DetectionFrameCoordinates.h" #include #include -#include #include #include @@ -32,27 +31,23 @@ namespace SourceXtractor { void MeasurementFrameRectangleTaskNoDetect::computeProperties(SourceInterface& source) const { auto measurement_frame_coordinates = source.getProperty(m_instance).getCoordinateSystem(); - auto reference_frame_coordinates = source.getProperty().getCoordinateSystem(); - const auto& measurement_frame_info = source.getProperty(m_instance); const auto& world_centroid = source.getProperty(); const auto& assoc_mode = source.getProperty(); + auto coord = world_centroid.getCentroid(); + bool bad_coordinates = false; ImageCoordinate coord1, coord2, coord3, coord4; try { - auto w = assoc_mode.getRefFramePixelWidth(); - auto h = assoc_mode.getRefFramePixelHeight(); + int w = assoc_mode.getRefFramePixelWidth(); + int h = assoc_mode.getRefFramePixelHeight(); - auto c = reference_frame_coordinates->worldToImage(world_centroid.getCentroid()); - coord1 = measurement_frame_coordinates->worldToImage( - reference_frame_coordinates->imageToWorld(ImageCoordinate(c.m_x - w, c.m_y - h))); - coord2 = measurement_frame_coordinates->worldToImage( - reference_frame_coordinates->imageToWorld(ImageCoordinate(c.m_x + w, c.m_y - h))); - coord3 = measurement_frame_coordinates->worldToImage( - reference_frame_coordinates->imageToWorld(ImageCoordinate(c.m_x - w, c.m_y + h))); - coord4 = measurement_frame_coordinates->worldToImage( - reference_frame_coordinates->imageToWorld(ImageCoordinate(c.m_x + w, c.m_y + h))); + auto c = measurement_frame_coordinates->worldToImage(coord); + coord1 = ImageCoordinate(c.m_x - w, c.m_y - h); + coord2 = ImageCoordinate(c.m_x + w, c.m_y - h); + coord3 = ImageCoordinate(c.m_x - w, c.m_y + h); + coord4 = ImageCoordinate(c.m_x + w, c.m_y + h); } catch (const InvalidCoordinatesException&) { bad_coordinates = true; @@ -82,8 +77,7 @@ void MeasurementFrameRectangleTaskNoDetect::computeProperties(SourceInterface& s max_coord.m_x = std::min(measurement_frame_info.getWidth() - 1, max_coord.m_x); max_coord.m_y = std::min(measurement_frame_info.getHeight() - 1, max_coord.m_y); - source.setIndexedProperty( - m_instance, min_coord, max_coord, ImageCoordinate(min_x, min_y), ImageCoordinate(max_x, max_y)); + source.setIndexedProperty(m_instance, min_coord, max_coord); } } From 634e46195924d81c7a4c5c5f3df4df0069ca239a Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Wed, 18 Mar 2026 15:26:22 +0100 Subject: [PATCH 20/40] propagate bad coords flag --- .../FlexibleModelFittingIterativeTask.cpp | 12 +++++++++++- .../MeasurementFrameRectangleTaskNoDetect.cpp | 7 ++++++- 2 files changed, 17 insertions(+), 2 deletions(-) diff --git a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp index 0c1425d5..24524e11 100644 --- a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp +++ b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp @@ -75,7 +75,13 @@ FlexibleModelFittingIterativeTask::~FlexibleModelFittingIterativeTask() { PixelRectangle FlexibleModelFittingIterativeTask::getUnclippedFittingRect(SourceInterface& source, int frame_index) const { ImageCoordinate min_coord, max_coord; - std::tie(min_coord, max_coord) = source.getProperty(frame_index).getImageRect(); + auto& measurement_frame_rectangle = source.getProperty(frame_index); + std::tie(min_coord, max_coord) = measurement_frame_rectangle.getImageRect(); + + if (measurement_frame_rectangle.badProjection()) { + return PixelRectangle(); + } + if (m_window_type == WindowType::ROTATED_ELLIPSE) { auto ellipse = getFittingEllipse(source, frame_index); @@ -135,6 +141,10 @@ PixelRectangle FlexibleModelFittingIterativeTask::getUnclippedFittingRect(Source PixelRectangle FlexibleModelFittingIterativeTask::clipFittingRect(PixelRectangle fitting_rect, SourceInterface& source, int frame_index) const { + if (fitting_rect.getWidth() <= 0 || fitting_rect.getHeight() <= 0) { + return PixelRectangle(); + } + const auto& frame_info = source.getProperty(frame_index); auto min = fitting_rect.getTopLeft(); diff --git a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp index 2cf63a2c..eaf08456 100644 --- a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp +++ b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp @@ -73,7 +73,7 @@ void MeasurementFrameRectangleTaskNoDetect::computeProperties(SourceInterface& s // The full boundaries may lie outside of the frame if (bad_coordinates || max_coord.m_x < 0 || max_coord.m_y < 0 || min_coord.m_x >= measurement_frame_info.getWidth() || min_coord.m_y >= measurement_frame_info.getHeight()) { - source.setIndexedProperty(m_instance, bad_coordinates); + source.setIndexedProperty(m_instance, true); } // Clip the coordinates to fit the available image else { @@ -82,6 +82,11 @@ void MeasurementFrameRectangleTaskNoDetect::computeProperties(SourceInterface& s max_coord.m_x = std::min(measurement_frame_info.getWidth() - 1, max_coord.m_x); max_coord.m_y = std::min(measurement_frame_info.getHeight() - 1, max_coord.m_y); + min_x = std::max(0.0, min_x); + min_y = std::max(0.0, min_y); + max_x = std::min(measurement_frame_info.getWidth() - 1.0, max_x); + max_y = std::min(measurement_frame_info.getHeight() - 1.0, max_y); + source.setIndexedProperty( m_instance, min_coord, max_coord, ImageCoordinate(min_x, min_y), ImageCoordinate(max_x, max_y)); } From 74f5e863355300d5508beb4000fb95e985ad7b20 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Thu, 19 Mar 2026 17:18:50 +0100 Subject: [PATCH 21/40] change pearson's formula constant for mode estimation --- .../src/lib/Background/SE/ImageMode.cpp | 13 ++----------- 1 file changed, 2 insertions(+), 11 deletions(-) diff --git a/SEImplementation/src/lib/Background/SE/ImageMode.cpp b/SEImplementation/src/lib/Background/SE/ImageMode.cpp index 857c1bb4..208aa513 100644 --- a/SEImplementation/src/lib/Background/SE/ImageMode.cpp +++ b/SEImplementation/src/lib/Background/SE/ImageMode.cpp @@ -106,17 +106,8 @@ std::tuple ImageMode::getBackGuess(const std::vector& data) const { // Sigma is 0 T mode; - if (std::abs(sigma) == 0) { - mode = mean; - } - // Not crowded: mean and median do not differ more than 30% - else if (std::abs((mean - median) / sigma) < 0.3) { - mode = 2.5 * median - 1.5 * mean; - } - // Crowded case: we use the median - else { - mode = median; - } + mode = 3.5 * median - 2.5 * mean; + return std::make_tuple(mode, sigma); } From f7bcc14b38336b2e2a97af969d2323fec1353729 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Fri, 20 Mar 2026 16:57:52 +0100 Subject: [PATCH 22/40] remove int coordinates completely --- .../MeasurementFrameRectangle.h | 23 ++++++++----------- .../MeasurementFrameRectangleTask.cpp | 20 ++++++---------- .../MeasurementFrameRectangleTaskNoDetect.cpp | 23 +++++-------------- 3 files changed, 22 insertions(+), 44 deletions(-) diff --git a/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h b/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h index d163e632..dab1ff02 100644 --- a/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h +++ b/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h @@ -39,23 +39,19 @@ class MeasurementFrameRectangle: public Property { virtual ~MeasurementFrameRectangle() = default; explicit MeasurementFrameRectangle(bool bad_projection): - m_min_coord{-1, -1}, m_max_coord{-1, -1}, m_min_coord_image{-1, -1}, - m_max_coord_image{-1, -1}, m_bad_projection{bad_projection} {} + m_min_coord{-1, -1}, m_max_coord{-1, -1}, m_bad_projection{bad_projection} {} - MeasurementFrameRectangle(PixelCoordinate min_coord, PixelCoordinate max_coord, - ImageCoordinate min_coord_image, ImageCoordinate max_coord_image): - m_min_coord(min_coord), m_max_coord(max_coord), - m_min_coord_image(min_coord_image), m_max_coord_image(max_coord_image), m_bad_projection{false} { + MeasurementFrameRectangle(ImageCoordinate min_coord, ImageCoordinate max_coord): + m_min_coord(min_coord), m_max_coord(max_coord), m_bad_projection{false} { assert(min_coord.m_x <= max_coord.m_x && min_coord.m_y <= max_coord.m_y); - assert(min_coord_image.m_x <= max_coord_image.m_x && min_coord_image.m_y <= max_coord_image.m_y); } - PixelCoordinate getTopLeft() const { + ImageCoordinate getTopLeft() const { assert(m_max_coord.m_x >= 0); return m_min_coord; } - PixelCoordinate getBottomRight() const { + ImageCoordinate getBottomRight() const { assert(m_max_coord.m_x >= 0); return m_max_coord; } @@ -63,17 +59,17 @@ class MeasurementFrameRectangle: public Property { int getWidth() const { if (m_max_coord.m_x < 0) return 0; - return m_max_coord.m_x - m_min_coord.m_x + 1; + return m_max_coord.m_x - m_min_coord.m_x; } int getHeight() const { if (m_max_coord.m_x < 0) return 0; - return m_max_coord.m_y - m_min_coord.m_y + 1; + return m_max_coord.m_y - m_min_coord.m_y; } std::tuple getImageRect() const { - return std::make_tuple(m_min_coord_image, m_max_coord_image); + return std::make_tuple(m_min_coord, m_max_coord); } bool badProjection() const { @@ -81,8 +77,7 @@ class MeasurementFrameRectangle: public Property { } private: - PixelCoordinate m_min_coord, m_max_coord; - ImageCoordinate m_min_coord_image, m_max_coord_image; + ImageCoordinate m_min_coord, m_max_coord; bool m_bad_projection; }; diff --git a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp index da88a8dc..4c3226a7 100644 --- a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp +++ b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp @@ -67,26 +67,20 @@ void MeasurementFrameRectangleTask::computeProperties(SourceInterface& source) c auto max_x = std::max(coord1.m_x, std::max(coord2.m_x, std::max(coord3.m_x, coord4.m_x))); auto max_y = std::max(coord1.m_y, std::max(coord2.m_y, std::max(coord3.m_y, coord4.m_y))); - PixelCoordinate min_coord, max_coord; - min_coord.m_x = int(min_x); - min_coord.m_y = int(min_y); - max_coord.m_x = int(max_x) + 1; - max_coord.m_y = int(max_y) + 1; - // The full boundaries may lie outside of the frame - if (bad_coordinates || max_coord.m_x < 0 || max_coord.m_y < 0 || - min_coord.m_x >= measurement_frame_info.getWidth() || min_coord.m_y >= measurement_frame_info.getHeight()) { + if (bad_coordinates || max_x < 0.0 || max_y < 0.0 || + int(min_x) >= measurement_frame_info.getWidth() || int(max_y) >= measurement_frame_info.getHeight()) { source.setIndexedProperty(m_instance, bad_coordinates); } // Clip the coordinates to fit the available image else { - min_coord.m_x = std::max(0, min_coord.m_x); - min_coord.m_y = std::max(0, min_coord.m_y); - max_coord.m_x = std::min(measurement_frame_info.getWidth() - 1, max_coord.m_x); - max_coord.m_y = std::min(measurement_frame_info.getHeight() - 1, max_coord.m_y); + min_x = std::max(0.0, min_x); + min_y = std::max(0.0, min_y); + max_x = std::min(static_cast(measurement_frame_info.getWidth()), max_x); + max_y = std::min(static_cast(measurement_frame_info.getHeight()), max_y); source.setIndexedProperty( - m_instance, min_coord, max_coord, ImageCoordinate(min_x, min_y), ImageCoordinate(max_x, max_y)); + m_instance, ImageCoordinate(min_x, min_y), ImageCoordinate(max_x, max_y)); } } diff --git a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp index eaf08456..fdbe4c46 100644 --- a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp +++ b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp @@ -64,31 +64,20 @@ void MeasurementFrameRectangleTaskNoDetect::computeProperties(SourceInterface& s auto max_x = std::max(coord1.m_x, std::max(coord2.m_x, std::max(coord3.m_x, coord4.m_x))); auto max_y = std::max(coord1.m_y, std::max(coord2.m_y, std::max(coord3.m_y, coord4.m_y))); - PixelCoordinate min_coord, max_coord; - min_coord.m_x = int(min_x); - min_coord.m_y = int(min_y); - max_coord.m_x = int(max_x) + 1; - max_coord.m_y = int(max_y) + 1; - // The full boundaries may lie outside of the frame - if (bad_coordinates || max_coord.m_x < 0 || max_coord.m_y < 0 || - min_coord.m_x >= measurement_frame_info.getWidth() || min_coord.m_y >= measurement_frame_info.getHeight()) { - source.setIndexedProperty(m_instance, true); + if (bad_coordinates || max_x < 0.0 || max_y < 0.0 || + int(min_x) >= measurement_frame_info.getWidth() || int(max_y) >= measurement_frame_info.getHeight()) { + source.setIndexedProperty(m_instance, bad_coordinates); } // Clip the coordinates to fit the available image else { - min_coord.m_x = std::max(0, min_coord.m_x); - min_coord.m_y = std::max(0, min_coord.m_y); - max_coord.m_x = std::min(measurement_frame_info.getWidth() - 1, max_coord.m_x); - max_coord.m_y = std::min(measurement_frame_info.getHeight() - 1, max_coord.m_y); - min_x = std::max(0.0, min_x); min_y = std::max(0.0, min_y); - max_x = std::min(measurement_frame_info.getWidth() - 1.0, max_x); - max_y = std::min(measurement_frame_info.getHeight() - 1.0, max_y); + max_x = std::min(static_cast(measurement_frame_info.getWidth()), max_x); + max_y = std::min(static_cast(measurement_frame_info.getHeight()), max_y); source.setIndexedProperty( - m_instance, min_coord, max_coord, ImageCoordinate(min_x, min_y), ImageCoordinate(max_x, max_y)); + m_instance, ImageCoordinate(min_x, min_y), ImageCoordinate(max_x, max_y)); } } From 5c8fde557966297cca102dd9b2723ab9cdc08a47 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Fri, 20 Mar 2026 17:06:40 +0100 Subject: [PATCH 23/40] Add isValid() to MeasurementFrameRectangle --- .../MeasurementFrameRectangle.h | 16 ++++++++++------ .../FlexibleModelFittingIterativeTask.cpp | 2 +- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h b/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h index dab1ff02..5753eced 100644 --- a/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h +++ b/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h @@ -47,23 +47,23 @@ class MeasurementFrameRectangle: public Property { } ImageCoordinate getTopLeft() const { - assert(m_max_coord.m_x >= 0); + assert(m_max_coord.m_x >= 0.0); return m_min_coord; } ImageCoordinate getBottomRight() const { - assert(m_max_coord.m_x >= 0); + assert(m_max_coord.m_x >= 0.0); return m_max_coord; } - int getWidth() const { - if (m_max_coord.m_x < 0) + double getWidth() const { + if (m_max_coord.m_x < 0.0) return 0; return m_max_coord.m_x - m_min_coord.m_x; } - int getHeight() const { - if (m_max_coord.m_x < 0) + double getHeight() const { + if (m_max_coord.m_x < 0.0) return 0; return m_max_coord.m_y - m_min_coord.m_y; } @@ -76,6 +76,10 @@ class MeasurementFrameRectangle: public Property { return m_bad_projection; } + bool isValid() const { + return !m_bad_projection && m_max_coord.m_x >= 0.0 && m_max_coord.m_y >= 0.0 && getWidth() > 0.0 && getHeight() > 0.0; + } + private: ImageCoordinate m_min_coord, m_max_coord; bool m_bad_projection; diff --git a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp index 24524e11..755834d8 100644 --- a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp +++ b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp @@ -78,7 +78,7 @@ PixelRectangle FlexibleModelFittingIterativeTask::getUnclippedFittingRect(Source auto& measurement_frame_rectangle = source.getProperty(frame_index); std::tie(min_coord, max_coord) = measurement_frame_rectangle.getImageRect(); - if (measurement_frame_rectangle.badProjection()) { + if (!measurement_frame_rectangle.isValid()) { return PixelRectangle(); } From 89d017443cfc17903615bbf9efc7cb52cfa20895 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Fri, 27 Mar 2026 16:36:27 +0100 Subject: [PATCH 24/40] stricter checks --- .../MeasurementFrameRectangle.h | 18 +++++++++++------- .../MeasurementFrameRectangleTask.cpp | 14 ++++++++++---- .../MeasurementFrameRectangleTaskNoDetect.cpp | 8 ++++---- .../SourceFlags/SourceFlagsSourceTask.cpp | 2 +- 4 files changed, 26 insertions(+), 16 deletions(-) diff --git a/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h b/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h index 5753eced..dbd27895 100644 --- a/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h +++ b/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h @@ -47,24 +47,28 @@ class MeasurementFrameRectangle: public Property { } ImageCoordinate getTopLeft() const { - assert(m_max_coord.m_x >= 0.0); + assert(m_min_coord.m_x >= 0.0); + assert(m_min_coord.m_y >= 0.0); return m_min_coord; } ImageCoordinate getBottomRight() const { assert(m_max_coord.m_x >= 0.0); + assert(m_max_coord.m_y >= 0.0); return m_max_coord; } double getWidth() const { - if (m_max_coord.m_x < 0.0) - return 0; + if (m_max_coord.m_x <= 0.0) { + return 0.0; + } return m_max_coord.m_x - m_min_coord.m_x; } double getHeight() const { - if (m_max_coord.m_x < 0.0) - return 0; + if (m_max_coord.m_y <= 0.0) { + return 0.0; + } return m_max_coord.m_y - m_min_coord.m_y; } @@ -72,12 +76,12 @@ class MeasurementFrameRectangle: public Property { return std::make_tuple(m_min_coord, m_max_coord); } - bool badProjection() const { + bool isBadProjection() const { return m_bad_projection; } bool isValid() const { - return !m_bad_projection && m_max_coord.m_x >= 0.0 && m_max_coord.m_y >= 0.0 && getWidth() > 0.0 && getHeight() > 0.0; + return !m_bad_projection && m_max_coord.m_x > 0.0 && m_max_coord.m_y > 0.0 && getWidth() > 0.0 && getHeight() > 0.0; } private: diff --git a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp index 4c3226a7..1a635308 100644 --- a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp +++ b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp @@ -44,9 +44,15 @@ void MeasurementFrameRectangleTask::computeProperties(SourceInterface& source) c auto width = detection_group_stamp.getWidth(); auto height = detection_group_stamp.getHeight(); + if (width <= 0 || height <= 0) { + // No valid stamp, set the property to an empty rectangle + source.setIndexedProperty(m_instance, false); + return; + } + // Transform the 4 corner coordinates from detection image ImageCoordinate coord1, coord2, coord3, coord4; - bool bad_coordinates = false; + bool bad_projection = false; try { coord1 = measurement_frame_coordinates->worldToImage( detection_frame_coordinates->imageToWorld(ImageCoordinate(stamp_top_left.m_x, stamp_top_left.m_y))); @@ -58,7 +64,7 @@ void MeasurementFrameRectangleTask::computeProperties(SourceInterface& source) c detection_frame_coordinates->imageToWorld(ImageCoordinate(stamp_top_left.m_x, stamp_top_left.m_y + height))); } catch (const InvalidCoordinatesException&) { - bad_coordinates = true; + bad_projection = true; } // Determine the min/max coordinates @@ -68,9 +74,9 @@ void MeasurementFrameRectangleTask::computeProperties(SourceInterface& source) c auto max_y = std::max(coord1.m_y, std::max(coord2.m_y, std::max(coord3.m_y, coord4.m_y))); // The full boundaries may lie outside of the frame - if (bad_coordinates || max_x < 0.0 || max_y < 0.0 || + if (bad_projection || max_x < 0.0 || max_y < 0.0 || int(min_x) >= measurement_frame_info.getWidth() || int(max_y) >= measurement_frame_info.getHeight()) { - source.setIndexedProperty(m_instance, bad_coordinates); + source.setIndexedProperty(m_instance, bad_projection); } // Clip the coordinates to fit the available image else { diff --git a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp index fdbe4c46..f5757b05 100644 --- a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp +++ b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp @@ -38,7 +38,7 @@ void MeasurementFrameRectangleTaskNoDetect::computeProperties(SourceInterface& s const auto& world_centroid = source.getProperty(); const auto& assoc_mode = source.getProperty(); - bool bad_coordinates = false; + bool bad_projection = false; ImageCoordinate coord1, coord2, coord3, coord4; try { auto w = assoc_mode.getRefFramePixelWidth(); @@ -55,7 +55,7 @@ void MeasurementFrameRectangleTaskNoDetect::computeProperties(SourceInterface& s reference_frame_coordinates->imageToWorld(ImageCoordinate(c.m_x + w, c.m_y + h))); } catch (const InvalidCoordinatesException&) { - bad_coordinates = true; + bad_projection = true; } // Determine the min/max coordinates @@ -65,9 +65,9 @@ void MeasurementFrameRectangleTaskNoDetect::computeProperties(SourceInterface& s auto max_y = std::max(coord1.m_y, std::max(coord2.m_y, std::max(coord3.m_y, coord4.m_y))); // The full boundaries may lie outside of the frame - if (bad_coordinates || max_x < 0.0 || max_y < 0.0 || + if (bad_projection || max_x < 0.0 || max_y < 0.0 || int(min_x) >= measurement_frame_info.getWidth() || int(max_y) >= measurement_frame_info.getHeight()) { - source.setIndexedProperty(m_instance, bad_coordinates); + source.setIndexedProperty(m_instance, bad_projection); } // Clip the coordinates to fit the available image else { diff --git a/SEImplementation/src/lib/Plugin/SourceFlags/SourceFlagsSourceTask.cpp b/SEImplementation/src/lib/Plugin/SourceFlags/SourceFlagsSourceTask.cpp index 00598a19..a86461a0 100644 --- a/SEImplementation/src/lib/Plugin/SourceFlags/SourceFlagsSourceTask.cpp +++ b/SEImplementation/src/lib/Plugin/SourceFlags/SourceFlagsSourceTask.cpp @@ -39,7 +39,7 @@ void SourceFlagsSourceTask::computeProperties(SourceInterface &source) const { source_flag |= Flags::SATURATED * source.getProperty(instance).getSaturateFlag(); source_flag |= Flags::BOUNDARY * source.getProperty(instance).getBoundaryFlag(); source_flag |= Flags::BLENDED * source.getProperty().getBlendedFlag(); - if (source.getProperty(instance).badProjection()) { + if (source.getProperty(instance).isBadProjection()) { source_flag |= Flags::BAD_PROJECTION; } } From 588bb42b6d8c3110737ea500fe87bde2e1f92dc4 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Tue, 31 Mar 2026 13:39:41 +0200 Subject: [PATCH 25/40] remove FindBoostPython.cmake, it should not be needed anymore --- SEImplementation/CMakeLists.txt | 7 +++---- cmake/modules/FindBoostPython.cmake | 2 +- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/SEImplementation/CMakeLists.txt b/SEImplementation/CMakeLists.txt index 1b52a93d..49981ff1 100644 --- a/SEImplementation/CMakeLists.txt +++ b/SEImplementation/CMakeLists.txt @@ -33,8 +33,7 @@ elements_depends_on_subdirs(Pyston) find_package(PythonInterp ${PYTHON_EXPLICIT_VERSION} REQUIRED) find_package(PythonLibs ${PYTHON_EXPLICIT_VERSION} REQUIRED) -find_package(BoostPython ${PYTHON_EXPLICIT_VERSION}) -find_package(Boost 1.53 REQUIRED) +find_package(Boost 1.53 REQUIRED COMPONENTS python${PYTHON_VERSION_MAJOR}${PYTHON_VERSION_MINOR}) find_package(OnnxRuntime) find_package(Log4CPP REQUIRED) @@ -101,10 +100,10 @@ elements_add_library(SEImplementation LINK_LIBRARIES ElementsKernel Configuration Table SEUtils SEFramework ModelFitting - BoostPython PythonLibs PythonInterp + PythonLibs PythonInterp Pyston ${OPT_LIBRARIES} - INCLUDE_DIRS ${BoostPython_INCLUDE_DIRS} ${PYTHON_INCLUDE_DIRS} ${OPT_INCLUDES} + INCLUDE_DIRS ${PYTHON_INCLUDE_DIRS} ${OPT_INCLUDES} PUBLIC_HEADERS SEImplementation) #=============================================================================== diff --git a/cmake/modules/FindBoostPython.cmake b/cmake/modules/FindBoostPython.cmake index 69fbcc59..dd19374b 100644 --- a/cmake/modules/FindBoostPython.cmake +++ b/cmake/modules/FindBoostPython.cmake @@ -1,6 +1,6 @@ # Normally, finding a boost component should rely on # find_package(Boost REQUIRED COMPONENTS python) -# But the packaging is now always the same: sometimes it will be python3, +# But the packaging is not always the same: sometimes it will be python3, # others python37, python2, python27, or just python # (i.e. python3 for Fedora < 30, but python37 for Fedora >= 30 and MacOSX via Homebrew) # We wrap all this in this module From c7559ec346771ae181eb66a48f5eac23db29fa31 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Tue, 31 Mar 2026 13:48:01 +0200 Subject: [PATCH 26/40] actually delete the file --- cmake/modules/FindBoostPython.cmake | 38 ----------------------------- 1 file changed, 38 deletions(-) delete mode 100644 cmake/modules/FindBoostPython.cmake diff --git a/cmake/modules/FindBoostPython.cmake b/cmake/modules/FindBoostPython.cmake deleted file mode 100644 index dd19374b..00000000 --- a/cmake/modules/FindBoostPython.cmake +++ /dev/null @@ -1,38 +0,0 @@ -# Normally, finding a boost component should rely on -# find_package(Boost REQUIRED COMPONENTS python) -# But the packaging is not always the same: sometimes it will be python3, -# others python37, python2, python27, or just python -# (i.e. python3 for Fedora < 30, but python37 for Fedora >= 30 and MacOSX via Homebrew) -# We wrap all this in this module -# -# Defines: -# BoostPython_FOUND -# BoostPython_INCLUDE_DIRS -# BoostPython_LIBRARIES - -if (NOT BoostPython_FOUND) - find_package(PythonInterp ${PYTHON_EXPLICIT_VERSION} REQUIRED) - - # Explicit versions - set (_BOOST_PYTHON_LIST "python${PYTHON_VERSION_MAJOR}${PYTHON_VERSION_MINOR}" "python${PYTHON_VERSION_MAJOR}") - - # Only for python2, without suffix - if (BoostPython_FIND_VERSION VERSION_LESS 3) - list (APPEND _BOOST_PYTHON_LIST "python") - endif () - - # Pick the most restrictive - foreach (_BOOST_PYTHON IN LISTS _BOOST_PYTHON_LIST) - find_package(Boost COMPONENTS "${_BOOST_PYTHON}") - if (Boost_FOUND) - set (BoostPython_INCLUDE_DIRS ${Boost_INCLUDE_DIRS}) - set (BoostPython_LIBRARIES ${Boost_LIBRARIES}) - set (BoostPython_FOUND TRUE) - break () - endif () - endforeach () - - include(FindPackageHandleStandardArgs) - find_package_handle_standard_args(BoostPython DEFAULT_MSG BoostPython_INCLUDE_DIRS) - mark_as_advanced(BoostPython_FOUND BoostPython_INCLUDE_DIRS BoostPython_LIBRARIES) -endif () From 4955ca140d428154909e8662b6627bbe2d4f2ad6 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Fri, 17 Apr 2026 14:56:11 +0200 Subject: [PATCH 27/40] changing MeasurementFrameRectangle.h --- .../MeasurementFrameRectangle.h | 25 ++++++------------- .../FlexibleModelFittingIterativeTask.cpp | 15 ++++++----- 2 files changed, 15 insertions(+), 25 deletions(-) diff --git a/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h b/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h index dbd27895..586b1f0b 100644 --- a/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h +++ b/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h @@ -39,54 +39,45 @@ class MeasurementFrameRectangle: public Property { virtual ~MeasurementFrameRectangle() = default; explicit MeasurementFrameRectangle(bool bad_projection): - m_min_coord{-1, -1}, m_max_coord{-1, -1}, m_bad_projection{bad_projection} {} + m_min_coord{-1, -1}, m_max_coord{-1, -1}, m_bad_projection{bad_projection}, m_is_valid{false} {} MeasurementFrameRectangle(ImageCoordinate min_coord, ImageCoordinate max_coord): - m_min_coord(min_coord), m_max_coord(max_coord), m_bad_projection{false} { + m_min_coord(min_coord), m_max_coord(max_coord), m_bad_projection{false}, m_is_valid{true} { assert(min_coord.m_x <= max_coord.m_x && min_coord.m_y <= max_coord.m_y); } ImageCoordinate getTopLeft() const { - assert(m_min_coord.m_x >= 0.0); - assert(m_min_coord.m_y >= 0.0); return m_min_coord; } ImageCoordinate getBottomRight() const { - assert(m_max_coord.m_x >= 0.0); - assert(m_max_coord.m_y >= 0.0); return m_max_coord; } double getWidth() const { - if (m_max_coord.m_x <= 0.0) { - return 0.0; - } return m_max_coord.m_x - m_min_coord.m_x; } double getHeight() const { - if (m_max_coord.m_y <= 0.0) { - return 0.0; - } return m_max_coord.m_y - m_min_coord.m_y; } - std::tuple getImageRect() const { - return std::make_tuple(m_min_coord, m_max_coord); - } - bool isBadProjection() const { return m_bad_projection; } bool isValid() const { - return !m_bad_projection && m_max_coord.m_x > 0.0 && m_max_coord.m_y > 0.0 && getWidth() > 0.0 && getHeight() > 0.0; + return m_is_valid; + } + + bool isEmpty() const { + return getWidth() <= 0.0 || getHeight() <= 0.0; } private: ImageCoordinate m_min_coord, m_max_coord; bool m_bad_projection; + bool m_is_valid; }; } // end SourceXtractor diff --git a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp index 755834d8..881ba79e 100644 --- a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp +++ b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp @@ -74,14 +74,13 @@ FlexibleModelFittingIterativeTask::~FlexibleModelFittingIterativeTask() { } PixelRectangle FlexibleModelFittingIterativeTask::getUnclippedFittingRect(SourceInterface& source, int frame_index) const { - ImageCoordinate min_coord, max_coord; auto& measurement_frame_rectangle = source.getProperty(frame_index); - std::tie(min_coord, max_coord) = measurement_frame_rectangle.getImageRect(); - - if (!measurement_frame_rectangle.isValid()) { + if (!measurement_frame_rectangle.isValid() || measurement_frame_rectangle.isEmpty()) { return PixelRectangle(); } - + + ImageCoordinate min_coord = measurement_frame_rectangle.getTopLeft(); + ImageCoordinate max_coord = measurement_frame_rectangle.getBottomRight(); if (m_window_type == WindowType::ROTATED_ELLIPSE) { auto ellipse = getFittingEllipse(source, frame_index); @@ -91,7 +90,7 @@ PixelRectangle FlexibleModelFittingIterativeTask::getUnclippedFittingRect(Source return getEllipseRect(ellipse); } - if (max_coord.m_x - min_coord.m_x <= 0 || max_coord.m_y - min_coord.m_y <= 0) { + if ((max_coord.m_x - min_coord.m_x <= 0.0) || (max_coord.m_y - min_coord.m_y <= 0.0)) { return PixelRectangle(); } else { auto min = min_coord; @@ -132,8 +131,8 @@ PixelRectangle FlexibleModelFittingIterativeTask::getUnclippedFittingRect(Source max.m_y = min.m_y + size; } - auto min_pc = PixelCoordinate(static_cast(min.m_x + 0.5), static_cast(min.m_y + 0.5)); - auto max_pc = PixelCoordinate(static_cast(max.m_x + 0.5), static_cast(max.m_y + 0.5)); + auto min_pc = PixelCoordinate(static_cast(min.m_x), static_cast(min.m_y)); + auto max_pc = PixelCoordinate(static_cast(max.m_x + 1.0), static_cast(max.m_y + 1.0)); return PixelRectangle(min_pc, max_pc); } From 3061074a51602295ed344524ac580618660cd500 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Fri, 17 Apr 2026 16:35:36 +0200 Subject: [PATCH 28/40] fix litmus test failing --- .../FlexibleModelFittingIterativeTask.cpp | 2 +- .../MeasurementFrameRectangleTask.cpp | 31 +++++++++++-------- .../MeasurementFrameRectangleTaskNoDetect.cpp | 31 +++++++++++-------- 3 files changed, 37 insertions(+), 27 deletions(-) diff --git a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp index 881ba79e..378f8b40 100644 --- a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp +++ b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp @@ -132,7 +132,7 @@ PixelRectangle FlexibleModelFittingIterativeTask::getUnclippedFittingRect(Source } auto min_pc = PixelCoordinate(static_cast(min.m_x), static_cast(min.m_y)); - auto max_pc = PixelCoordinate(static_cast(max.m_x + 1.0), static_cast(max.m_y + 1.0)); + auto max_pc = PixelCoordinate(static_cast(max.m_x), static_cast(max.m_y)); return PixelRectangle(min_pc, max_pc); } diff --git a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp index 1a635308..6fd6bc44 100644 --- a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp +++ b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp @@ -52,7 +52,6 @@ void MeasurementFrameRectangleTask::computeProperties(SourceInterface& source) c // Transform the 4 corner coordinates from detection image ImageCoordinate coord1, coord2, coord3, coord4; - bool bad_projection = false; try { coord1 = measurement_frame_coordinates->worldToImage( detection_frame_coordinates->imageToWorld(ImageCoordinate(stamp_top_left.m_x, stamp_top_left.m_y))); @@ -64,7 +63,8 @@ void MeasurementFrameRectangleTask::computeProperties(SourceInterface& source) c detection_frame_coordinates->imageToWorld(ImageCoordinate(stamp_top_left.m_x, stamp_top_left.m_y + height))); } catch (const InvalidCoordinatesException&) { - bad_projection = true; + source.setIndexedProperty(m_instance, true); + return; } // Determine the min/max coordinates @@ -73,21 +73,26 @@ void MeasurementFrameRectangleTask::computeProperties(SourceInterface& source) c auto max_x = std::max(coord1.m_x, std::max(coord2.m_x, std::max(coord3.m_x, coord4.m_x))); auto max_y = std::max(coord1.m_y, std::max(coord2.m_y, std::max(coord3.m_y, coord4.m_y))); - // The full boundaries may lie outside of the frame - if (bad_projection || max_x < 0.0 || max_y < 0.0 || - int(min_x) >= measurement_frame_info.getWidth() || int(max_y) >= measurement_frame_info.getHeight()) { - source.setIndexedProperty(m_instance, bad_projection); + // test for nans + if (std::isnan(min_x) || std::isnan(min_y) || std::isnan(max_x) || std::isnan(max_y)) { + source.setIndexedProperty(m_instance, true); + return; } + // Clip the coordinates to fit the available image - else { - min_x = std::max(0.0, min_x); - min_y = std::max(0.0, min_y); - max_x = std::min(static_cast(measurement_frame_info.getWidth()), max_x); - max_y = std::min(static_cast(measurement_frame_info.getHeight()), max_y); + min_x = std::max(0.0, min_x); + min_y = std::max(0.0, min_y); + max_x = std::min(static_cast(measurement_frame_info.getWidth() - 1), max_x); + max_y = std::min(static_cast(measurement_frame_info.getHeight() - 1), max_y); - source.setIndexedProperty( - m_instance, ImageCoordinate(min_x, min_y), ImageCoordinate(max_x, max_y)); + if (max_x <= min_x || max_y <= min_y) { + // No valid rectangle, set the property to an empty rectangle + source.setIndexedProperty(m_instance, false); + return; } + + source.setIndexedProperty( + m_instance, ImageCoordinate(min_x, min_y), ImageCoordinate(max_x, max_y)); } } // SEImplementation namespace diff --git a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp index f5757b05..b6fcbf29 100644 --- a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp +++ b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp @@ -38,7 +38,6 @@ void MeasurementFrameRectangleTaskNoDetect::computeProperties(SourceInterface& s const auto& world_centroid = source.getProperty(); const auto& assoc_mode = source.getProperty(); - bool bad_projection = false; ImageCoordinate coord1, coord2, coord3, coord4; try { auto w = assoc_mode.getRefFramePixelWidth(); @@ -55,7 +54,8 @@ void MeasurementFrameRectangleTaskNoDetect::computeProperties(SourceInterface& s reference_frame_coordinates->imageToWorld(ImageCoordinate(c.m_x + w, c.m_y + h))); } catch (const InvalidCoordinatesException&) { - bad_projection = true; + source.setIndexedProperty(m_instance, true); + return; } // Determine the min/max coordinates @@ -64,21 +64,26 @@ void MeasurementFrameRectangleTaskNoDetect::computeProperties(SourceInterface& s auto max_x = std::max(coord1.m_x, std::max(coord2.m_x, std::max(coord3.m_x, coord4.m_x))); auto max_y = std::max(coord1.m_y, std::max(coord2.m_y, std::max(coord3.m_y, coord4.m_y))); - // The full boundaries may lie outside of the frame - if (bad_projection || max_x < 0.0 || max_y < 0.0 || - int(min_x) >= measurement_frame_info.getWidth() || int(max_y) >= measurement_frame_info.getHeight()) { - source.setIndexedProperty(m_instance, bad_projection); + // test for nans + if (std::isnan(min_x) || std::isnan(min_y) || std::isnan(max_x) || std::isnan(max_y)) { + source.setIndexedProperty(m_instance, true); + return; } + // Clip the coordinates to fit the available image - else { - min_x = std::max(0.0, min_x); - min_y = std::max(0.0, min_y); - max_x = std::min(static_cast(measurement_frame_info.getWidth()), max_x); - max_y = std::min(static_cast(measurement_frame_info.getHeight()), max_y); + min_x = std::max(0.0, min_x); + min_y = std::max(0.0, min_y); + max_x = std::min(static_cast(measurement_frame_info.getWidth() - 1), max_x); + max_y = std::min(static_cast(measurement_frame_info.getHeight() - 1), max_y); - source.setIndexedProperty( - m_instance, ImageCoordinate(min_x, min_y), ImageCoordinate(max_x, max_y)); + if (max_x <= min_x || max_y <= min_y) { + // No valid rectangle, set the property to an empty rectangle + source.setIndexedProperty(m_instance, false); + return; } + + source.setIndexedProperty( + m_instance, ImageCoordinate(min_x, min_y), ImageCoordinate(max_x, max_y)); } } // SEImplementation namespace From a7eeab80224849e764e4134dd96ab1afefdc30c8 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Thu, 23 Apr 2026 12:53:59 +0200 Subject: [PATCH 29/40] add framevalid check when building deblend image --- .../FlexibleModelFittingIterativeTask.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp index 378f8b40..71fc5859 100644 --- a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp +++ b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp @@ -532,6 +532,11 @@ std::shared_ptr> FlexibleModelFittingIterativeTask::createD index = 0; for (auto& src : group) { if (index != source_index) { + // check if frame is valid + if (!isFrameValid(src, frame->getFrameNb())) { + index++; + continue; + } auto frame_model = createFrameModel(src, pixel_scale, parameter_manager, frame, rect); auto final_stamp = frame_model.getImage(); From 8149c6dfe3170995e0a88aaddcf7c4373fdcf1b2 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Thu, 23 Apr 2026 13:33:55 +0200 Subject: [PATCH 30/40] more infinite tests --- .../FlexibleModelFittingIterativeTask.cpp | 19 +++++++------------ .../MeasurementFrameRectangleTask.cpp | 6 ++++++ .../MeasurementFrameRectangleTaskNoDetect.cpp | 6 ++++++ 3 files changed, 19 insertions(+), 12 deletions(-) diff --git a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp index 71fc5859..278af1b1 100644 --- a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp +++ b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp @@ -531,20 +531,15 @@ std::shared_ptr> FlexibleModelFittingIterativeTask::createD auto deblend_image = VectorImage::create(rect.getWidth(), rect.getHeight()); index = 0; for (auto& src : group) { - if (index != source_index) { - // check if frame is valid - if (!isFrameValid(src, frame->getFrameNb())) { - index++; - continue; - } - auto frame_model = createFrameModel(src, pixel_scale, parameter_manager, frame, rect); - auto final_stamp = frame_model.getImage(); + if (index != source_index && isFrameValid(src, frame->getFrameNb())) { + auto frame_model = createFrameModel(src, pixel_scale, parameter_manager, frame, rect); + auto final_stamp = frame_model.getImage(); - for (int y = 0; y < final_stamp->getHeight(); ++y) { - for (int x = 0; x < final_stamp->getWidth(); ++x) { - deblend_image->at(x, y) += final_stamp->at(x, y); - } + for (int y = 0; y < final_stamp->getHeight(); ++y) { + for (int x = 0; x < final_stamp->getWidth(); ++x) { + deblend_image->at(x, y) += final_stamp->at(x, y); } + } } index++; } diff --git a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp index 6fd6bc44..7594eee6 100644 --- a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp +++ b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp @@ -79,6 +79,12 @@ void MeasurementFrameRectangleTask::computeProperties(SourceInterface& source) c return; } + // test for infinite values + if (std::isinf(min_x) || std::isinf(min_y) || std::isinf(max_x) || std::isinf(max_y)) { + source.setIndexedProperty(m_instance, true); + return; + } + // Clip the coordinates to fit the available image min_x = std::max(0.0, min_x); min_y = std::max(0.0, min_y); diff --git a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp index b6fcbf29..3d7778bf 100644 --- a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp +++ b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp @@ -70,6 +70,12 @@ void MeasurementFrameRectangleTaskNoDetect::computeProperties(SourceInterface& s return; } + // test for infinite values + if (std::isinf(min_x) || std::isinf(min_y) || std::isinf(max_x) || std::isinf(max_y)) { + source.setIndexedProperty(m_instance, true); + return; + } + // Clip the coordinates to fit the available image min_x = std::max(0.0, min_x); min_y = std::max(0.0, min_y); From 2ee3de39038b29d0ca70543072a1e0292089db7e Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Thu, 23 Apr 2026 14:07:41 +0200 Subject: [PATCH 31/40] remove infiite in case of error --- SEFramework/src/lib/CoordinateSystem/WCS.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/SEFramework/src/lib/CoordinateSystem/WCS.cpp b/SEFramework/src/lib/CoordinateSystem/WCS.cpp index ac14c061..e6ff5142 100644 --- a/SEFramework/src/lib/CoordinateSystem/WCS.cpp +++ b/SEFramework/src/lib/CoordinateSystem/WCS.cpp @@ -255,8 +255,10 @@ ImageCoordinate WCS::worldToImage(WorldCoordinate world_coordinate) const { int ret_val = wcss2p(&wcs_copy, 1, 1, wc_array, &phi, &theta, ic_array, pc_array, &status); if (ret_val != WCSERR_SUCCESS) { logger.warn() << "Bad worldToImage from RA/Dec: " << wc_array[0] << "/" << wc_array[1]; - pc_array[0] = -std::numeric_limits::infinity(); - pc_array[1] = -std::numeric_limits::infinity(); + // pc_array[0] = -std::numeric_limits::infinity(); + // pc_array[1] = -std::numeric_limits::infinity(); + pc_array[0] = -10000000.0; + pc_array[1] = -10000000.0; } wcsfree(&wcs_copy); return ImageCoordinate(pc_array[0] - 1, pc_array[1] - 1); // -1 as fits standard coordinates start at 1 From e95897cb9f6ea4071e2632e41d192f0e9d3e99df Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Fri, 24 Apr 2026 13:43:35 +0200 Subject: [PATCH 32/40] remove commented out code --- SEFramework/src/lib/CoordinateSystem/WCS.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/SEFramework/src/lib/CoordinateSystem/WCS.cpp b/SEFramework/src/lib/CoordinateSystem/WCS.cpp index e6ff5142..b3c19d36 100644 --- a/SEFramework/src/lib/CoordinateSystem/WCS.cpp +++ b/SEFramework/src/lib/CoordinateSystem/WCS.cpp @@ -255,8 +255,6 @@ ImageCoordinate WCS::worldToImage(WorldCoordinate world_coordinate) const { int ret_val = wcss2p(&wcs_copy, 1, 1, wc_array, &phi, &theta, ic_array, pc_array, &status); if (ret_val != WCSERR_SUCCESS) { logger.warn() << "Bad worldToImage from RA/Dec: " << wc_array[0] << "/" << wc_array[1]; - // pc_array[0] = -std::numeric_limits::infinity(); - // pc_array[1] = -std::numeric_limits::infinity(); pc_array[0] = -10000000.0; pc_array[1] = -10000000.0; } From 984c21b2167105e64359524ce152df9271ca4708 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Thu, 30 Apr 2026 16:25:42 +0200 Subject: [PATCH 33/40] Reapply "Changed fitting window calculation to use doubles" This reverts commit bcb72281967f6a4acc60bf359d7291718aaf92d7. --- .../CoordinateSystem/CoordinateSystem.h | 24 +++++++++++++++++ .../MeasurementFrameRectangle.h | 18 +++++++++---- .../FlexibleModelFittingIterativeTask.cpp | 16 +++++++----- .../MeasurementFrameRectangleTask.cpp | 3 ++- .../MeasurementFrameRectangleTaskNoDetect.cpp | 26 ++++++++++++------- 5 files changed, 65 insertions(+), 22 deletions(-) diff --git a/SEFramework/SEFramework/CoordinateSystem/CoordinateSystem.h b/SEFramework/SEFramework/CoordinateSystem/CoordinateSystem.h index d43e26a4..30605da8 100644 --- a/SEFramework/SEFramework/CoordinateSystem/CoordinateSystem.h +++ b/SEFramework/SEFramework/CoordinateSystem/CoordinateSystem.h @@ -45,6 +45,30 @@ struct ImageCoordinate { ImageCoordinate() : m_x(0), m_y(0) {} ImageCoordinate(double x, double y) : m_x(x), m_y(y) {} + + ImageCoordinate operator*(double scalar) const { + return ImageCoordinate(m_x * scalar, m_y * scalar); + } + + ImageCoordinate operator+(const ImageCoordinate& other) const { + return ImageCoordinate(m_x + other.m_x, m_y + other.m_y); + } + + ImageCoordinate& operator+=(const ImageCoordinate& other) { + m_x += other.m_x; + m_y += other.m_y; + return *this; + } + + ImageCoordinate operator-(const ImageCoordinate& other) const { + return ImageCoordinate(m_x - other.m_x, m_y - other.m_y); + } + + ImageCoordinate& operator-=(const ImageCoordinate& other) { + m_x -= other.m_x; + m_y -= other.m_y; + return *this; + } }; class CoordinateSystem { diff --git a/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h b/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h index fb7ae567..d163e632 100644 --- a/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h +++ b/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h @@ -24,7 +24,10 @@ #ifndef _SEIMPLEMENTATION_PLUGIN_MEASUREMENTFRAMERECTANGLE_MEASUREMENTFRAMERECTANGLE_H_ #define _SEIMPLEMENTATION_PLUGIN_MEASUREMENTFRAMERECTANGLE_MEASUREMENTFRAMERECTANGLE_H_ +#include + #include "SEUtils/PixelRectangle.h" +#include "SEFramework/CoordinateSystem/CoordinateSystem.h" #include "SEFramework/Property/Property.h" #include "SEFramework/Image/Image.h" @@ -36,11 +39,15 @@ class MeasurementFrameRectangle: public Property { virtual ~MeasurementFrameRectangle() = default; explicit MeasurementFrameRectangle(bool bad_projection): - m_min_coord{-1, -1}, m_max_coord{-1, -1}, m_bad_projection{bad_projection}{} + m_min_coord{-1, -1}, m_max_coord{-1, -1}, m_min_coord_image{-1, -1}, + m_max_coord_image{-1, -1}, m_bad_projection{bad_projection} {} - MeasurementFrameRectangle(PixelCoordinate min_coord, PixelCoordinate max_coord): - m_min_coord{min_coord}, m_max_coord{max_coord}, m_bad_projection{false} { + MeasurementFrameRectangle(PixelCoordinate min_coord, PixelCoordinate max_coord, + ImageCoordinate min_coord_image, ImageCoordinate max_coord_image): + m_min_coord(min_coord), m_max_coord(max_coord), + m_min_coord_image(min_coord_image), m_max_coord_image(max_coord_image), m_bad_projection{false} { assert(min_coord.m_x <= max_coord.m_x && min_coord.m_y <= max_coord.m_y); + assert(min_coord_image.m_x <= max_coord_image.m_x && min_coord_image.m_y <= max_coord_image.m_y); } PixelCoordinate getTopLeft() const { @@ -65,8 +72,8 @@ class MeasurementFrameRectangle: public Property { return m_max_coord.m_y - m_min_coord.m_y + 1; } - PixelRectangle getRect() const { - return PixelRectangle(m_min_coord, m_max_coord); + std::tuple getImageRect() const { + return std::make_tuple(m_min_coord_image, m_max_coord_image); } bool badProjection() const { @@ -75,6 +82,7 @@ class MeasurementFrameRectangle: public Property { private: PixelCoordinate m_min_coord, m_max_coord; + ImageCoordinate m_min_coord_image, m_max_coord_image; bool m_bad_projection; }; diff --git a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp index 4d7e91bd..5a669b5b 100644 --- a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp +++ b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp @@ -74,7 +74,8 @@ FlexibleModelFittingIterativeTask::~FlexibleModelFittingIterativeTask() { } PixelRectangle FlexibleModelFittingIterativeTask::getUnclippedFittingRect(SourceInterface& source, int frame_index) const { - auto fitting_rect = source.getProperty(frame_index).getRect(); + ImageCoordinate min_coord, max_coord; + std::tie(min_coord, max_coord) = source.getProperty(frame_index).getImageRect(); if (m_window_type == WindowType::ROTATED_ELLIPSE) { auto ellipse = getFittingEllipse(source, frame_index); @@ -84,14 +85,14 @@ PixelRectangle FlexibleModelFittingIterativeTask::getUnclippedFittingRect(Source return getEllipseRect(ellipse); } - if (fitting_rect.getWidth() <= 0 || fitting_rect.getHeight() <= 0) { + if (max_coord.m_x - min_coord.m_x <= 0 || max_coord.m_y - min_coord.m_y <= 0) { return PixelRectangle(); } else { - auto min = fitting_rect.getTopLeft(); - auto max = fitting_rect.getBottomRight(); + auto min = min_coord; + auto max = max_coord; // FIXME temporary, for now just enlarge the area by a fixed amount of pixels - PixelCoordinate border = (max - min) * .8 + PixelCoordinate(2, 2); + ImageCoordinate border = (max - min) * .8 + ImageCoordinate(2.0, 2.0); min -= border; max += border; @@ -125,7 +126,10 @@ PixelRectangle FlexibleModelFittingIterativeTask::getUnclippedFittingRect(Source max.m_y = min.m_y + size; } - return PixelRectangle(min, max); + auto min_pc = PixelCoordinate(static_cast(min.m_x + 0.5), static_cast(min.m_y + 0.5)); + auto max_pc = PixelCoordinate(static_cast(max.m_x + 0.5), static_cast(max.m_y + 0.5)); + + return PixelRectangle(min_pc, max_pc); } } diff --git a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp index f57d56fb..da88a8dc 100644 --- a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp +++ b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp @@ -85,7 +85,8 @@ void MeasurementFrameRectangleTask::computeProperties(SourceInterface& source) c max_coord.m_x = std::min(measurement_frame_info.getWidth() - 1, max_coord.m_x); max_coord.m_y = std::min(measurement_frame_info.getHeight() - 1, max_coord.m_y); - source.setIndexedProperty(m_instance, min_coord, max_coord); + source.setIndexedProperty( + m_instance, min_coord, max_coord, ImageCoordinate(min_x, min_y), ImageCoordinate(max_x, max_y)); } } diff --git a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp index 68caedfc..2cf63a2c 100644 --- a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp +++ b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTaskNoDetect.cpp @@ -23,6 +23,7 @@ #include "SEImplementation/Plugin/DetectionFrameCoordinates/DetectionFrameCoordinates.h" #include #include +#include #include #include @@ -31,23 +32,27 @@ namespace SourceXtractor { void MeasurementFrameRectangleTaskNoDetect::computeProperties(SourceInterface& source) const { auto measurement_frame_coordinates = source.getProperty(m_instance).getCoordinateSystem(); + auto reference_frame_coordinates = source.getProperty().getCoordinateSystem(); + const auto& measurement_frame_info = source.getProperty(m_instance); const auto& world_centroid = source.getProperty(); const auto& assoc_mode = source.getProperty(); - auto coord = world_centroid.getCentroid(); - bool bad_coordinates = false; ImageCoordinate coord1, coord2, coord3, coord4; try { - int w = assoc_mode.getRefFramePixelWidth(); - int h = assoc_mode.getRefFramePixelHeight(); + auto w = assoc_mode.getRefFramePixelWidth(); + auto h = assoc_mode.getRefFramePixelHeight(); - auto c = measurement_frame_coordinates->worldToImage(coord); - coord1 = ImageCoordinate(c.m_x - w, c.m_y - h); - coord2 = ImageCoordinate(c.m_x + w, c.m_y - h); - coord3 = ImageCoordinate(c.m_x - w, c.m_y + h); - coord4 = ImageCoordinate(c.m_x + w, c.m_y + h); + auto c = reference_frame_coordinates->worldToImage(world_centroid.getCentroid()); + coord1 = measurement_frame_coordinates->worldToImage( + reference_frame_coordinates->imageToWorld(ImageCoordinate(c.m_x - w, c.m_y - h))); + coord2 = measurement_frame_coordinates->worldToImage( + reference_frame_coordinates->imageToWorld(ImageCoordinate(c.m_x + w, c.m_y - h))); + coord3 = measurement_frame_coordinates->worldToImage( + reference_frame_coordinates->imageToWorld(ImageCoordinate(c.m_x - w, c.m_y + h))); + coord4 = measurement_frame_coordinates->worldToImage( + reference_frame_coordinates->imageToWorld(ImageCoordinate(c.m_x + w, c.m_y + h))); } catch (const InvalidCoordinatesException&) { bad_coordinates = true; @@ -77,7 +82,8 @@ void MeasurementFrameRectangleTaskNoDetect::computeProperties(SourceInterface& s max_coord.m_x = std::min(measurement_frame_info.getWidth() - 1, max_coord.m_x); max_coord.m_y = std::min(measurement_frame_info.getHeight() - 1, max_coord.m_y); - source.setIndexedProperty(m_instance, min_coord, max_coord); + source.setIndexedProperty( + m_instance, min_coord, max_coord, ImageCoordinate(min_x, min_y), ImageCoordinate(max_x, max_y)); } } From c2b1d44d21f85afa21489446a4966ab1094f69a0 Mon Sep 17 00:00:00 2001 From: Marc Schefer <31961290+marcschefer@users.noreply.github.com> Date: Wed, 3 Jun 2026 15:03:37 +0200 Subject: [PATCH 34/40] Update CMakeLists.txt Bump Alexandria version --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 265b89ed..73d2679b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -6,4 +6,4 @@ find_package(ElementsProject) #--------------------------------------------------------------- # Declare project name and version -elements_project(SourceXtractorPlusPlus 1.1.0 USE Alexandria 2.32.0 DESCRIPTION "SourceXtractor++, the next generation SExtractor") +elements_project(SourceXtractorPlusPlus 1.1.0 USE Alexandria 2.32.2 DESCRIPTION "SourceXtractor++, the next generation SExtractor") From 0df24e56b67045a54aa8131ac6b7f9c0b20ec889 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Fri, 5 Jun 2026 11:45:47 +0200 Subject: [PATCH 35/40] bump cmake version --- CMakeLists.txt | 2 +- SEBenchmarks/CMakeLists.txt | 2 -- SEFramework/CMakeLists.txt | 2 -- SEImplementation/CMakeLists.txt | 2 -- SEMain/CMakeLists.txt | 2 -- SEUtils/CMakeLists.txt | 2 -- 6 files changed, 1 insertion(+), 11 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 265b89ed..29e20607 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,4 +1,4 @@ -CMAKE_MINIMUM_REQUIRED(VERSION 2.8.5) +CMAKE_MINIMUM_REQUIRED(VERSION 3.5) #--------------------------------------------------------------- # Load macros and functions for Elements-based projects diff --git a/SEBenchmarks/CMakeLists.txt b/SEBenchmarks/CMakeLists.txt index 0b5b7dd3..d779a967 100644 --- a/SEBenchmarks/CMakeLists.txt +++ b/SEBenchmarks/CMakeLists.txt @@ -1,5 +1,3 @@ -CMAKE_MINIMUM_REQUIRED(VERSION 2.8.5) - #=============================================================================== # Load elements_subdir macro here # Examples: diff --git a/SEFramework/CMakeLists.txt b/SEFramework/CMakeLists.txt index 77555254..fcf4a327 100644 --- a/SEFramework/CMakeLists.txt +++ b/SEFramework/CMakeLists.txt @@ -1,5 +1,3 @@ -CMAKE_MINIMUM_REQUIRED(VERSION 2.8.5) - #=============================================================================== # Load elements_subdir macro here # Examples: diff --git a/SEImplementation/CMakeLists.txt b/SEImplementation/CMakeLists.txt index 1b52a93d..75b0df5e 100644 --- a/SEImplementation/CMakeLists.txt +++ b/SEImplementation/CMakeLists.txt @@ -1,5 +1,3 @@ -CMAKE_MINIMUM_REQUIRED(VERSION 2.8.5) - #=============================================================================== # Load elements_subdir macro here # Examples: diff --git a/SEMain/CMakeLists.txt b/SEMain/CMakeLists.txt index 8fdb92a8..b904bc4f 100644 --- a/SEMain/CMakeLists.txt +++ b/SEMain/CMakeLists.txt @@ -1,5 +1,3 @@ -CMAKE_MINIMUM_REQUIRED(VERSION 2.8.5) - #=============================================================================== # Load elements_subdir macro here # Examples: diff --git a/SEUtils/CMakeLists.txt b/SEUtils/CMakeLists.txt index 64bcc685..8203f66f 100644 --- a/SEUtils/CMakeLists.txt +++ b/SEUtils/CMakeLists.txt @@ -1,5 +1,3 @@ -CMAKE_MINIMUM_REQUIRED(VERSION 2.8.5) - #=============================================================================== # Load elements_subdir macro here # Examples: From 4896489430adfd944ef390702f6ca54f99e71438 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Fri, 5 Jun 2026 17:16:17 +0200 Subject: [PATCH 36/40] add pearson factor parameter --- .../Background/BackgroundAnalyzerFactory.h | 2 ++ .../SEImplementation/Background/SE/ImageMode.h | 4 ++-- .../Background/SE/SEBackgroundLevelAnalyzer.h | 3 ++- .../SEImplementation/Configuration/BackgroundConfig.h | 6 ++++++ .../src/lib/Background/BackgroundAnalyzerFactory.cpp | 7 +++++-- SEImplementation/src/lib/Background/SE/ImageMode.cpp | 6 +++--- .../src/lib/Background/SE/SEBackgroundLevelAnalyzer.cpp | 7 ++++--- .../src/lib/Configuration/BackgroundConfig.cpp | 9 ++++++++- 8 files changed, 32 insertions(+), 12 deletions(-) diff --git a/SEImplementation/SEImplementation/Background/BackgroundAnalyzerFactory.h b/SEImplementation/SEImplementation/Background/BackgroundAnalyzerFactory.h index 96883fcc..ca20b183 100644 --- a/SEImplementation/SEImplementation/Background/BackgroundAnalyzerFactory.h +++ b/SEImplementation/SEImplementation/Background/BackgroundAnalyzerFactory.h @@ -26,6 +26,7 @@ #include "Configuration/Configuration.h" +#include "SEImplementation/Configuration/BackgroundConfig.h" #include "SEImplementation/Configuration/SE2BackgroundConfig.h" #include "SEImplementation/Configuration/WeightImageConfig.h" @@ -49,6 +50,7 @@ class BackgroundAnalyzerFactory : public Euclid::Configuration::Configuration { std::vector m_cell_size; std::vector m_smoothing_box; WeightImageConfig::WeightType m_weight_type; + SeFloat m_pearson_factor; }; } diff --git a/SEImplementation/SEImplementation/Background/SE/ImageMode.h b/SEImplementation/SEImplementation/Background/SE/ImageMode.h index 57f9c045..f4f2ddaf 100644 --- a/SEImplementation/SEImplementation/Background/SE/ImageMode.h +++ b/SEImplementation/SEImplementation/Background/SE/ImageMode.h @@ -68,7 +68,7 @@ class ImageMode { ImageMode(const std::shared_ptr>& image, const std::shared_ptr>& variance, int cell_w, int cell_h, T invalid_value, T kappa1 = 2, T kappa2 = 5, T kappa3 = 3, - T rtol = 1e-4, size_t max_iter = 100); + T rtol = 1e-4, size_t max_iter = 100, T pearson_factor = 2.5); /** * Destructor @@ -106,7 +106,7 @@ class ImageMode { std::shared_ptr> m_mode, m_sigma; std::shared_ptr> m_var_mode, m_var_sigma; int m_cell_w, m_cell_h; - T m_invalid, m_kappa1, m_kappa2, m_kappa3, m_rtol; + T m_invalid, m_kappa1, m_kappa2, m_kappa3, m_rtol, m_pearson_factor; size_t m_max_iter; std::tuple getBackGuess(const std::vector &data) const; diff --git a/SEImplementation/SEImplementation/Background/SE/SEBackgroundLevelAnalyzer.h b/SEImplementation/SEImplementation/Background/SE/SEBackgroundLevelAnalyzer.h index 1fe32982..1e811743 100644 --- a/SEImplementation/SEImplementation/Background/SE/SEBackgroundLevelAnalyzer.h +++ b/SEImplementation/SEImplementation/Background/SE/SEBackgroundLevelAnalyzer.h @@ -28,7 +28,7 @@ namespace SourceXtractor { class SEBackgroundLevelAnalyzer : public BackgroundAnalyzer { public: SEBackgroundLevelAnalyzer(const std::vector& cell_size, const std::vector& smoothing_box, - const WeightImageConfig::WeightType weight_type); + const WeightImageConfig::WeightType weight_type, SeFloat pearson_factor = 2.5); virtual ~SEBackgroundLevelAnalyzer() = default; @@ -41,6 +41,7 @@ class SEBackgroundLevelAnalyzer : public BackgroundAnalyzer { std::array m_smoothing_box; WeightImageConfig::WeightType m_weight_type; + SeFloat m_pearson_factor; }; } // end of namespace SourceXtractor diff --git a/SEImplementation/SEImplementation/Configuration/BackgroundConfig.h b/SEImplementation/SEImplementation/Configuration/BackgroundConfig.h index 2f1445e5..772b9495 100644 --- a/SEImplementation/SEImplementation/Configuration/BackgroundConfig.h +++ b/SEImplementation/SEImplementation/Configuration/BackgroundConfig.h @@ -55,12 +55,18 @@ class BackgroundConfig : public Euclid::Configuration::Configuration { return m_detection_threshold_absolute; } + SeFloat getPearsonFactor() const { + return m_pearson_factor; + } + private: bool m_background_level_absolute; SeFloat m_background_level; bool m_detection_threshold_absolute; SeFloat m_detection_threshold; + + SeFloat m_pearson_factor; }; } /* namespace SourceXtractor */ diff --git a/SEImplementation/src/lib/Background/BackgroundAnalyzerFactory.cpp b/SEImplementation/src/lib/Background/BackgroundAnalyzerFactory.cpp index 3c5d31ff..f7a19e19 100644 --- a/SEImplementation/src/lib/Background/BackgroundAnalyzerFactory.cpp +++ b/SEImplementation/src/lib/Background/BackgroundAnalyzerFactory.cpp @@ -37,7 +37,7 @@ std::shared_ptr BackgroundAnalyzerFactory::createBackgroundA WeightImageConfig::WeightType weight_type) const { // make a SE2 background if cell size and smoothing box are given if (m_cell_size.size() > 0 && m_smoothing_box.size() > 0) { - return std::make_shared(m_cell_size, m_smoothing_box, weight_type); + return std::make_shared(m_cell_size, m_smoothing_box, weight_type, m_pearson_factor); } else { // make a simple background return std::make_shared(); @@ -45,17 +45,20 @@ std::shared_ptr BackgroundAnalyzerFactory::createBackgroundA } BackgroundAnalyzerFactory::BackgroundAnalyzerFactory(long manager_id) - : Configuration(manager_id), m_weight_type(WeightImageConfig::WeightType::WEIGHT_TYPE_NONE) { + : Configuration(manager_id), m_weight_type(WeightImageConfig::WeightType::WEIGHT_TYPE_NONE), m_pearson_factor(2.5) { declareDependency(); + declareDependency(); declareDependency(); } void BackgroundAnalyzerFactory::initialize(const UserValues&) { + auto background_config = getDependency(); auto se2background_config = getDependency(); auto weight_image_config = getDependency(); m_cell_size = se2background_config.getCellSize(); m_smoothing_box = se2background_config.getSmoothingBox(); m_weight_type = weight_image_config.getWeightType(); + m_pearson_factor = background_config.getPearsonFactor(); } } diff --git a/SEImplementation/src/lib/Background/SE/ImageMode.cpp b/SEImplementation/src/lib/Background/SE/ImageMode.cpp index 208aa513..fa858d08 100644 --- a/SEImplementation/src/lib/Background/SE/ImageMode.cpp +++ b/SEImplementation/src/lib/Background/SE/ImageMode.cpp @@ -30,11 +30,11 @@ template ImageMode::ImageMode(const std::shared_ptr>& image, const std::shared_ptr>& variance, int cell_w, int cell_h, T invalid_value, T kappa1, T kappa2, T kappa3, - T rtol, size_t max_iter): m_image(image), + T rtol, size_t max_iter, T pearson_factor): m_image(image), 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_rtol(rtol), m_max_iter(max_iter), m_pearson_factor(pearson_factor) { auto hist_width = std::div(image->getWidth(), m_cell_w); if (hist_width.rem) ++hist_width.quot; @@ -106,7 +106,7 @@ std::tuple ImageMode::getBackGuess(const std::vector& data) const { // Sigma is 0 T mode; - mode = 3.5 * median - 2.5 * mean; + mode = m_pearson_factor * median - (m_pearson_factor - 1) * mean; return std::make_tuple(mode, sigma); } diff --git a/SEImplementation/src/lib/Background/SE/SEBackgroundLevelAnalyzer.cpp b/SEImplementation/src/lib/Background/SE/SEBackgroundLevelAnalyzer.cpp index 3c5edac1..3debb953 100644 --- a/SEImplementation/src/lib/Background/SE/SEBackgroundLevelAnalyzer.cpp +++ b/SEImplementation/src/lib/Background/SE/SEBackgroundLevelAnalyzer.cpp @@ -31,8 +31,9 @@ namespace SourceXtractor { SEBackgroundLevelAnalyzer::SEBackgroundLevelAnalyzer(const std::vector& cell_size, const std::vector& smoothing_box, - const WeightImageConfig::WeightType weight_type) - : m_weight_type(weight_type) { + const WeightImageConfig::WeightType weight_type, + SeFloat pearson_factor) + : m_weight_type(weight_type), m_pearson_factor(pearson_factor) { assert(cell_size.size() > 0 && cell_size.size() < 3); assert(smoothing_box.size() > 0 && smoothing_box.size() < 3); m_cell_size[0] = cell_size.front(); @@ -114,7 +115,7 @@ BackgroundModel SEBackgroundLevelAnalyzer::analyzeBackground( } // Create histogram model for the image - ImageMode histo(image, variance_map, m_cell_size[0], m_cell_size[1], mask_value, 2, 5, 3); + ImageMode histo(image, variance_map, m_cell_size[0], m_cell_size[1], mask_value, 2, 5, 3, m_pearson_factor); auto mode = histo.getModeImage(); auto var = histo.getSigmaImage(); diff --git a/SEImplementation/src/lib/Configuration/BackgroundConfig.cpp b/SEImplementation/src/lib/Configuration/BackgroundConfig.cpp index 5574fa24..3fb87473 100644 --- a/SEImplementation/src/lib/Configuration/BackgroundConfig.cpp +++ b/SEImplementation/src/lib/Configuration/BackgroundConfig.cpp @@ -29,13 +29,15 @@ namespace SourceXtractor { static const std::string BACKGROUND_VALUE {"background-value" }; static const std::string THRESHOLD_VALUE {"detection-threshold" }; +static const std::string PEARSON_FACTOR_VALUE {"pearson-factor" }; BackgroundConfig::BackgroundConfig(long manager_id) : Configuration(manager_id), m_background_level_absolute(false), m_background_level(0), m_detection_threshold_absolute(false), - m_detection_threshold(0) { + m_detection_threshold(0), + m_pearson_factor(2.5) { } std::map BackgroundConfig::getProgramOptions() { @@ -44,6 +46,8 @@ std::map BackgroundConfig::ge "Background value to be subtracted from the detection image."}, {THRESHOLD_VALUE.c_str(), po::value()->default_value((1.5)), "Detection threshold above the background."}, + {PEARSON_FACTOR_VALUE.c_str(), po::value()->default_value((2.5)), + "Pearson factor for background mode estimation."}, }}}; } @@ -56,6 +60,9 @@ void BackgroundConfig::initialize(const UserValues& args) { m_detection_threshold_absolute = true; m_detection_threshold = args.find(THRESHOLD_VALUE)->second.as(); } + if (args.find(PEARSON_FACTOR_VALUE) != args.end()) { + m_pearson_factor = args.find(PEARSON_FACTOR_VALUE)->second.as(); + } } } // SourceXtractor namespace From a7f090cfca72c8f629ea45dfbdcdec9a93fea242 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Tue, 9 Jun 2026 14:34:01 +0200 Subject: [PATCH 37/40] fix: pass the right parameter --- .../src/lib/Background/SE/SEBackgroundLevelAnalyzer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SEImplementation/src/lib/Background/SE/SEBackgroundLevelAnalyzer.cpp b/SEImplementation/src/lib/Background/SE/SEBackgroundLevelAnalyzer.cpp index 3debb953..e7fefe2a 100644 --- a/SEImplementation/src/lib/Background/SE/SEBackgroundLevelAnalyzer.cpp +++ b/SEImplementation/src/lib/Background/SE/SEBackgroundLevelAnalyzer.cpp @@ -115,7 +115,7 @@ BackgroundModel SEBackgroundLevelAnalyzer::analyzeBackground( } // Create histogram model for the image - ImageMode histo(image, variance_map, m_cell_size[0], m_cell_size[1], mask_value, 2, 5, 3, m_pearson_factor); + ImageMode histo(image, variance_map, m_cell_size[0], m_cell_size[1], mask_value, 2, 5, 3, 1e-4, 100, m_pearson_factor); auto mode = histo.getModeImage(); auto var = histo.getSigmaImage(); From 8e60b6047c5b43d8e8cbc561b189785a2b98dcb8 Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Tue, 9 Jun 2026 16:54:15 +0200 Subject: [PATCH 38/40] make sure the mode stays 100% identical when using the default value, use double --- .../Background/BackgroundAnalyzerFactory.h | 2 +- .../SEImplementation/Background/SE/ImageMode.h | 5 +++-- .../Configuration/BackgroundConfig.h | 4 ++-- .../src/lib/Background/SE/ImageMode.cpp | 15 ++++++++++++--- 4 files changed, 18 insertions(+), 8 deletions(-) diff --git a/SEImplementation/SEImplementation/Background/BackgroundAnalyzerFactory.h b/SEImplementation/SEImplementation/Background/BackgroundAnalyzerFactory.h index ca20b183..eefa5dff 100644 --- a/SEImplementation/SEImplementation/Background/BackgroundAnalyzerFactory.h +++ b/SEImplementation/SEImplementation/Background/BackgroundAnalyzerFactory.h @@ -50,7 +50,7 @@ class BackgroundAnalyzerFactory : public Euclid::Configuration::Configuration { std::vector m_cell_size; std::vector m_smoothing_box; WeightImageConfig::WeightType m_weight_type; - SeFloat m_pearson_factor; + double m_pearson_factor; }; } diff --git a/SEImplementation/SEImplementation/Background/SE/ImageMode.h b/SEImplementation/SEImplementation/Background/SE/ImageMode.h index f4f2ddaf..9041319d 100644 --- a/SEImplementation/SEImplementation/Background/SE/ImageMode.h +++ b/SEImplementation/SEImplementation/Background/SE/ImageMode.h @@ -68,7 +68,7 @@ class ImageMode { ImageMode(const std::shared_ptr>& image, const std::shared_ptr>& variance, int cell_w, int cell_h, T invalid_value, T kappa1 = 2, T kappa2 = 5, T kappa3 = 3, - T rtol = 1e-4, size_t max_iter = 100, T pearson_factor = 2.5); + T rtol = 1e-4, size_t max_iter = 100, double pearson_factor = 2.5); /** * Destructor @@ -106,7 +106,8 @@ class ImageMode { std::shared_ptr> m_mode, m_sigma; std::shared_ptr> m_var_mode, m_var_sigma; int m_cell_w, m_cell_h; - T m_invalid, m_kappa1, m_kappa2, m_kappa3, m_rtol, m_pearson_factor; + T m_invalid, m_kappa1, m_kappa2, m_kappa3, m_rtol; + double m_pearson_factor; size_t m_max_iter; std::tuple getBackGuess(const std::vector &data) const; diff --git a/SEImplementation/SEImplementation/Configuration/BackgroundConfig.h b/SEImplementation/SEImplementation/Configuration/BackgroundConfig.h index 772b9495..ac659779 100644 --- a/SEImplementation/SEImplementation/Configuration/BackgroundConfig.h +++ b/SEImplementation/SEImplementation/Configuration/BackgroundConfig.h @@ -55,7 +55,7 @@ class BackgroundConfig : public Euclid::Configuration::Configuration { return m_detection_threshold_absolute; } - SeFloat getPearsonFactor() const { + double getPearsonFactor() const { return m_pearson_factor; } @@ -66,7 +66,7 @@ class BackgroundConfig : public Euclid::Configuration::Configuration { bool m_detection_threshold_absolute; SeFloat m_detection_threshold; - SeFloat m_pearson_factor; + double m_pearson_factor; }; } /* namespace SourceXtractor */ diff --git a/SEImplementation/src/lib/Background/SE/ImageMode.cpp b/SEImplementation/src/lib/Background/SE/ImageMode.cpp index fa858d08..65ed307f 100644 --- a/SEImplementation/src/lib/Background/SE/ImageMode.cpp +++ b/SEImplementation/src/lib/Background/SE/ImageMode.cpp @@ -30,7 +30,7 @@ template ImageMode::ImageMode(const std::shared_ptr>& image, const std::shared_ptr>& variance, int cell_w, int cell_h, T invalid_value, T kappa1, T kappa2, T kappa3, - T rtol, size_t max_iter, T pearson_factor): m_image(image), + T rtol, size_t max_iter, double pearson_factor): m_image(image), m_cell_w(cell_w), m_cell_h(cell_h), m_invalid(invalid_value), m_kappa1(kappa1), m_kappa2(kappa2), m_kappa3(kappa3), @@ -106,8 +106,17 @@ std::tuple ImageMode::getBackGuess(const std::vector& data) const { // Sigma is 0 T mode; - mode = m_pearson_factor * median - (m_pearson_factor - 1) * mean; - + if (std::abs(sigma) == 0) { + mode = mean; + } + // Not crowded: mean and median do not differ more than 30% + else if (std::abs((mean - median) / sigma) < 0.3) { + mode = m_pearson_factor * median - (m_pearson_factor - 1.0) * mean; + } + // Crowded case: we use the median + else { + mode = median; + } return std::make_tuple(mode, sigma); } From 3659a3a7fc76a3fa032780ba127e48c1a1ec255a Mon Sep 17 00:00:00 2001 From: Marc Schefer Date: Fri, 19 Jun 2026 16:49:46 +0200 Subject: [PATCH 39/40] Use Alexandria version 2.33.0 Elements 6.3.7 --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index e2ac2962..a94dde9a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -6,4 +6,4 @@ find_package(ElementsProject) #--------------------------------------------------------------- # Declare project name and version -elements_project(SourceXtractorPlusPlus 1.1.0 USE Alexandria 2.32.2 DESCRIPTION "SourceXtractor++, the next generation SExtractor") +elements_project(SourceXtractorPlusPlus 1.1.0 USE Alexandria 2.33.0 DESCRIPTION "SourceXtractor++, the next generation SExtractor") From 17643d941df53674576378135a25c30a45517d3d Mon Sep 17 00:00:00 2001 From: Marc Schefer <31961290+marcschefer@users.noreply.github.com> Date: Tue, 7 Jul 2026 12:25:27 +0200 Subject: [PATCH 40/40] Update README.md --- README.md | 46 ++++------------------------------------------ 1 file changed, 4 insertions(+), 42 deletions(-) diff --git a/README.md b/README.md index 3caad664..91b62cab 100644 --- a/README.md +++ b/README.md @@ -17,13 +17,13 @@ Documentation: https://astrorama.github.io/SourceXtractorPlusPlus/ SourceXtractor++ is available on [Anaconda Cloud for Linux and MacOSX](https://anaconda.org/astrorama/sourcextractor) ```bash -conda install -c conda-forge -c astrorama sourcextractor==1.0.3 +conda install -c conda-forge -c astrorama sourcextractor==1.1.0 ``` We would recommend, however, to install it into its own environment. ```bash -conda create -n sourcex -c astrorama -c conda-forge sourcextractor==1.0.3 +conda create -n sourcex -c astrorama -c conda-forge sourcextractor==1.1.0 conda activate sourcex ``` @@ -36,47 +36,9 @@ conda activate sourcex-dev We recommend always using an explicit version number as there is no guarantee that `conda` will install the latest version otherwise. -## Fedora +Fedora and RPM -SourceXtractor++ is shipped directly in Fedora. It can be installed simply as follows: - -```bash -sudo dnf install sourcextractor++ -``` - -[We also build our own RPMs as part of our CI setup](http://repository.astro.unige.ch/euclid/) -They can be installed adding a `.repo` file under `/etc/yum.repos.d` with the following content: - -```ini -# Stable and tagged versions -# This is the recommended way to install sourcextractor++ -[Astrorama] -name=Astrorama -baseurl=http://repository.astro.unige.ch/euclid/fedora/$releasever/$basearch -enabled=1 -gpgcheck=0 - -# Enable this to get the stable (untagged) version -[Astrorama-Stable] -name=Astrorama Stable -baseurl=http://repository.astro.unige.ch/euclid/devel/stable/fedora/$releasever/$basearch -enabled=0 -gpgcheck=0 - -# Enable this to get the development version -# Use with care! -[Astrorama-Develop] -name=Astrorama Develop -baseurl=http://repository.astro.unige.ch/euclid/devel/develop/fedora/$releasever/$basearch -enabled=0 -gpgcheck=0 -``` - -And install with - -```bash -dnf install SourceXtractorPlusPlus -``` +RPM packages are temporarily not available due to technical problems. ## Test Data Set