diff --git a/.gitignore b/.gitignore index 8b9266d5a..7b4781cce 100644 --- a/.gitignore +++ b/.gitignore @@ -11,3 +11,4 @@ __pycache__ /packages /doc/build/ /build/ +.vscode diff --git a/CMakeLists.txt b/CMakeLists.txt index e658ef8f8..a94dde9a7 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 @@ -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.33.0 DESCRIPTION "SourceXtractor++, the next generation SExtractor") diff --git a/ModelFitting/ModelFitting/Parameters/DependentParameter.h b/ModelFitting/ModelFitting/Parameters/DependentParameter.h index c759a110b..8dc9e8c38 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 */ diff --git a/README.md b/README.md index b04dbdd08..91b62cabe 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.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.2 +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 diff --git a/SEBenchmarks/CMakeLists.txt b/SEBenchmarks/CMakeLists.txt index 8c45f204c..d779a9678 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: @@ -56,6 +54,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 000000000..e3de2eb86 --- /dev/null +++ b/SEBenchmarks/src/program/BenchVariablePsfStack.cpp @@ -0,0 +1,136 @@ +/** 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(3), "Number of timing measurements to take") + ("fits-file", po::value()->default_value(""), "FITS file containing PSF stack"); + 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(); + return Elements::ExitCode::DATAERR; + } + } else { + logger.error() << "No FITS file provided"; + return Elements::ExitCode::USAGE; + } + + 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(); + return Elements::ExitCode::DATAERR; + } + + return Elements::ExitCode::OK; + } +}; + +MAIN_FOR(BenchVariablePsfStack) diff --git a/SEFramework/CMakeLists.txt b/SEFramework/CMakeLists.txt index 775552542..fcf4a3271 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/SEFramework/SEFramework/CoordinateSystem/CoordinateSystem.h b/SEFramework/SEFramework/CoordinateSystem/CoordinateSystem.h index d43e26a42..30605da8b 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/SEFramework/SEFramework/FITS/FitsImageSource.h b/SEFramework/SEFramework/FITS/FitsImageSource.h index 5686df98c..a66107d63 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/SEFramework/Output/OutputRegistry.h b/SEFramework/SEFramework/Output/OutputRegistry.h index 83c5b0ae0..1a6a0d889 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/SEFramework/Psf/VariablePsfStack.h b/SEFramework/SEFramework/Psf/VariablePsfStack.h index 836b9cf11..c71fa0456 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/CoordinateSystem/WCS.cpp b/SEFramework/src/lib/CoordinateSystem/WCS.cpp index ac14c0612..b3c19d364 100644 --- a/SEFramework/src/lib/CoordinateSystem/WCS.cpp +++ b/SEFramework/src/lib/CoordinateSystem/WCS.cpp @@ -255,8 +255,8 @@ 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; } wcsfree(&wcs_copy); return ImageCoordinate(pc_array[0] - 1, pc_array[1] - 1); // -1 as fits standard coordinates start at 1 diff --git a/SEFramework/src/lib/FITS/FitsImageSource.cpp b/SEFramework/src/lib/FITS/FitsImageSource.cpp index f68b06297..776ac420b 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/SEFramework/src/lib/Output/OutputRegistry.cpp b/SEFramework/src/lib/Output/OutputRegistry.cpp index 6cef389b1..cf284cdf4 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/SEFramework/src/lib/Psf/VariablePsf.cpp b/SEFramework/src/lib/Psf/VariablePsf.cpp index c8367c350..6e66e98a9 100644 --- a/SEFramework/src/lib/Psf/VariablePsf.cpp +++ b/SEFramework/src/lib/Psf/VariablePsf.cpp @@ -21,10 +21,14 @@ * Author: Alejandro Álvarez Ayllón */ -#include #include +#include +#include +#include "SEUtils/IsNan.h" + #include "SEFramework/Psf/VariablePsf.h" +static auto stack_logger = Elements::Logging::getLogger("PSFExPsf"); namespace SourceXtractor { @@ -130,6 +134,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 (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/src/lib/Psf/VariablePsfStack.cpp b/SEFramework/src/lib/Psf/VariablePsfStack.cpp index 2276e55d7..9b65da752 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/SEImplementation/CMakeLists.txt b/SEImplementation/CMakeLists.txt index 1b52a93d3..02cc4bafc 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: @@ -33,8 +31,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 +98,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/SEImplementation/SEImplementation/Background/BackgroundAnalyzerFactory.h b/SEImplementation/SEImplementation/Background/BackgroundAnalyzerFactory.h index 96883fcc7..eefa5dffa 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; + double m_pearson_factor; }; } diff --git a/SEImplementation/SEImplementation/Background/SE/ImageMode.h b/SEImplementation/SEImplementation/Background/SE/ImageMode.h index 57f9c045c..9041319d4 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, double pearson_factor = 2.5); /** * Destructor @@ -107,6 +107,7 @@ class ImageMode { 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; + double 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 1fe32982b..1e8117433 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 2f1445e52..ac659779f 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; } + double 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; + + double m_pearson_factor; }; } /* namespace SourceXtractor */ diff --git a/SEImplementation/SEImplementation/Configuration/FileManagerConfig.h b/SEImplementation/SEImplementation/Configuration/FileManagerConfig.h new file mode 100644 index 000000000..08f1e7d80 --- /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/SEImplementation/Measurement/MultithreadedMeasurement.h b/SEImplementation/SEImplementation/Measurement/MultithreadedMeasurement.h index 7e46a1390..3524303d9 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/SEImplementation/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.h b/SEImplementation/SEImplementation/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.h index 9c0ecc3cc..839c2ef51 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/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h b/SEImplementation/SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h index fb7ae567c..586b1f0b5 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,46 +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(PixelCoordinate min_coord, PixelCoordinate max_coord): - m_min_coord{min_coord}, m_max_coord{max_coord}, 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}, m_is_valid{true} { assert(min_coord.m_x <= max_coord.m_x && min_coord.m_y <= max_coord.m_y); } - PixelCoordinate getTopLeft() const { - assert(m_max_coord.m_x >= 0); + ImageCoordinate getTopLeft() const { return m_min_coord; } - PixelCoordinate getBottomRight() const { - assert(m_max_coord.m_x >= 0); + ImageCoordinate getBottomRight() const { return m_max_coord; } - int getWidth() const { - if (m_max_coord.m_x < 0) - return 0; - return m_max_coord.m_x - m_min_coord.m_x + 1; + double getWidth() const { + 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; + double getHeight() const { + return m_max_coord.m_y - m_min_coord.m_y; } - PixelRectangle getRect() const { - return PixelRectangle(m_min_coord, m_max_coord); + bool isBadProjection() const { + return m_bad_projection; } - bool badProjection() const { - return m_bad_projection; + bool isValid() const { + return m_is_valid; + } + + bool isEmpty() const { + return getWidth() <= 0.0 || getHeight() <= 0.0; } private: - PixelCoordinate m_min_coord, m_max_coord; + ImageCoordinate m_min_coord, m_max_coord; bool m_bad_projection; + bool m_is_valid; }; } // end SourceXtractor diff --git a/SEImplementation/src/lib/Background/BackgroundAnalyzerFactory.cpp b/SEImplementation/src/lib/Background/BackgroundAnalyzerFactory.cpp index 3c5d31ff4..f7a19e19a 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 857c1bb42..65ed307f8 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, 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), - 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; @@ -111,7 +111,7 @@ std::tuple ImageMode::getBackGuess(const std::vector& data) const { } // 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; + mode = m_pearson_factor * median - (m_pearson_factor - 1.0) * mean; } // Crowded case: we use the median else { diff --git a/SEImplementation/src/lib/Background/SE/SEBackgroundLevelAnalyzer.cpp b/SEImplementation/src/lib/Background/SE/SEBackgroundLevelAnalyzer.cpp index 3c5edac12..e7fefe2a2 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, 1e-4, 100, 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 5574fa240..3fb874739 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 diff --git a/SEImplementation/src/lib/Configuration/FileManagerConfig.cpp b/SEImplementation/src/lib/Configuration/FileManagerConfig.cpp new file mode 100644 index 000000000..e37c683c9 --- /dev/null +++ b/SEImplementation/src/lib/Configuration/FileManagerConfig.cpp @@ -0,0 +1,58 @@ +/** 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() { + 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/SEImplementation/src/lib/Measurement/MultithreadedMeasurement.cpp b/SEImplementation/src/lib/Measurement/MultithreadedMeasurement.cpp index b29fdecff..b0d7046aa 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; } @@ -118,23 +121,20 @@ 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)); - } + { + 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)); } } diff --git a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp index 1e99d9eae..19b64a5eb 100644 --- a/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp +++ b/SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp @@ -74,7 +74,13 @@ FlexibleModelFittingIterativeTask::~FlexibleModelFittingIterativeTask() { } PixelRectangle FlexibleModelFittingIterativeTask::getUnclippedFittingRect(SourceInterface& source, int frame_index) const { - auto fitting_rect = source.getProperty(frame_index).getRect(); + auto& measurement_frame_rectangle = source.getProperty(frame_index); + 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); @@ -84,14 +90,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.0) || (max_coord.m_y - min_coord.m_y <= 0.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,12 +131,19 @@ PixelRectangle FlexibleModelFittingIterativeTask::getUnclippedFittingRect(Source max.m_y = min.m_y + size; } - return PixelRectangle(min, max); + auto min_pc = PixelCoordinate(static_cast(min.m_x), static_cast(min.m_y)); + auto max_pc = PixelCoordinate(static_cast(max.m_x), static_cast(max.m_y)); + + return PixelRectangle(min_pc, max_pc); } } 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(); @@ -413,6 +426,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)); } @@ -488,45 +503,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) { + if (index != source_index && isFrameValid(src, frame->getFrameNb())) { + // 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++; } @@ -537,7 +534,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); @@ -548,8 +545,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)); @@ -704,17 +701,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 @@ -724,7 +723,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; } @@ -741,14 +740,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; @@ -760,49 +759,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); @@ -842,6 +826,7 @@ void FlexibleModelFittingIterativeTask::updateCheckImages(SourceGroupInterface& } } + index++; } } diff --git a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp index f57d56fb4..7594eee65 100644 --- a/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp +++ b/SEImplementation/src/lib/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangleTask.cpp @@ -44,9 +44,14 @@ 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; try { coord1 = measurement_frame_coordinates->worldToImage( detection_frame_coordinates->imageToWorld(ImageCoordinate(stamp_top_left.m_x, stamp_top_left.m_y))); @@ -58,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_coordinates = true; + source.setIndexedProperty(m_instance, true); + return; } // Determine the min/max coordinates @@ -67,26 +73,32 @@ 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()) { - source.setIndexedProperty(m_instance, bad_coordinates); + // 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; } + + // 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 - 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() - 1), max_x); + max_y = std::min(static_cast(measurement_frame_info.getHeight() - 1), max_y); - source.setIndexedProperty(m_instance, min_coord, max_coord); + 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 68caedfc5..3d7778bf1 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,26 +32,30 @@ 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 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 w = assoc_mode.getRefFramePixelWidth(); + auto 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))); } catch (const InvalidCoordinatesException&) { - bad_coordinates = true; + source.setIndexedProperty(m_instance, true); + return; } // Determine the min/max coordinates @@ -59,26 +64,32 @@ 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, bad_coordinates); + // 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; } + + // 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 - 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() - 1), max_x); + max_y = std::min(static_cast(measurement_frame_info.getHeight() - 1), max_y); - source.setIndexedProperty(m_instance, min_coord, max_coord); + 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/Psf/PsfPluginConfig.cpp b/SEImplementation/src/lib/Plugin/Psf/PsfPluginConfig.cpp index 462835f93..e0c61a9da 100644 --- a/SEImplementation/src/lib/Plugin/Psf/PsfPluginConfig.cpp +++ b/SEImplementation/src/lib/Plugin/Psf/PsfPluginConfig.cpp @@ -124,10 +124,13 @@ static std::shared_ptr readPsfEx(std::unique_ptr &pFi } return std::make_shared(pixel_sampling, components, group_degrees, coefficients); - } catch (CCfits::FITS::NoSuchHDU&) { + } 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 (...) { + logger.error() << "Failed loading a psf file: " << pFits->name(); + throw; } } diff --git a/SEImplementation/src/lib/Plugin/SourceFlags/SourceFlagsSourceTask.cpp b/SEImplementation/src/lib/Plugin/SourceFlags/SourceFlagsSourceTask.cpp index 00598a196..a86461a0b 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; } } diff --git a/SEMain/CMakeLists.txt b/SEMain/CMakeLists.txt index 8fdb92a86..b904bc4f4 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/SEMain/src/program/SourceXtractor.cpp b/SEMain/src/program/SourceXtractor.cpp index 49b2d4937..e79ae077e 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(), diff --git a/SEUtils/CMakeLists.txt b/SEUtils/CMakeLists.txt index 9f3c453aa..8203f66f7 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: @@ -66,6 +64,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/IsNan.h b/SEUtils/SEUtils/IsNan.h new file mode 100644 index 000000000..c4dea6479 --- /dev/null +++ b/SEUtils/SEUtils/IsNan.h @@ -0,0 +1,80 @@ +#ifndef _SEUTILS_ISNAN_H +#define _SEUTILS_ISNAN_H + +#include + +namespace SourceXtractor { + +/* + * From: + * https://github.com/searchivarius/BlogCode/tree/2d947d8c42ab0e10f1f28a2ad036eee10389997a/2017/6/1 + * + * IEEE 754 compliant simple functions to test for NANs and INFs. + * This functions are necessary, b/c isnan doesn't work with -Ofast -ffast-math flags + * see http://searchivarius.org/blog/gcc-disables-isnan-and-isinf-when-compiling-ffast-math-flag + * + * See also test files. The file ./regular is compiled without -ffast-math flag and it + * checks for a large number of values that the output of my functions is the same + * as the output of standard functions. For single precision numbers, i.e., for floats + * these checks are exhaustive. That is, I go over the set of all 4B+ possible float values. + * For doubles, this is not possible, so I use tests where the lower 32 bits of mantissa + * are set to zero. + * + * Copyright (c) 2017 Leonid Boytsov + * + * This code is released under the + * Apache License Version 2.0 http://www.apache.org/licenses/. + * + */ + +// A mask to extract an exponent from the single-precision floating point number +// 01111111100000000000000000000000 +const unsigned FLOAT_EXP_MASK = 0x7F800000; +// A mask to extract a mantissa/fractional part from the single-precision floating point number +// 00000000011111111111111111111111 +const unsigned FLOAT_FRAC_PART_MASK = 0x7FFFFF; + +inline bool fastmath_isnan(float x) { + union { + uint32_t u; + float f; + } conv; + conv.f = x; + return ((conv.u & FLOAT_EXP_MASK) == FLOAT_EXP_MASK) && ((conv.u & FLOAT_FRAC_PART_MASK) != 0); +}; + +inline bool fastmath_isinf(float x) { + union { + uint32_t u; + float f; + } conv; + conv.f = x; + return ((conv.u & FLOAT_EXP_MASK) == FLOAT_EXP_MASK) && ((conv.u & FLOAT_FRAC_PART_MASK) == 0); +}; + +// 0111 1111 1111 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 +const uint64_t DOUBLE_EXP_MASK = 0x7FF0000000000000ul; +// 0000 0000 0000 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 +const uint64_t DOUBLE_FRAC_PART_MASK = 0x000FFFFFFFFFFFFFul; + +inline bool fastmath_isnan(double x) { + union { + uint64_t u; + double f; + } conv; + conv.f = x; + return ((conv.u & DOUBLE_EXP_MASK) == DOUBLE_EXP_MASK) && ((conv.u & DOUBLE_FRAC_PART_MASK) != 0); +}; + +inline bool fastmath_isinf(double x) { + union { + uint64_t u; + double f; + } conv; + conv.f = x; + return ((conv.u & DOUBLE_EXP_MASK) == DOUBLE_EXP_MASK) && ((conv.u & DOUBLE_FRAC_PART_MASK) == 0); +}; + +} // namespace SourceXtractor + +#endif diff --git a/SEUtils/SEUtils/KdTree.h b/SEUtils/SEUtils/KdTree.h index 5c4fc3546..3e739c2eb 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 606b939cf..ce45e4d22 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 000000000..5c5cf9e1b --- /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 () + +} diff --git a/cmake/modules/FindBoostPython.cmake b/cmake/modules/FindBoostPython.cmake deleted file mode 100644 index 69fbcc596..000000000 --- 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 now 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 ()