From 73541d6582aed7f66e2d95f7cc03ac1a79ef0dec Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Thu, 14 Aug 2025 11:21:08 +0200 Subject: [PATCH 01/29] Basic libasdf detection at cmake-time and minimal stub for AsdfFile to test Nothing that does anything here yet, just bare minimum to test libasdf detection and compile something with it. --- SEFramework/CMakeLists.txt | 24 ++++++- SEFramework/SEFramework/ASDF/AsdfFile.h | 59 +++++++++++++++ SEFramework/auxdir/basic.asdf | Bin 0 -> 824 bytes SEFramework/src/lib/ASDF/AsdfFile.cpp | 72 +++++++++++++++++++ SEFramework/tests/src/ASDF/AsdfFile_test.cpp | 56 +++++++++++++++ cmake/modules/FindASDF.cmake | 41 +++++++++++ 6 files changed, 250 insertions(+), 2 deletions(-) create mode 100644 SEFramework/SEFramework/ASDF/AsdfFile.h create mode 100644 SEFramework/auxdir/basic.asdf create mode 100644 SEFramework/src/lib/ASDF/AsdfFile.cpp create mode 100644 SEFramework/tests/src/ASDF/AsdfFile_test.cpp create mode 100644 cmake/modules/FindASDF.cmake diff --git a/SEFramework/CMakeLists.txt b/SEFramework/CMakeLists.txt index 775552542..a4ef0f751 100644 --- a/SEFramework/CMakeLists.txt +++ b/SEFramework/CMakeLists.txt @@ -31,10 +31,22 @@ elements_depends_on_subdirs(ModelFitting) find_package(GMock) find_package(CCfits) find_package(BoostDLL) +find_package(ASDF) find_package(FFTW COMPONENTS single double) find_package(WCSLIB REQUIRED) +#=============================================================================== +# ASDF support is optional--don't compile it if libasdf was not found +#=============================================================================== +if(NOT ASDF_FOUND) + message("libasdf not found, compiling without ASDF support!") +else() + file(GLOB ASDF_SRC src/lib/ASDF/*.cpp) +endif() + + + #=============================================================================== # Declare the library dependencies here # Example: @@ -57,10 +69,12 @@ elements_add_library(SEFramework src/lib/FITS/*.cpp src/lib/Frame/*.cpp src/lib/CoordinateSystem/*.cpp + ${ASDF_SRC} LINK_LIBRARIES ElementsKernel SEUtils Table Configuration MathUtils FilePool - WCSLIB ${CCFITS_LIBRARIES} ${FFTW_LIBRARIES} - INCLUDE_DIRS SEUtils ${CCFITS_INCLUDE_DIRS} ${BoostDLL_INCLUDE_DIRS} ${FFTW_INCLUDE_DIRS} + WCSLIB ${CCFITS_LIBRARIES} ${FFTW_LIBRARIES} ${ASDF_LIBRARIES} + INCLUDE_DIRS SEUtils ${CCFITS_INCLUDE_DIRS} ${BoostDLL_INCLUDE_DIRS} + ${FFTW_INCLUDE_DIRS} ${ASDF_INCLUDE_DIRS} PUBLIC_HEADERS SEFramework) #=============================================================================== @@ -187,6 +201,12 @@ elements_add_unit_test(TemporaryFitsSource_test tests/src/FITS/TemporaryFitsSour elements_add_unit_test(WCS_test tests/src/CoordinateSystem/WCS_test.cpp LINK_LIBRARIES SEFramework TYPE Boost) + +if(ASDF_FOUND) + elements_add_unit_test(AsdfFile_test tests/src/ASDF/AsdfFile_test.cpp + LINK_LIBRARIES SEFramework + TYPE Boost) +endif() #=============================================================================== # Declare the Python programs here # Examples : diff --git a/SEFramework/SEFramework/ASDF/AsdfFile.h b/SEFramework/SEFramework/ASDF/AsdfFile.h new file mode 100644 index 000000000..3546b97f9 --- /dev/null +++ b/SEFramework/SEFramework/ASDF/AsdfFile.h @@ -0,0 +1,59 @@ +/** Copyright © 2019 Université de Genève, LMU Munich - Faculty of Physics, IAP-CNRS/Sorbonne Université + * + * This library is free software; you can redistribute it and/or modify it under + * the terms of the GNU Lesser General Public License as published by the Free + * Software Foundation; either version 3.0 of the License, or (at your option) + * any later version. + * + * This library is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + */ + +/* + * AsdfFile.h + * + * Created on: Aug 14, 2025 + * Author: embray + */ + +#ifndef _SEFRAMEWORK_ASDF_ASDFFILE_H_ +#define _SEFRAMEWORK_ASDF_ASDFFILE_H_ + +#include +#include + +#include "SEFramework/Image/ImageSourceWithMetadata.h" + +namespace SourceXtractor { + +/** + * @class AsdfFile + * @brief represents access to a whole ASDF file + * + */ +class AsdfFile { +public: + AsdfFile(const boost::filesystem::path& path); + + AsdfFile(AsdfFile&&) = default; + + virtual ~AsdfFile(); + + asdf_file_t* getAsdfFilePtr(); + +private: + boost::filesystem::path m_path; + std::unique_ptr m_asdf_ptr; + + void open(); +}; + +} // namespace SourceXtractor + +#endif /* _SEFRAMEWORK_ASDF_ASDFFILE_H_ */ diff --git a/SEFramework/auxdir/basic.asdf b/SEFramework/auxdir/basic.asdf new file mode 100644 index 0000000000000000000000000000000000000000..b39b0f5a7a44fb60a25baf1bb9552f2bd27be268 GIT binary patch literal 824 zcma)5&2H2%5N@GRDX+j~s;cy|adxG$LOGF=QdLchME5|2s>qvoHDNl!6R_vzzc9;oH)yZgw$N-$KT9+{>IsjF5-7EVna6ccCDXYOrrBR`VdocOv$~= z=qw?BjeSB>rFGL%PBqud&>b5!!yB#z zL1-e}@&*f~Lzol-S|+i$7NuGUZ3a-NN-X$Xgs@*AF2W#~ODvix%T(13nkijzZ2kR# z2T=2h(01?uuBnGqbqEKR&7D_BgNm3COngRB+e2fcp86`ARswm>k#{T!vR0C2B{!xU z!d4|Lt#~bG!XP_y&QU+)pZZ42uA1)Cs?}YBDLoHb30F-`8xsk;2TP)#k@hnFPu4#= zqy~ANbDggVV%IaeYdb!^(OSI|)GE{yup=~Iy3GR^)#!{?VR>5|E*BKZ8pj8enJ!T% zoeNrgDG^ISIa6>8`SkkAV~ly<4wKW3;l|(J8yBwZtIOme{qo}bkLdp8*W}mF9=_?E rC(b@~cFWmk&OUc`+u8OHfWzed+fz6>kB>gCJ%4)Ro)2TXfaUia>;UQY literal 0 HcmV?d00001 diff --git a/SEFramework/src/lib/ASDF/AsdfFile.cpp b/SEFramework/src/lib/ASDF/AsdfFile.cpp new file mode 100644 index 000000000..5aeccdc26 --- /dev/null +++ b/SEFramework/src/lib/ASDF/AsdfFile.cpp @@ -0,0 +1,72 @@ +/** Copyright © 2019 Université de Genève, LMU Munich - Faculty of Physics, IAP-CNRS/Sorbonne Université + * + * This library is free software; you can redistribute it and/or modify it under + * the terms of the GNU Lesser General Public License as published by the Free + * Software Foundation; either version 3.0 of the License, or (at your option) + * any later version. + * + * This library is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + */ + +/* + * AsdfFile.cpp + * + * Created on: Aug 14, 2025 + * Author: embray + */ + +#include + +#include + +#include + +#include "ElementsKernel/Exception.h" +#include "ElementsKernel/Logging.h" + +#include "SEFramework/ASDF/AsdfFile.h" + +namespace SourceXtractor { + +static Elements::Logging logger = Elements::Logging::getLogger("AsdfFile"); + + +AsdfFile::AsdfFile(const boost::filesystem::path& path) + : m_path(path), m_asdf_ptr(nullptr, asdf_close) { + + open(); +} + +AsdfFile::~AsdfFile() {} + +asdf_file_t* AsdfFile::getAsdfFilePtr() { + return m_asdf_ptr.get(); +} + +void AsdfFile::open() { + asdf_file_t* ptr = asdf_open_file(m_path.native().c_str(), "r"); + + // TODO: Actually kinda stinks there's not enough options for retrieving file open errors + // yet in libasdf; this needs improvement on the libasdf side. + if (ptr == nullptr) { + throw Elements::Exception() + << "Can't open ASDF file: "; + } else { + // Check if the file was created but has an error condition set + const char *error_message = asdf_error(ptr); + + if (error_message != nullptr) { + throw Elements::Exception() + << "Can't open ASDF file: " << m_path << " reason: " << error_message; + } + } +} + +} // namespace SourceXtractor diff --git a/SEFramework/tests/src/ASDF/AsdfFile_test.cpp b/SEFramework/tests/src/ASDF/AsdfFile_test.cpp new file mode 100644 index 000000000..faacaea08 --- /dev/null +++ b/SEFramework/tests/src/ASDF/AsdfFile_test.cpp @@ -0,0 +1,56 @@ +/** Copyright © 2019 Université de Genève, LMU Munich - Faculty of Physics, IAP-CNRS/Sorbonne Université + * + * This library is free software; you can redistribute it and/or modify it under + * the terms of the GNU Lesser General Public License as published by the Free + * Software Foundation; either version 3.0 of the License, or (at your option) + * any later version. + * + * This library is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + */ +/** + * @file tests/src/FitsReader_test.cpp + * @date 06/14/16 + * @author nikoapos + */ + +#include + +#include +#include + +#include "SEFramework/ASDF/AsdfFile.h" + +using namespace SourceXtractor; + +//----------------------------------------------------------------------------- + +BOOST_AUTO_TEST_SUITE (AsdfFile_test) + + +//----------------------------------------------------------------------------- + +BOOST_AUTO_TEST_CASE( open_file ) { + BOOST_CHECK_NO_THROW({ + AsdfFile asdf_file(Elements::getAuxiliaryPath("basic.asdf").native()); + }); +} + +//----------------------------------------------------------------------------- + +BOOST_AUTO_TEST_CASE( missing_file ) { + BOOST_CHECK_THROW({ + AsdfFile asdf_file("/does/not/exist"); + }, Elements::Exception); +} + +//----------------------------------------------------------------------------- + + +BOOST_AUTO_TEST_SUITE_END () diff --git a/cmake/modules/FindASDF.cmake b/cmake/modules/FindASDF.cmake new file mode 100644 index 000000000..28351906a --- /dev/null +++ b/cmake/modules/FindASDF.cmake @@ -0,0 +1,41 @@ +# - Locate the libasdf library +# Defines: +# +# ASDF_FOUND +# ASDF_INCLUDE_DIR +# ASDF_INCLUDE_DIRS (not cached) +# ASDF_LIBRARY +# ASDF_LIBRARIES (not cached) + +if(NOT ASDF_FOUND) + + # Look for the header + find_path(ASDF_INCLUDE_DIR asdf.h + HINTS ENV ASDF_ROOT_DIR ASDF_INSTALL_DIR + PATH_SUFFIXES include) + + # Look for the library + find_library(ASDF_LIBRARY asdf + HINTS ENV ASDF_ROOT_DIR ASDF_INSTALL_DIR + PATH_SUFFIXES lib) + + # Look for libfyaml which is a dependency of libasdf + find_library(FYAML_LIBRARY fyaml + HINTS ENV ASDF_ROOT_DIR ASDF_INSTALL_DIR + PATH_SUFFIXES lib) + + # You can add other deps here if libasdf needs them at link time + set(ASDF_LIBRARIES ${ASDF_LIBRARY} ${FYAML_LIBRARY}) + set(ASDF_INCLUDE_DIRS ${ASDF_INCLUDE_DIR}) + + include(FindPackageHandleStandardArgs) + find_package_handle_standard_args(ASDF + DEFAULT_MSG + ASDF_INCLUDE_DIRS ASDF_LIBRARIES) + + mark_as_advanced(ASDF_FOUND ASDF_INCLUDE_DIRS ASDF_LIBRARIES) + + list(REMOVE_DUPLICATES ASDF_LIBRARIES) + list(REMOVE_DUPLICATES ASDF_INCLUDE_DIRS) + +endif() From 72f0d1b60922cac884bbb9aef67c3e80f5edd990 Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Thu, 14 Aug 2025 11:23:49 +0200 Subject: [PATCH 02/29] fix file attribution --- SEFramework/tests/src/ASDF/AsdfFile_test.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/SEFramework/tests/src/ASDF/AsdfFile_test.cpp b/SEFramework/tests/src/ASDF/AsdfFile_test.cpp index faacaea08..acf4456e2 100644 --- a/SEFramework/tests/src/ASDF/AsdfFile_test.cpp +++ b/SEFramework/tests/src/ASDF/AsdfFile_test.cpp @@ -15,9 +15,9 @@ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ /** - * @file tests/src/FitsReader_test.cpp - * @date 06/14/16 - * @author nikoapos + * @file tests/src/AsdfFile_test.cpp + * @date 08/14/25 + * @author embray */ #include From 6fb744de6cd182d4fc166c70342633735e27fcf4 Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Mon, 1 Sep 2025 14:33:48 +0200 Subject: [PATCH 03/29] First prototype of ImageFileReader interface Extend the existing (barely used) FitsReader class to implement ImageFileReader Extend FitsImageSource to also allow selecting FITS HDUs by extname (should probably support extver as well, but this is less common) --- .../SEFramework/FITS/FitsImageSource.h | 13 +- SEFramework/SEFramework/FITS/FitsReader.h | 22 +- .../SEFramework/Image/ImageFileReader.h | 194 +++++++++++++++++ SEFramework/src/lib/FITS/FitsImageSource.cpp | 26 ++- SEFramework/src/lib/FITS/FitsReader.cpp | 72 +++++++ SEFramework/src/lib/Image/ImageFileReader.cpp | 198 ++++++++++++++++++ .../tests/src/FITS/FitsReader_test.cpp | 25 ++- .../lib/Configuration/SegmentationConfig.cpp | 4 +- 8 files changed, 537 insertions(+), 17 deletions(-) create mode 100644 SEFramework/SEFramework/Image/ImageFileReader.h create mode 100644 SEFramework/src/lib/FITS/FitsReader.cpp create mode 100644 SEFramework/src/lib/Image/ImageFileReader.cpp diff --git a/SEFramework/SEFramework/FITS/FitsImageSource.h b/SEFramework/SEFramework/FITS/FitsImageSource.h index a66107d63..bedf630f1 100644 --- a/SEFramework/SEFramework/FITS/FitsImageSource.h +++ b/SEFramework/SEFramework/FITS/FitsImageSource.h @@ -58,7 +58,13 @@ 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 = getFileManager()); + std::shared_ptr manager = getFileManager()) : + FitsImageSource(filename, hdu_number, std::nullopt, image_type, std::move(manager)) {} + + explicit FitsImageSource(const std::string& filename, const std::string& extname, + ImageTile::ImageType image_type = ImageTile::AutoType, + std::shared_ptr manager = getFileManager()) : + FitsImageSource(filename, 0, extname, image_type, std::move(manager)) {} FitsImageSource(const std::string& filename, int width, int height, ImageTile::ImageType image_type, @@ -120,7 +126,12 @@ class FitsImageSource : public ImageSource, public std::enable_shared_from_this< void setMetadata(const std::string& key, const MetadataEntry& value) override; private: + FitsImageSource(const std::string& filename, int hdu_number, + std::optional extname, + ImageTile::ImageType image_type, + std::shared_ptr manager); void switchHdu(fitsfile *fptr, int hdu_number) const; + void switchHdu(fitsfile *fptr, const std::string& extname) const; int getDataType() const; diff --git a/SEFramework/SEFramework/FITS/FitsReader.h b/SEFramework/SEFramework/FITS/FitsReader.h index 5ca785741..f6195c5be 100644 --- a/SEFramework/SEFramework/FITS/FitsReader.h +++ b/SEFramework/SEFramework/FITS/FitsReader.h @@ -16,14 +16,15 @@ */ /** * @file SEFramework/Image/FitsReader.h - * @date 06/14/16 - * @author nikoapos + * @date 09/01/25 + * @author embray */ #ifndef _SEFRAMEWORK_IMAGE_FITSREADER_H #define _SEFRAMEWORK_IMAGE_FITSREADER_H #include "SEFramework/Image/BufferedImage.h" +#include "SEFramework/Image/ImageFileReader.h" #include "SEFramework/FITS/FitsImageSource.h" namespace SourceXtractor { @@ -33,16 +34,17 @@ namespace SourceXtractor { * @brief * */ -template -class FitsReader { +class FitsReader : public ImageFileReader { public: - - /** - * @brief Destructor - */ - virtual ~FitsReader() = default; - + using ImageFileReader::ImageFileReader; + + static bool test(std::istream& stream); + std::shared_ptr get() override; + std::shared_ptr get(int image_index) override; + std::shared_ptr get(const std::string& image_path) override; + // TODO: Integrate this interface into the ImageFileReader base class + template static std::shared_ptr> readFile(const std::string& filename) { auto image_source = std::make_shared(filename, 0, ImageTile::getTypeValue(T())); return BufferedImage::create(image_source); diff --git a/SEFramework/SEFramework/Image/ImageFileReader.h b/SEFramework/SEFramework/Image/ImageFileReader.h new file mode 100644 index 000000000..f1fa2713c --- /dev/null +++ b/SEFramework/SEFramework/Image/ImageFileReader.h @@ -0,0 +1,194 @@ +/** Copyright © 2019 Université de Genève, LMU Munich - Faculty of Physics, IAP-CNRS/Sorbonne Université + * + * This library is free software; you can redistribute it and/or modify it under + * the terms of the GNU Lesser General Public License as published by the Free + * Software Foundation; either version 3.0 of the License, or (at your option) + * any later version. + * + * This library is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + */ +/** + * @file ImageFileReader.h + * @date Aug 14, 2025 + * @author embray + */ + +#ifndef _SEFRAMEWORK_IMAGE_IMAGEFILEREADER_H_ +#define _SEFRAMEWORK_IMAGE_IMAGEFILEREADER_H_ + +#include + +#include "SEFramework/Image/ImageSource.h" + +namespace SourceXtractor { + +/** + * @class ImageFileReader + * + * @brief + * Registry of supported file types from which ImageSources can be loaded + * + * Subclasses should implement the methods for returning / iterating individual ImageSources + * from a file. + * + * @details + * File type-specific ImageFileReaders (e.g. FitsImageFileReader) can register themselves at + * compilation time using StaticFileType: This accepts a name for the file format, a vector + * of supported filename extensions (this is used only as a hint for prioritizing which file + * type to assume), and a test function that should be able to rapidly check the file type based on + * its contents, and finally a factory function for constructing the ImageFileReader itself. + * + * The test method simply receives an open read-only file stream and should return a boolean + * if the test passes. + * + * The factory method receives the filename which may include an optional [...] suffix containing + * either an integer or a string. For example, image.fits[1] receives an image_index 1, which + * in the case of a FITS file indicates loading the image from the first HDU. + * + * Currently only supports opening files for reading, though writing will be a future extension. + */ +class ImageFileReader { +public: + ImageFileReader(const std::string& filename); + ImageFileReader(const std::string& filename, const std::string& image_path); + ImageFileReader(const std::string& filename, int image_index); + + virtual ~ImageFileReader() = default; + + virtual std::shared_ptr get(); + virtual std::shared_ptr get(int image_index) = 0; + virtual std::shared_ptr get(const std::string& image_path) = 0; + + /* ImageFileReader iterator interface */ + class Iterator { + public: + using iterator_category = std::input_iterator_tag; + using value_type = std::shared_ptr; + using difference_type = std::ptrdiff_t; + using pointer = std::shared_ptr*; + using reference = std::shared_ptr&; + + Iterator(ImageFileReader* reader, int index) + : m_reader(reader), m_index(index) { + advance(); + } + + std::shared_ptr operator*() const { return m_current; } + + Iterator& operator++() { + ++m_index; + advance(); + return *this; + } + + bool operator==(const Iterator& other) const { + return m_current == other.m_current; + } + + bool operator!=(const Iterator& other) const { + return !(*this == other); + } + + private: + void advance() { + try { + m_current = m_reader->get(m_index); + } catch (...) { + m_current = nullptr; + } + } + + ImageFileReader* m_reader; + int m_index; + std::shared_ptr m_current; + }; + + Iterator begin() { + return Iterator(this, 0); + } + + Iterator end() { + return Iterator(this, -1); + } + + /** + * Test methods receive an open handle to the file (if the file existed) and should read + * the file contents (e.g. for magic bytes) to check if the file type matches expecations. + */ + using TestMethod = std::function; + + /** + * Factory methods must return a shared pointer to an ImageFileReader, and receive a string + * containing the filename + */ + using FactoryMethod = std::function(const std::string&)>; + + /** + * @return A list of known file types. + */ + static std::vector getFileTypes(); + + /** + * @return The default file type to use (currently "fits", when in doubt) + */ + static std::string getDefault(); + + /** + * Create an instance of an `ImageFileReader` from the filename (which may contain an optional + * extension suffix) + * @param filename + * The filename of the file + * @return + * A new instance of the ImageSource. An Elements::Exception is thrown if the file type + * cannot be determined. + */ + static std::unique_ptr create(const std::string &filename); + + static bool test(std::istream& stream); + + struct ImageFileType { + std::string name; + std::vector extensions; + ImageFileReader::TestMethod test; + ImageFileReader::FactoryMethod factory; + }; + + /** + * Register a new file type reader + * @param name + * The name of the file type (e.g. "fits"). Case insensitive. + * @param extensions + * Vector of filename extensions (without the .), e.g. "fits", "fit". If the filename + * matches one of these extensions, then this file type will be checked first (this is + * only a hint though). + */ + template + struct StaticFileType : ImageFileType { + StaticFileType(const std::string& name, + const std::vector& extensions) + : ImageFileType{ + name, extensions, &Reader::test, + [](const std::string& filename) { return std::make_unique(filename); } + } { + registerFileType(*this); + } + }; + + static void registerFileType(ImageFileType file_type); + +protected: + std::string m_filename; + std::optional m_image_path; + int m_image_index; +}; + +} // end of namespace SourceXtractor + +#endif /* _SEFRAMEWORK_IMAGE_IMAGEFILEREADER_H_ */ diff --git a/SEFramework/src/lib/FITS/FitsImageSource.cpp b/SEFramework/src/lib/FITS/FitsImageSource.cpp index 776ac420b..72bd4a597 100644 --- a/SEFramework/src/lib/FITS/FitsImageSource.cpp +++ b/SEFramework/src/lib/FITS/FitsImageSource.cpp @@ -78,6 +78,7 @@ ImageTile::ImageType convertImageType(int bitpix) { } FitsImageSource::FitsImageSource(const std::string& filename, int hdu_number, + std::optional extname, ImageTile::ImageType image_type, std::shared_ptr manager) : m_filename(filename) @@ -91,12 +92,16 @@ FitsImageSource::FitsImageSource(const std::string& filename, int hdu_number, auto acc = m_handler->getAccessor(); auto fptr = acc->m_fd.getFitsFilePtr(); - if (m_hdu_number <= 0) { + if (m_hdu_number <= 0 && !extname) { if (fits_get_hdu_num(fptr, &m_hdu_number) < 0) { throw Elements::Exception() << "Can't get the active HDU from the FITS file: " << filename; } - } - else { + } else if (extname) { + switchHdu(fptr, *extname); + if (fits_get_hdu_num(fptr, &m_hdu_number) < 0) { + throw Elements::Exception() << "Can't get the active HDU from the FITS file: " << filename; + } + } else { switchHdu(fptr, m_hdu_number); } @@ -278,6 +283,21 @@ void FitsImageSource::switchHdu(fitsfile *fptr, int hdu_number) const { } } + +void FitsImageSource::switchHdu(fitsfile *fptr, const std::string& extname) const { + int status = 0; + + fits_movnam_hdu(fptr, IMAGE_HDU, const_cast(extname.c_str()), 0, &status); + + if (status != 0) { + char error_message[32]; + fits_get_errstatus(status, error_message); + throw Elements::Exception() << "Could not switch to HDU with EXTNAME " << extname << " in file " << m_filename + << " status: " << status << " = " << error_message; + } +} + + void FitsImageSource::setLayer(int layer) { if (layer < 0 && layer >= m_depth) { throw Elements::Exception() << "Trying to access an inexistent data cube layer (" << layer << ") in " << m_filename; diff --git a/SEFramework/src/lib/FITS/FitsReader.cpp b/SEFramework/src/lib/FITS/FitsReader.cpp new file mode 100644 index 000000000..5d6a1230d --- /dev/null +++ b/SEFramework/src/lib/FITS/FitsReader.cpp @@ -0,0 +1,72 @@ +/** Copyright © 2019 Université de Genève, LMU Munich - Faculty of Physics, IAP-CNRS/Sorbonne Université + * + * This library is free software; you can redistribute it and/or modify it under + * the terms of the GNU Lesser General Public License as published by the Free + * Software Foundation; either version 3.0 of the License, or (at your option) + * any later version. + * + * This library is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + */ +/** + * @file SEFramework/Image/FitsReader.cpp + * @date 09/01/25 + * @author embray + */ + +#include + +#include "SEFramework/FITS/FitsReader.h" +#include "SEFramework/Image/ImageFileReader.h" +#include "SEFramework/Image/ImageSource.h" + + +namespace SourceXtractor { + + +bool FitsReader::test(std::istream& stream) { + char header[9] = {0}; + if (!stream.read(header, 9)) { + return false; + } + + return std::strncmp(header, "SIMPLE =", 9) == 0; +} + + +std::shared_ptr FitsReader::get() { + if (0 == m_image_index) { + // If the primary HDU was not an image try the first extension HDU + try { + return get(1); + } catch (...) { + return get(2); + } + } + + return ImageFileReader::get(); +} + + +std::shared_ptr FitsReader::get(int hdu_num) { + return std::make_shared(m_filename, hdu_num); +} + + +std::shared_ptr FitsReader::get(const std::string& extname) { + return std::make_shared(m_filename, extname); +} + + +static ImageFileReader::StaticFileType fits_reader{ + "fits", + {"fits", "fit"} +}; + +} diff --git a/SEFramework/src/lib/Image/ImageFileReader.cpp b/SEFramework/src/lib/Image/ImageFileReader.cpp new file mode 100644 index 000000000..338bd1037 --- /dev/null +++ b/SEFramework/src/lib/Image/ImageFileReader.cpp @@ -0,0 +1,198 @@ +/** Copyright © 2019 Université de Genève, LMU Munich - Faculty of Physics, IAP-CNRS/Sorbonne Université + * + * This library is free software; you can redistribute it and/or modify it under + * the terms of the GNU Lesser General Public License as published by the Free + * Software Foundation; either version 3.0 of the License, or (at your option) + * any later version. + * + * This library is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + */ +/** + * @file ImageFileReader.cpp + * @date Aug 14, 2025 + * @author embray + */ + +#include + +#include +#include +#include + +#include + +#include "SEFramework/Image/ImageFileReader.h" + +namespace SourceXtractor { + +// We use this idiom to make sure the map is initialized on the first need +// Otherwise, this may be initialized *after* the concrete implementations try to register themselves +static std::map& getFileTypeMap() { + static std::map file_type_map; + return file_type_map; +} + + +static std::map& getExtensionMap() { + static std::map extension_map; + return extension_map; +} + + +ImageFileReader::ImageFileReader(const std::string& filename) { + std::string base_filename = filename; + std::optional image_path; + + // Check for optional [image_path] + boost::regex suffix_regex(R"((^.*)\[(.*)\]$)"); + boost::smatch match; + if (boost::regex_match(filename, match, suffix_regex)) { + base_filename = match[1]; + image_path = match[2]; + } + + if (image_path) { + m_filename = base_filename; + + try { + size_t idx; + int image_index = std::stoi(*image_path, &idx); + if (idx == image_path->size()) { + m_image_path = std::nullopt; + m_image_index = image_index; + return; + } + } catch (...) { + // Not an integer, treat as string + } + + m_image_path = image_path; + m_image_index = -1; + } else { + m_filename = base_filename; + m_image_path = std::nullopt; + m_image_index = -1; + } +} + + +ImageFileReader::ImageFileReader(const std::string& filename, const std::string& image_path) + : m_filename(filename), m_image_path(image_path), m_image_index(-1) {} + +ImageFileReader::ImageFileReader(const std::string& filename, int image_index) + : m_filename(filename), m_image_index(image_index) {} + + +std::shared_ptr ImageFileReader::get() { + if (m_image_index >= 0) { + return get(m_image_index); + } else if (m_image_path) { + return get(*m_image_path); + } + + return get(0); +} + + +void ImageFileReader::registerFileType(ImageFileType file_type) { + auto lower_name = boost::algorithm::to_lower_copy(file_type.name); + auto& type_map = getFileTypeMap(); + if (type_map.find(lower_name) != type_map.end()) { + throw std::runtime_error("File type already registered: " + file_type.name); + } + + type_map.emplace(lower_name, file_type); + + auto& ext_map = getExtensionMap(); + // Also populate extension lookup map + for (auto& extension : file_type.extensions) { + auto lower_ext = boost::algorithm::to_lower_copy(extension); + if (ext_map.find(lower_ext) != ext_map.end()) { + throw std::runtime_error("File extension already registered: " + extension); + } + ext_map.emplace(extension, file_type); + } +} + + +std::vector ImageFileReader::getFileTypes() { + std::vector keys; + for (auto &e : getFileTypeMap()) { + keys.emplace_back(e.first); + } + return keys; +} + + +std::string ImageFileReader::getDefault() { + return "fits"; +} + + +static std::unique_ptr tryCreateFromFile( + const std::string& base_filename, const std::string& filename, + ImageFileReader::ImageFileType& file_type) { + std::ifstream ifs(base_filename, std::ios::binary); + if (ifs && file_type.test(ifs)) { + return file_type.factory(filename); + } + return nullptr; +} + + +std::unique_ptr ImageFileReader::create(const std::string& filename) { + // Get file extension without leading dot + std::string base_filename = filename; + + // Check for optional [image_path] + boost::regex suffix_regex(R"((^.*)\[(.*)\]$)"); + boost::smatch match; + if (boost::regex_match(filename, match, suffix_regex)) { + base_filename = match[1]; + } + + std::string ext = boost::filesystem::path(base_filename).extension().string(); + if (!ext.empty() && ext[0] == '.') { + ext.erase(0, 1); + } + std::string ext_lower = boost::algorithm::to_lower_copy(ext); + + // Try match by filename extension first + auto& ext_map = getExtensionMap(); + auto entry = ext_map.find(ext_lower); + if (entry != ext_map.end()) { + if (auto img = tryCreateFromFile(base_filename, filename, entry->second)) { + return img; + } + } + + auto& type_map = getFileTypeMap(); + + // Try the default type + std::string default_type = getDefault(); + entry = type_map.find(default_type); + if (entry != ext_map.end()) { + if (auto img = tryCreateFromFile(base_filename, filename, entry->second)) { + return img; + } + } + + // Try all known file types + for (auto& it : type_map) { + if (auto img = tryCreateFromFile(base_filename, filename, it.second)) { + return img; + } + } + + throw Elements::Exception() << "File type of " << base_filename << " could not be determined"; +} + +} // end namespace SourceXtractor + diff --git a/SEFramework/tests/src/FITS/FitsReader_test.cpp b/SEFramework/tests/src/FITS/FitsReader_test.cpp index f667773ec..e1c59ed08 100644 --- a/SEFramework/tests/src/FITS/FitsReader_test.cpp +++ b/SEFramework/tests/src/FITS/FitsReader_test.cpp @@ -23,10 +23,12 @@ #include #include +#include #include #include #include "SEFramework/FITS/FitsReader.h" +#include "SEFramework/Image/ImageFileReader.h" #include "1px.fits.h" @@ -49,7 +51,7 @@ BOOST_AUTO_TEST_SUITE (FitsReader_test) //----------------------------------------------------------------------------- BOOST_FIXTURE_TEST_CASE( read_file, FitsReaderFixture ) { - auto img = FitsReader::readFile(m_tmp_fits.path().native()); + auto img = FitsReader::readFile(m_tmp_fits.path().native()); BOOST_CHECK_EQUAL(img->getWidth(), 1); BOOST_CHECK_EQUAL(img->getHeight(), 1); BOOST_CHECK_EQUAL(img->getChunk(0, 0, 1, 1)->getValue(0, 0), 42); @@ -64,10 +66,29 @@ BOOST_FIXTURE_TEST_CASE ( image_source, FitsReaderFixture ) { BOOST_CHECK_EQUAL(naxis, 2); } + +BOOST_AUTO_TEST_CASE ( detect_file_type ) { + auto reader = ImageFileReader::create(Elements::getAuxiliaryPath("with_primary.fits").native()); + auto* fits_reader = dynamic_cast(reader.get()); + BOOST_CHECK(fits_reader != nullptr); +} + + +BOOST_AUTO_TEST_CASE ( open_with_extname ) { + auto reader = ImageFileReader::create( + Elements::getAuxiliaryPath("multiple_hdu.fits").native() + "[IMAGE2]"); + auto* fits_reader = dynamic_cast(reader.get()); + BOOST_CHECK(fits_reader != nullptr); + auto img = reader->get(); + BOOST_CHECK_EQUAL(img->getWidth(), 1); + BOOST_CHECK_EQUAL(img->getHeight(), 1); +} + + //----------------------------------------------------------------------------- BOOST_AUTO_TEST_CASE( missing_file ) { - BOOST_CHECK_THROW(FitsReader::readFile("/not/existing/path"), Elements::Exception); + BOOST_CHECK_THROW(FitsReader::readFile("/not/existing/path"), Elements::Exception); } //----------------------------------------------------------------------------- diff --git a/SEImplementation/src/lib/Configuration/SegmentationConfig.cpp b/SEImplementation/src/lib/Configuration/SegmentationConfig.cpp index 04a6438f9..4eb6a4338 100644 --- a/SEImplementation/src/lib/Configuration/SegmentationConfig.cpp +++ b/SEImplementation/src/lib/Configuration/SegmentationConfig.cpp @@ -136,6 +136,8 @@ std::shared_ptr SegmentationConfig::getDefault } std::shared_ptr SegmentationConfig::loadFilter(const std::string& filename) const { + // TODO: Replace this to use ImageFileReader, don't need to assume the convolution kernel is + // in a FITS file. // check for the extension ".fits" std::string fits_ending(".fits"); if (filename.length() >= fits_ending.length() @@ -152,7 +154,7 @@ std::shared_ptr SegmentationConfig::loadFilter std::shared_ptr SegmentationConfig::loadFITSFilter(const std::string& filename) const { // read in the FITS file - auto convolution_kernel = FitsReader::readFile(filename); + auto convolution_kernel = FitsReader::readFile(filename); // give some feedback on the filter segConfigLogger.info() << "Loaded segmentation filter: " << filename << " height: " << convolution_kernel->getHeight() << " width: " << convolution_kernel->getWidth(); From fe4a03e144eea07ae736f321e2893a61c938e73c Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Mon, 1 Sep 2025 15:33:57 +0200 Subject: [PATCH 04/29] Add a WCS::identity static constructor This will be useful for creating a DetectionImage from some non-FITS file type; mostly as a placeholder until reading a WCS from ASDF is implemented. --- .../SEFramework/CoordinateSystem/WCS.h | 6 +++++ SEFramework/src/lib/CoordinateSystem/WCS.cpp | 24 +++++++++++++++++++ 2 files changed, 30 insertions(+) diff --git a/SEFramework/SEFramework/CoordinateSystem/WCS.h b/SEFramework/SEFramework/CoordinateSystem/WCS.h index 16f1ac8f7..5f0524e42 100644 --- a/SEFramework/SEFramework/CoordinateSystem/WCS.h +++ b/SEFramework/SEFramework/CoordinateSystem/WCS.h @@ -41,6 +41,9 @@ class WCS : public CoordinateSystem { virtual ~WCS(); + // Create a trivial WCS for a given number of axes + static WCS identity(int naxis); + WorldCoordinate imageToWorld(ImageCoordinate image_coordinate) const override; ImageCoordinate worldToImage(WorldCoordinate world_coordinate) const override; @@ -49,6 +52,9 @@ class WCS : public CoordinateSystem { void addOffset(PixelCoordinate pc); private: + explicit WCS(std::unique_ptr> wcs) + : m_wcs(std::move(wcs)) {} + void init(char* headers, int number_of_records); std::unique_ptr> m_wcs; diff --git a/SEFramework/src/lib/CoordinateSystem/WCS.cpp b/SEFramework/src/lib/CoordinateSystem/WCS.cpp index ac14c0612..fd5997701 100644 --- a/SEFramework/src/lib/CoordinateSystem/WCS.cpp +++ b/SEFramework/src/lib/CoordinateSystem/WCS.cpp @@ -219,6 +219,30 @@ void WCS::init(char* headers, int number_of_records) { WCS::~WCS() { } + +WCS WCS::identity(int naxis) { + auto wcs = std::unique_ptr>( + new wcsprm, [](wcsprm* w) { + if (w) { + wcsfree(w); + delete w; + } + } + ); + + wcsini(1, naxis, wcs.get()); + wcs->flag = -1; + for (int i = 0; i < naxis; i++) { + wcs->crpix[i] = 1.0; + wcs->crval[i] = 0.0; + wcs->cdelt[i] = 1.0; + std::strncpy(wcs->ctype[i], "LINEAR", 72); + } + wcsset(wcs.get()); + return WCS(std::move(wcs)); +} + + WorldCoordinate WCS::imageToWorld(ImageCoordinate image_coordinate) const { // wcsprm is in/out wcsprm wcs_copy; From 8641432a5b2939b8ad6915157ab80138c2d6f9fb Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Mon, 1 Sep 2025 16:36:59 +0200 Subject: [PATCH 05/29] Refactor DetectionImageConfig to use ImageFileReader... ...and prepare for the possibility of non-FITS images. --- .../Configuration/DetectionImageConfig.h | 21 ++- .../Configuration/DetectionImageConfig.cpp | 178 ++++++++++-------- 2 files changed, 113 insertions(+), 86 deletions(-) diff --git a/SEImplementation/SEImplementation/Configuration/DetectionImageConfig.h b/SEImplementation/SEImplementation/Configuration/DetectionImageConfig.h index 54241f27d..90efe9cfa 100644 --- a/SEImplementation/SEImplementation/Configuration/DetectionImageConfig.h +++ b/SEImplementation/SEImplementation/Configuration/DetectionImageConfig.h @@ -24,8 +24,9 @@ #define _SEIMPLEMENTATION_DETECTIONIMAGECONFIG_H #include "Configuration/Configuration.h" +#include "SEFramework/FITS/FitsImageSource.h" #include "SEFramework/Image/Image.h" -#include "SEFramework/Image/ImageSourceWithMetadata.h" +#include "SEFramework/Image/ImageSource.h" #include "SEFramework/CoordinateSystem/CoordinateSystem.h" namespace SourceXtractor { @@ -47,7 +48,7 @@ class DetectionImageConfig : public Euclid::Configuration::Configuration { explicit DetectionImageConfig(long manager_id); std::map getProgramOptions() override; - + void initialize(const UserValues& args) override; std::string getDetectionImagePath() const; @@ -55,7 +56,7 @@ class DetectionImageConfig : public Euclid::Configuration::Configuration { std::shared_ptr getDetectionImage(size_t index = 0) const; std::shared_ptr getCoordinateSystem(size_t index = 0) const; bool isReferenceImage() const { return m_is_reference_image; } - + double getGain(size_t index = 0) const { return m_extensions.at(index).m_gain; } double getSaturation(size_t index = 0) const { return m_extensions.at(index).m_saturation; } int getInterpolationGap(size_t index = 0) const { return m_extensions.at(index).m_interpolation_gap; } @@ -90,6 +91,20 @@ class DetectionImageConfig : public Euclid::Configuration::Configuration { double m_flux_scale {1.0}; int m_interpolation_gap {0}; + + DetectionImageExtension() = default; + + DetectionImageExtension(std::shared_ptr image_source, double gain, + double saturation, double flux_scale, int interpolation_gap); + DetectionImageExtension(std::shared_ptr fits_image_source, double gain, + double saturation, double flux_scale, int interpolation_gap); + + static DetectionImageExtension create(std::shared_ptr, const UserValues& args); + + private: + void init(std::shared_ptr image_source, double gain, double saturation, + double flux_scale, int interpolation_gap); + void rescale(); }; std::vector m_extensions; diff --git a/SEImplementation/src/lib/Configuration/DetectionImageConfig.cpp b/SEImplementation/src/lib/Configuration/DetectionImageConfig.cpp index de46e9678..8f0d6ed9f 100644 --- a/SEImplementation/src/lib/Configuration/DetectionImageConfig.cpp +++ b/SEImplementation/src/lib/Configuration/DetectionImageConfig.cpp @@ -26,8 +26,9 @@ using boost::regex; using boost::regex_match; using boost::smatch; -#include -#include +#include "SEFramework/Image/BufferedImage.h" +#include "SEFramework/Image/ImageFileReader.h" +#include "SEFramework/Image/ProcessedImage.h" #include "SEFramework/FITS/FitsImageSource.h" #include "SEFramework/CoordinateSystem/WCS.h" @@ -95,109 +96,120 @@ void DetectionImageConfig::initialize(const UserValues& args) { m_detection_image_path = args.find(DETECTION_IMAGE)->second.as(); - boost::regex hdu_regex(".*\\[[0-9]*\\]$"); + auto image_reader = ImageFileReader::create(m_detection_image_path); - for (int i=0;; i++) { - DetectionImageExtension extension; + for (const auto& img_source: *image_reader) { + DetectionImageExtension extension = DetectionImageExtension::create(img_source, args); + m_extensions.emplace_back(std::move(extension)); + } +} - std::shared_ptr fits_image_source; - if (boost::regex_match(m_detection_image_path, hdu_regex)) { - if (i==0) { - fits_image_source = std::make_shared(m_detection_image_path, 0, ImageTile::FloatImage); - } else { - break; - } - } else { - try { - fits_image_source = std::make_shared(m_detection_image_path, i+1, ImageTile::FloatImage); - } catch (...) { - if (i==0) { - // Skip past primary HDU if it doesn't have an image - continue; - } else { - if (m_extensions.size() == 0) { - throw; - } - break; - } - } - } - extension.m_image_source = fits_image_source; - extension.m_detection_image = BufferedImage::create(fits_image_source); - extension.m_coordinate_system = std::make_shared(*fits_image_source); - - double detection_image_gain = 0, detection_image_saturate = 0; - auto img_metadata = fits_image_source->getMetadata(); - - if (img_metadata.count("GAIN")){ - // read the keyword GAIN from the metadata - if (double* double_gain = boost::get(&img_metadata.at("GAIN").m_value)){ - detection_image_gain = *double_gain; - } else if (int64_t *int64_gain = boost::get(&img_metadata.at("GAIN").m_value)){ - detection_image_gain = (double) *int64_gain; - } - else { - throw Elements::Exception() << "Keyword GAIN must be either float or int!"; - } - } +DetectionImageConfig::DetectionImageExtension::DetectionImageExtension( + std::shared_ptr image_source, double gain, double saturation, + double flux_scale, int interpolation_gap) { + init(image_source, gain, saturation, flux_scale, interpolation_gap); + m_coordinate_system = std::make_shared(WCS::identity(2)); + rescale(); +} - if (img_metadata.count("SATURATE")){ - // read the keyword SATURATE from the metadata - if (double* double_saturate = boost::get(&img_metadata.at("SATURATE").m_value)){ - detection_image_saturate = *double_saturate; - } else if (int64_t *int64_saturate = boost::get(&img_metadata.at("SATURATE").m_value)){ - detection_image_saturate = (double) *int64_saturate; - } - else { - throw Elements::Exception() << "Keyword SATURATE must be either float or int!"; - } - } - if (args.find(DETECTION_IMAGE_FLUX_SCALE) != args.end()) { - extension.m_flux_scale = args.find(DETECTION_IMAGE_FLUX_SCALE)->second.as(); - } - else if (img_metadata.count("FLXSCALE")) { - // read the keyword FLXSCALE from the metadata - if (double* f_scale = boost::get(&img_metadata.at("FLXSCALE").m_value)){ - extension.m_flux_scale = *f_scale; - } else if (int64_t *int64_f_scale = boost::get(&img_metadata.at("FLXSCALE").m_value)){ - extension.m_flux_scale = (double) *int64_f_scale; +DetectionImageConfig::DetectionImageExtension::DetectionImageExtension( + std::shared_ptr fits_image_source, double gain, double saturation, + double flux_scale, int interpolation_gap) { + init(fits_image_source, gain, saturation, flux_scale, interpolation_gap); + m_coordinate_system = std::make_shared(*fits_image_source); + auto img_metadata = fits_image_source->getMetadata(); + + if (img_metadata.count("GAIN")){ + // read the keyword GAIN from the metadata + if (double* double_gain = boost::get(&img_metadata.at("GAIN").m_value)){ + m_gain = *double_gain; + } else if (int64_t *int64_gain = boost::get(&img_metadata.at("GAIN").m_value)){ + m_gain = (double) *int64_gain; } else { - throw Elements::Exception() << "Keyword FLXSCALE must be either float or int!"; + throw Elements::Exception() << "Keyword GAIN must be either float or int!"; } - } + } - if (args.find(DETECTION_IMAGE_GAIN) != args.end()) { - extension.m_gain = args.find(DETECTION_IMAGE_GAIN)->second.as(); + if (img_metadata.count("FLXSCALE")) { + // read the keyword FLXSCALE from the metadata + if (double* f_scale = boost::get(&img_metadata.at("FLXSCALE").m_value)){ + m_flux_scale = *f_scale; + } else if (int64_t *int64_f_scale = boost::get(&img_metadata.at("FLXSCALE").m_value)){ + m_flux_scale = (double) *int64_f_scale; } else { - extension.m_gain = detection_image_gain; + throw Elements::Exception() << "Keyword FLXSCALE must be either float or int!"; } + } - if (args.find(DETECTION_IMAGE_SATURATION) != args.end()) { - extension.m_saturation = args.find(DETECTION_IMAGE_SATURATION)->second.as(); + if (img_metadata.count("SATURATE")){ + // read the keyword SATURATE from the metadata + if (double* double_saturate = boost::get(&img_metadata.at("SATURATE").m_value)){ + m_saturation = *double_saturate; + } else if (int64_t *int64_saturate = boost::get(&img_metadata.at("SATURATE").m_value)){ + m_saturation = (double) *int64_saturate; } else { - extension.m_saturation = detection_image_saturate; + throw Elements::Exception() << "Keyword SATURATE must be either float or int!"; } + } - extension.m_interpolation_gap = args.find(DETECTION_IMAGE_INTERPOLATION)->second.as() ? - std::max(0, args.find(DETECTION_IMAGE_INTERPOLATION_GAP)->second.as()) : 0; + rescale(); +} - // Adapt image and parameters to take flux_scale into consideration - if (extension.m_flux_scale != 1.0) { - extension.m_detection_image = - MultiplyImage::create(extension.m_detection_image, extension.m_flux_scale); - extension.m_gain /= extension.m_flux_scale; - extension.m_saturation *= extension.m_flux_scale; - } - m_extensions.emplace_back(std::move(extension)); +void DetectionImageConfig::DetectionImageExtension::init( + std::shared_ptr image_source, double gain, double saturation, + double flux_scale, int interpolation_gap) { + m_image_source = image_source; + m_detection_image = BufferedImage::create(image_source); + m_gain = gain; + m_saturation = saturation; + m_flux_scale = flux_scale; + m_interpolation_gap = interpolation_gap; +} + + +DetectionImageConfig::DetectionImageExtension DetectionImageConfig::DetectionImageExtension::create( + std::shared_ptr image_source, const UserValues& args) { + double gain = (args.find(DETECTION_IMAGE_GAIN) != args.end()) + ? args.find(DETECTION_IMAGE_GAIN)->second.as() + : 0.0; + + double saturation = (args.find(DETECTION_IMAGE_SATURATION) != args.end()) + ? args.find(DETECTION_IMAGE_SATURATION)->second.as() + : 0.0; + + double flux_scale = (args.find(DETECTION_IMAGE_FLUX_SCALE) != args.end()) + ? args.find(DETECTION_IMAGE_FLUX_SCALE)->second.as() + : 1.0; + + int interpolation_gap = args.find(DETECTION_IMAGE_INTERPOLATION)->second.as() + ? std::max(0, args.find(DETECTION_IMAGE_INTERPOLATION_GAP)->second.as()) + : 0; + + if (auto fits_src = std::dynamic_pointer_cast(image_source)) { + return DetectionImageExtension(fits_src, gain, saturation, flux_scale, interpolation_gap); + } + + return DetectionImageExtension(image_source, gain, saturation, flux_scale, interpolation_gap); +} + + +void DetectionImageConfig::DetectionImageExtension::rescale() { + // Adapt image and parameters to take flux_scale into consideration + if (m_flux_scale != 1.0) { + m_detection_image = + MultiplyImage::create(m_detection_image, m_flux_scale); + m_gain /= m_flux_scale; + m_saturation *= m_flux_scale; } } + std::string DetectionImageConfig::getDetectionImagePath() const { return m_detection_image_path; } From 77cc3863d8e1230b7209ab9c9b4a6cb56edeaff5 Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Tue, 2 Sep 2025 12:46:19 +0200 Subject: [PATCH 06/29] Fix: make FitsReader iteration more sane--it can iterate over all image HDUs (even skipping over non-image HDUs) Consequently make FitsReader::get(int) more sane: it returns the N-th image HDU in the file (0-indexed) rather than FITSish 1-indexing. This enables more consistent use of ImageFileReader::get regardless of the underlying file type. --- .../SEFramework/FITS/FitsImageSource.h | 13 +++++ SEFramework/SEFramework/FITS/FitsReader.h | 18 +++++- .../SEFramework/Image/ImageFileReader.h | 10 +++- SEFramework/src/lib/FITS/FitsImageSource.cpp | 16 +++--- SEFramework/src/lib/FITS/FitsReader.cpp | 55 +++++++++++++------ .../tests/src/FITS/FitsReader_test.cpp | 17 +++++- 6 files changed, 102 insertions(+), 27 deletions(-) diff --git a/SEFramework/SEFramework/FITS/FitsImageSource.h b/SEFramework/SEFramework/FITS/FitsImageSource.h index bedf630f1..042d567e6 100644 --- a/SEFramework/SEFramework/FITS/FitsImageSource.h +++ b/SEFramework/SEFramework/FITS/FitsImageSource.h @@ -31,6 +31,8 @@ #include +#include + #include "FilePool/FileManager.h" #include "SEFramework/CoordinateSystem/CoordinateSystem.h" #include "SEFramework/Image/ImageSourceWithMetadata.h" @@ -45,9 +47,20 @@ using Euclid::FilePool::FileHandler; class FitsImageSource : public ImageSource, public std::enable_shared_from_this { public: + /** + * Exception thrown when opening the file with or switching to an HDU that does not + * exist. + */ + class UnknownHduException : public Elements::Exception {}; static std::shared_ptr getFileManager(unsigned int max_open_files = 500); + /** + * Exception thrown when opening the file with or switching to an HDU that does not contain + * an image array. + */ + class InvalidHduTypeException : public Elements::Exception {}; + /** * Constructor * @param filename diff --git a/SEFramework/SEFramework/FITS/FitsReader.h b/SEFramework/SEFramework/FITS/FitsReader.h index f6195c5be..6f0128867 100644 --- a/SEFramework/SEFramework/FITS/FitsReader.h +++ b/SEFramework/SEFramework/FITS/FitsReader.h @@ -40,9 +40,22 @@ class FitsReader : public ImageFileReader { using ImageFileReader::ImageFileReader; static bool test(std::istream& stream); - std::shared_ptr get() override; + + using ImageFileReader::get; + /** + * Get the N-th image HDU from the FITS file + * + * NOTE: For consistency's sake across ImageFileReaders this is 0-indexed, and *not* the + * 1-indexed FORTRAN/FITS convention. It also does not correspond to the HDU number itself, + * but rather the count of HDUs containing valid image arrays. + * + * Use FitsReader.getHdu to get by FITS HDU number + */ std::shared_ptr get(int image_index) override; std::shared_ptr get(const std::string& image_path) override; + + std::shared_ptr getHdu(int hdu_num); + // TODO: Integrate this interface into the ImageFileReader base class template static std::shared_ptr> readFile(const std::string& filename) { @@ -50,6 +63,9 @@ class FitsReader : public ImageFileReader { return BufferedImage::create(image_source); } +private: + std::map m_image_hdu_map; + }; /* End of FitsReader class */ } /* namespace SourceXtractor */ diff --git a/SEFramework/SEFramework/Image/ImageFileReader.h b/SEFramework/SEFramework/Image/ImageFileReader.h index f1fa2713c..abce18f95 100644 --- a/SEFramework/SEFramework/Image/ImageFileReader.h +++ b/SEFramework/SEFramework/Image/ImageFileReader.h @@ -75,6 +75,8 @@ class ImageFileReader { using pointer = std::shared_ptr*; using reference = std::shared_ptr&; + Iterator() noexcept : m_reader(nullptr), m_index(0) {} + Iterator(ImageFileReader* reader, int index) : m_reader(reader), m_index(index) { advance(); @@ -98,6 +100,12 @@ class ImageFileReader { private: void advance() { + // If m_reader is null the iterator is considered "done" + if (nullptr == m_reader) { + m_current = nullptr; + return; + } + try { m_current = m_reader->get(m_index); } catch (...) { @@ -115,7 +123,7 @@ class ImageFileReader { } Iterator end() { - return Iterator(this, -1); + return Iterator(); } /** diff --git a/SEFramework/src/lib/FITS/FitsImageSource.cpp b/SEFramework/src/lib/FITS/FitsImageSource.cpp index 72bd4a597..ea3586a1e 100644 --- a/SEFramework/src/lib/FITS/FitsImageSource.cpp +++ b/SEFramework/src/lib/FITS/FitsImageSource.cpp @@ -69,7 +69,7 @@ ImageTile::ImageType convertImageType(int bitpix) { image_type = ImageTile::LongLongImage; break; default: - throw Elements::Exception() << "Unsupported FITS image type: " << bitpix; + throw FitsImageSource::InvalidHduTypeException() << "Unsupported FITS image type: " << bitpix; } return image_type; @@ -94,12 +94,12 @@ FitsImageSource::FitsImageSource(const std::string& filename, int hdu_number, if (m_hdu_number <= 0 && !extname) { if (fits_get_hdu_num(fptr, &m_hdu_number) < 0) { - throw Elements::Exception() << "Can't get the active HDU from the FITS file: " << filename; + throw UnknownHduException() << "Can't get the active HDU from the FITS file: " << filename; } } else if (extname) { switchHdu(fptr, *extname); if (fits_get_hdu_num(fptr, &m_hdu_number) < 0) { - throw Elements::Exception() << "Can't get the active HDU from the FITS file: " << filename; + throw UnknownHduException() << "Can't get the active HDU from the FITS file: " << filename; } } else { switchHdu(fptr, m_hdu_number); @@ -109,7 +109,7 @@ FitsImageSource::FitsImageSource(const std::string& filename, int hdu_number, if (status != 0 || (naxis != 2 && naxis != 3)) { char error_message[32]; fits_get_errstatus(status, error_message); - throw Elements::Exception() + throw InvalidHduTypeException() << "Can't find 2D image or data cube in FITS file: " << filename << "[" << m_hdu_number << "]" << " status: " << status << " = " << error_message; } @@ -166,7 +166,7 @@ FitsImageSource::FitsImageSource(const std::string& filename, int width, int hei if (fits_get_hdu_num(fptr, &m_hdu_number) < 0) { char error_message[32]; fits_get_errstatus(status, error_message); - throw Elements::Exception() << "Can't get the active HDU from the FITS file: " << filename + throw UnknownHduException() << "Can't get the active HDU from the FITS file: " << filename << " status: " << status << " = " << error_message; } @@ -275,11 +275,11 @@ void FitsImageSource::switchHdu(fitsfile *fptr, int hdu_number) const { if (status != 0) { char error_message[32]; fits_get_errstatus(status, error_message); - throw Elements::Exception() << "Could not switch to HDU # " << hdu_number << " in file " << m_filename + throw UnknownHduException() << "Could not switch to HDU # " << hdu_number << " in file " << m_filename << " status: " << status << " = " << error_message; } if (hdu_type != IMAGE_HDU) { - throw Elements::Exception() << "Trying to access non-image HDU in file " << m_filename; + throw InvalidHduTypeException() << "Trying to access non-image HDU in file " << m_filename; } } @@ -292,7 +292,7 @@ void FitsImageSource::switchHdu(fitsfile *fptr, const std::string& extname) cons if (status != 0) { char error_message[32]; fits_get_errstatus(status, error_message); - throw Elements::Exception() << "Could not switch to HDU with EXTNAME " << extname << " in file " << m_filename + throw UnknownHduException() << "Could not switch to HDU with EXTNAME " << extname << " in file " << m_filename << " status: " << status << " = " << error_message; } } diff --git a/SEFramework/src/lib/FITS/FitsReader.cpp b/SEFramework/src/lib/FITS/FitsReader.cpp index 5d6a1230d..f07c0c699 100644 --- a/SEFramework/src/lib/FITS/FitsReader.cpp +++ b/SEFramework/src/lib/FITS/FitsReader.cpp @@ -31,36 +31,59 @@ namespace SourceXtractor { bool FitsReader::test(std::istream& stream) { - char header[9] = {0}; - if (!stream.read(header, 9)) { - return false; - } + char header[9] = {0}; + if (!stream.read(header, 9)) { + return false; + } - return std::strncmp(header, "SIMPLE =", 9) == 0; + return std::strncmp(header, "SIMPLE =", 9) == 0; } -std::shared_ptr FitsReader::get() { - if (0 == m_image_index) { - // If the primary HDU was not an image try the first extension HDU +std::shared_ptr FitsReader::get(int image_index) { + auto it = m_image_hdu_map.find(image_index); + + if (it != m_image_hdu_map.end()) { + return std::make_shared(m_filename, it->second); + } + + // Try loading HDUs from the file until we find the next image HDU + int hdu_num = 0; + int known_index = -1; + + if (!m_image_hdu_map.empty()) { + auto rit = m_image_hdu_map.rbegin(); + known_index = rit->first; + hdu_num = rit->second; + } + + while (known_index < image_index) { try { - return get(1); - } catch (...) { - return get(2); + auto image_source = getHdu(++hdu_num); + m_image_hdu_map[++known_index] = hdu_num; + + if (known_index == image_index) { + return image_source; + } + } catch (const FitsImageSource::InvalidHduTypeException&) { + continue; + } catch (const FitsImageSource::UnknownHduException&) { + throw; } } - return ImageFileReader::get(); + // Shouldn't get here + throw FitsImageSource::UnknownHduException(); } -std::shared_ptr FitsReader::get(int hdu_num) { - return std::make_shared(m_filename, hdu_num); +std::shared_ptr FitsReader::get(const std::string& extname) { + return std::make_shared(m_filename, extname); } -std::shared_ptr FitsReader::get(const std::string& extname) { - return std::make_shared(m_filename, extname); +std::shared_ptr FitsReader::getHdu(int hdu_num) { + return std::make_shared(m_filename, hdu_num); } diff --git a/SEFramework/tests/src/FITS/FitsReader_test.cpp b/SEFramework/tests/src/FITS/FitsReader_test.cpp index e1c59ed08..56c0e5c68 100644 --- a/SEFramework/tests/src/FITS/FitsReader_test.cpp +++ b/SEFramework/tests/src/FITS/FitsReader_test.cpp @@ -66,6 +66,7 @@ BOOST_FIXTURE_TEST_CASE ( image_source, FitsReaderFixture ) { BOOST_CHECK_EQUAL(naxis, 2); } +//----------------------------------------------------------------------------- BOOST_AUTO_TEST_CASE ( detect_file_type ) { auto reader = ImageFileReader::create(Elements::getAuxiliaryPath("with_primary.fits").native()); @@ -73,6 +74,21 @@ BOOST_AUTO_TEST_CASE ( detect_file_type ) { BOOST_CHECK(fits_reader != nullptr); } +//----------------------------------------------------------------------------- + +BOOST_AUTO_TEST_CASE ( iterate ) { + auto reader = ImageFileReader::create(Elements::getAuxiliaryPath("with_primary.fits").native()); + // Should just return one image, from the primary HDU (the next HDU in this file is a table) + int image_count = 0; + for (const auto& img_source: *reader) { + auto fits_source = std::dynamic_pointer_cast(img_source); + BOOST_CHECK(fits_source != nullptr); + image_count++; + } + BOOST_CHECK_EQUAL(image_count, 1); +} + +//----------------------------------------------------------------------------- BOOST_AUTO_TEST_CASE ( open_with_extname ) { auto reader = ImageFileReader::create( @@ -84,7 +100,6 @@ BOOST_AUTO_TEST_CASE ( open_with_extname ) { BOOST_CHECK_EQUAL(img->getHeight(), 1); } - //----------------------------------------------------------------------------- BOOST_AUTO_TEST_CASE( missing_file ) { From d30f61fe97994229e04d4268fa9a84f64c063b87 Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Wed, 3 Sep 2025 11:54:05 +0200 Subject: [PATCH 07/29] TODO: Fix the redundant code mentioned in this comment later --- SEFramework/src/lib/FITS/FitsReader.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/SEFramework/src/lib/FITS/FitsReader.cpp b/SEFramework/src/lib/FITS/FitsReader.cpp index f07c0c699..3486fec0c 100644 --- a/SEFramework/src/lib/FITS/FitsReader.cpp +++ b/SEFramework/src/lib/FITS/FitsReader.cpp @@ -41,6 +41,10 @@ bool FitsReader::test(std::istream& stream) { std::shared_ptr FitsReader::get(int image_index) { + // TODO: It turns out none of this is necessary: It's already implemented in + // the FitsFile class, which actually loops over all the HDUs when opening the file + // and gets the HDU numbers of images (can be accessed with FitsFile::getImageHdus) + // so should just use that instead. auto it = m_image_hdu_map.find(image_index); if (it != m_image_hdu_map.end()) { From 01f8cb06dfcb624ecc1a5e298989a75e4a04eb04 Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Fri, 5 Sep 2025 13:20:16 +0200 Subject: [PATCH 08/29] Add some new capabilities to AsdfFile to allow getting some primitive Ndarray objects out of it. There is a bit of a foot-gun here that AsdfFile has to be wrapped in shared_ptr when using this. I may try to get rid of that limitation later. --- SEFramework/SEFramework/ASDF/AsdfFile.h | 88 +++++++++++++++++++- SEFramework/src/lib/ASDF/AsdfFile.cpp | 79 +++++++++++++++++- SEFramework/tests/src/ASDF/AsdfFile_test.cpp | 28 +++++++ 3 files changed, 191 insertions(+), 4 deletions(-) diff --git a/SEFramework/SEFramework/ASDF/AsdfFile.h b/SEFramework/SEFramework/ASDF/AsdfFile.h index 3546b97f9..77b0a556b 100644 --- a/SEFramework/SEFramework/ASDF/AsdfFile.h +++ b/SEFramework/SEFramework/ASDF/AsdfFile.h @@ -25,10 +25,15 @@ #ifndef _SEFRAMEWORK_ASDF_ASDFFILE_H_ #define _SEFRAMEWORK_ASDF_ASDFFILE_H_ +#include + #include #include +#include + #include "SEFramework/Image/ImageSourceWithMetadata.h" +#include "SEFramework/Image/ImageTile.h" namespace SourceXtractor { @@ -37,9 +42,21 @@ namespace SourceXtractor { * @brief represents access to a whole ASDF file * */ -class AsdfFile { +class AsdfFile : public std::enable_shared_from_this { public: - AsdfFile(const boost::filesystem::path& path); + /** + * Exception thrown when trying to access a value that does not exist in the ASDF tree + */ + class AsdfValueNotFoundException : public Elements::Exception {}; + + /** + * Exception thrown when trying to access a value of the wrong type from the ASDF tree + */ + class AsdfValueTypeMismatchException : public Elements::Exception {}; + + AsdfFile(const boost::filesystem::path& path, bool writeable); + + AsdfFile(const boost::filesystem::path& path) : AsdfFile(path, false) {} AsdfFile(AsdfFile&&) = default; @@ -47,7 +64,74 @@ class AsdfFile { asdf_file_t* getAsdfFilePtr(); + /** + * Wrapper around asdf_ndarray_t + */ + class Ndarray { + public: + friend class AsdfFile; + + ~Ndarray() { + asdf_ndarray_destroy(m_ndarray_ptr); + } + + uint32_t ndim() const { return m_ndarray_ptr->ndim; } + + std::vector shape() const { + return std::vector(m_ndarray_ptr->shape, + m_ndarray_ptr->shape + m_ndarray_ptr->ndim); + } + + /** + * Given a shared pointer to an already allocated ImageTile (i.e. from ImageTile::create) + * copy a tile from the raw ndarray data into the ImageTile. + */ + void fillImageTile(const std::shared_ptr image_tile); + + private: + /** + * Private constructor for creating the `Ndarray` wrapper from a raw asdf_value_t * + */ + explicit Ndarray(std::shared_ptr file, asdf_value_t *ptr); + /** + * Private constructor for creating the `Ndarray` wrapper from a raw asdf_ndarray_t * + */ + explicit Ndarray(std::shared_ptr file, asdf_ndarray_t *ptr) + : m_file(file), m_ndarray_ptr(ptr) {} + + std::shared_ptr m_file; + asdf_ndarray_t* m_ndarray_ptr; + }; + + /** + * Return the N-th ndarray from the top-level of the ASDF tree iterating the top-level + * keys in order. + */ + Ndarray getNdarray(int index); + + /** + * Return any ndarray from any path in the ASDF tree + * + * If the given path does not exist an AsdfValueNotFoundException is thrown. + * If the given path exists but is not an ndarray, an AsdfValueTypeMismatchException is + * thrown. + */ + Ndarray getNdarray(const std::string& path); + + /* TODO: More general methods for reading metadata from the ASDF tree; for the first version + * not needed though. */ + private: + struct AsdfValueDestroy { + void operator()(asdf_value_t *v) const noexcept { + asdf_value_destroy(v); + } + }; + + using AsdfValuePtr = std::unique_ptr; + + AsdfValuePtr getValue(const std::string& path); + boost::filesystem::path m_path; std::unique_ptr m_asdf_ptr; diff --git a/SEFramework/src/lib/ASDF/AsdfFile.cpp b/SEFramework/src/lib/ASDF/AsdfFile.cpp index 5aeccdc26..80af0aeba 100644 --- a/SEFramework/src/lib/ASDF/AsdfFile.cpp +++ b/SEFramework/src/lib/ASDF/AsdfFile.cpp @@ -38,14 +38,20 @@ namespace SourceXtractor { static Elements::Logging logger = Elements::Logging::getLogger("AsdfFile"); -AsdfFile::AsdfFile(const boost::filesystem::path& path) +AsdfFile::AsdfFile(const boost::filesystem::path& path, bool writeable) : m_path(path), m_asdf_ptr(nullptr, asdf_close) { + if (writeable) { + throw Elements::Exception() << "ASDF files are not currently supported for writing"; + } + open(); } + AsdfFile::~AsdfFile() {} + asdf_file_t* AsdfFile::getAsdfFilePtr() { return m_asdf_ptr.get(); } @@ -64,9 +70,78 @@ void AsdfFile::open() { if (error_message != nullptr) { throw Elements::Exception() - << "Can't open ASDF file: " << m_path << " reason: " << error_message; + << "Can't open ASDF file: " << m_path.native() << " reason: " << error_message; + } + } + + m_asdf_ptr.reset(ptr); +} + +// TODO: Support negative indexing as well +AsdfFile::Ndarray AsdfFile::getNdarray(int index) { + AsdfValuePtr root = getValue("/"); + + if (!asdf_value_is_mapping(root.get())) { + throw AsdfValueNotFoundException() << "ASDF tree root is not a mapping in file: " + << m_path.native(); + } + + asdf_mapping_iter_t iter = asdf_mapping_iter_init(); + asdf_mapping_item_t* item; + int count = 0; + + while ((item = asdf_mapping_iter(root.get(), &iter)) != nullptr) { + asdf_value_t* value = asdf_mapping_item_value(item); + if (asdf_value_is_ndarray(value)) { + if (count == index) { + return Ndarray(shared_from_this(), value); + } else if (count > index) { + break; + } + count++; + } + } + + throw AsdfValueNotFoundException() << "No ndarray at index " << index << " in file: " + << m_path.native(); +} + +AsdfFile::Ndarray AsdfFile::getNdarray(const std::string &path) { + AsdfValuePtr value = getValue(path); + return Ndarray(shared_from_this(), value.get()); +} + + +AsdfFile::AsdfValuePtr AsdfFile::getValue(const std::string& path) { + asdf_value_t *v = asdf_get_value(m_asdf_ptr.get(), path.c_str()); + if (!v) { + throw AsdfValueNotFoundException() << "No value at path " << path << " in file: " + << m_path.native(); + } + + return AsdfValuePtr(v); +} + +AsdfFile::Ndarray::Ndarray(std::shared_ptr file, asdf_value_t *value) + : Ndarray(file, (asdf_ndarray_t*)nullptr) { + asdf_ndarray_t *ndarray_ptr = nullptr; + asdf_value_err_t err = asdf_value_as_ndarray(value, &ndarray_ptr); + switch (err) { + case ASDF_VALUE_OK: + // Value exists and is an ndarray: OK + break; + case ASDF_VALUE_ERR_TYPE_MISMATCH: { + const char* path = asdf_value_path(value); + throw AsdfValueTypeMismatchException() << "Value at " << path << " is not an ndarray: " + << file->m_path.native(); + } + default: { + const char *error_message = asdf_error(file->getAsdfFilePtr()); + throw AsdfValueNotFoundException() << "An error occurred reading the ASDF file " + << file->m_path.native() << ": " << error_message; } } + m_ndarray_ptr = ndarray_ptr; } } // namespace SourceXtractor diff --git a/SEFramework/tests/src/ASDF/AsdfFile_test.cpp b/SEFramework/tests/src/ASDF/AsdfFile_test.cpp index acf4456e2..99fa83a12 100644 --- a/SEFramework/tests/src/ASDF/AsdfFile_test.cpp +++ b/SEFramework/tests/src/ASDF/AsdfFile_test.cpp @@ -52,5 +52,33 @@ BOOST_AUTO_TEST_CASE( missing_file ) { //----------------------------------------------------------------------------- +BOOST_AUTO_TEST_CASE( get_ndarray_by_index ) { + auto asdf_file = std::make_shared(Elements::getAuxiliaryPath("basic.asdf").native()); + auto ndarray = asdf_file->getNdarray(0); + BOOST_CHECK_EQUAL(ndarray.ndim(), 1); + auto shape = ndarray.shape(); + std::vector expected_shape{8}; + BOOST_CHECK_EQUAL_COLLECTIONS( + shape.begin(), shape.end(), + expected_shape.begin(), expected_shape.end() + ); +} + +//----------------------------------------------------------------------------- + +BOOST_AUTO_TEST_CASE( get_ndarray_by_path ) { + auto asdf_file = std::make_shared(Elements::getAuxiliaryPath("basic.asdf").native()); + auto ndarray = asdf_file->getNdarray("data"); + BOOST_CHECK_EQUAL(ndarray.ndim(), 1); + auto shape = ndarray.shape(); + std::vector expected_shape{8}; + BOOST_CHECK_EQUAL_COLLECTIONS( + shape.begin(), shape.end(), + expected_shape.begin(), expected_shape.end() + ); +} + +//----------------------------------------------------------------------------- + BOOST_AUTO_TEST_SUITE_END () From e836139eb51218c5b89658f26d11ab4dbe35ff8d Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Wed, 10 Sep 2025 17:35:10 +0200 Subject: [PATCH 09/29] Implement AsdfFile::Ndarray::fileImageTile with asdf_ndarray_read_tile_2d This will be used to implement AsdfImageSource::getImageTile --- SEFramework/SEFramework/ASDF/AsdfFile.h | 2 +- SEFramework/src/lib/ASDF/AsdfFile.cpp | 31 +++++++++++++++++++++++++ 2 files changed, 32 insertions(+), 1 deletion(-) diff --git a/SEFramework/SEFramework/ASDF/AsdfFile.h b/SEFramework/SEFramework/ASDF/AsdfFile.h index 77b0a556b..83550ef95 100644 --- a/SEFramework/SEFramework/ASDF/AsdfFile.h +++ b/SEFramework/SEFramework/ASDF/AsdfFile.h @@ -86,7 +86,7 @@ class AsdfFile : public std::enable_shared_from_this { * Given a shared pointer to an already allocated ImageTile (i.e. from ImageTile::create) * copy a tile from the raw ndarray data into the ImageTile. */ - void fillImageTile(const std::shared_ptr image_tile); + void fillImageTile(const std::shared_ptr image_tile, int layer); private: /** diff --git a/SEFramework/src/lib/ASDF/AsdfFile.cpp b/SEFramework/src/lib/ASDF/AsdfFile.cpp index 80af0aeba..47e6cac00 100644 --- a/SEFramework/src/lib/ASDF/AsdfFile.cpp +++ b/SEFramework/src/lib/ASDF/AsdfFile.cpp @@ -122,6 +122,7 @@ AsdfFile::AsdfValuePtr AsdfFile::getValue(const std::string& path) { return AsdfValuePtr(v); } + AsdfFile::Ndarray::Ndarray(std::shared_ptr file, asdf_value_t *value) : Ndarray(file, (asdf_ndarray_t*)nullptr) { asdf_ndarray_t *ndarray_ptr = nullptr; @@ -144,4 +145,34 @@ AsdfFile::Ndarray::Ndarray(std::shared_ptr file, asdf_value_t *value) m_ndarray_ptr = ndarray_ptr; } + +void AsdfFile::Ndarray::fillImageTile(const std::shared_ptr image_tile, int layer) { + uint64_t plane_origin = layer; + void* data = image_tile->getDataPtr(); + asdf_ndarray_err_t err = asdf_ndarray_read_tile_2d( + m_ndarray_ptr, + image_tile->getPosX(), + image_tile->getPosY(), + image_tile->getWidth(), + image_tile->getHeight(), + &plane_origin, + &data + ); + + switch (err) { + case ASDF_NDARRAY_OK: + break; + case ASDF_NDARRAY_ERR_OUT_OF_BOUNDS: + throw Elements::Exception() << "requested image tile (" << image_tile->getPosX() + << ", " << image_tile->getPosY() << ") -> (" + << image_tile->getPosX() + image_tile->getWidth() << ", " + << image_tile->getPosY() + image_tile->getHeight() << ") is out of bounds"; + case ASDF_NDARRAY_ERR_OOM: + throw Elements::Exception() << "out of memory copying image tile"; + case ASDF_NDARRAY_ERR_INVAL: + default: + throw Elements::Exception() << "invalid argument to asdf_ndarray_read_tile_2d"; + } +} + } // namespace SourceXtractor From 2d37d2e9a02c84918ac4ed06e8c665e55a0e0202 Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Mon, 15 Sep 2025 12:02:29 +0200 Subject: [PATCH 10/29] complete implementation of the basic AsdfImageSource and AsdfReader classes Drop use of enable_shared_from_this on AsdfFile which was misguided. Should clarify that the lifetime of an `AsdfFile::Ndarray` is concurrent with the `AsdfFile` itself (which is kept alive so long as we have a pointer to its `FileAccessor`) --- SEFramework/SEFramework/ASDF/AsdfFile.h | 43 ++++-- .../SEFramework/ASDF/AsdfImageSource.h | 135 ++++++++++++++++++ SEFramework/SEFramework/ASDF/AsdfReader.h | 56 ++++++++ SEFramework/SEFramework/FITS/FitsReader.h | 11 +- SEFramework/src/lib/ASDF/AsdfFile.cpp | 72 ++++++++-- SEFramework/src/lib/ASDF/AsdfImageSource.cpp | 118 +++++++++++++++ SEFramework/src/lib/ASDF/AsdfReader.cpp | 65 +++++++++ SEFramework/src/lib/FITS/FitsReader.cpp | 2 +- SEFramework/tests/src/ASDF/AsdfFile_test.cpp | 16 +-- 9 files changed, 480 insertions(+), 38 deletions(-) create mode 100644 SEFramework/SEFramework/ASDF/AsdfImageSource.h create mode 100644 SEFramework/SEFramework/ASDF/AsdfReader.h create mode 100644 SEFramework/src/lib/ASDF/AsdfImageSource.cpp create mode 100644 SEFramework/src/lib/ASDF/AsdfReader.cpp diff --git a/SEFramework/SEFramework/ASDF/AsdfFile.h b/SEFramework/SEFramework/ASDF/AsdfFile.h index 83550ef95..b7594922c 100644 --- a/SEFramework/SEFramework/ASDF/AsdfFile.h +++ b/SEFramework/SEFramework/ASDF/AsdfFile.h @@ -42,7 +42,7 @@ namespace SourceXtractor { * @brief represents access to a whole ASDF file * */ -class AsdfFile : public std::enable_shared_from_this { +class AsdfFile { public: /** * Exception thrown when trying to access a value that does not exist in the ASDF tree @@ -54,6 +54,12 @@ class AsdfFile : public std::enable_shared_from_this { */ class AsdfValueTypeMismatchException : public Elements::Exception {}; + /** + * Exception thrown when an ASDF ndarray contains a datatype not (currently) supported by + * SE++ + */ + class AsdfUnsupportedDatatypeException: public Elements::Exception {}; + AsdfFile(const boost::filesystem::path& path, bool writeable); AsdfFile(const boost::filesystem::path& path) : AsdfFile(path, false) {} @@ -62,7 +68,7 @@ class AsdfFile : public std::enable_shared_from_this { virtual ~AsdfFile(); - asdf_file_t* getAsdfFilePtr(); + asdf_file_t* getAsdfFilePtr() const { return m_asdf_ptr.get(); } /** * Wrapper around asdf_ndarray_t @@ -77,11 +83,25 @@ class AsdfFile : public std::enable_shared_from_this { uint32_t ndim() const { return m_ndarray_ptr->ndim; } - std::vector shape() const { - return std::vector(m_ndarray_ptr->shape, - m_ndarray_ptr->shape + m_ndarray_ptr->ndim); + std::vector& shape() const { + if (!m_shape) { + m_shape = std::make_unique>( + m_ndarray_ptr->shape, + m_ndarray_ptr->shape + m_ndarray_ptr->ndim + ); + } + return *m_shape; } + ImageTile::ImageType getImageType() const { return m_image_type; } + + /** + * Convenience method to test whether the ndarray can be read as an image by SE++ + * + * That is, it has ndim == 2 || ndim == 3 (for data cubes) and a supported datatype. + */ + bool isSupportedImage() const; + /** * Given a shared pointer to an already allocated ImageTile (i.e. from ImageTile::create) * copy a tile from the raw ndarray data into the ImageTile. @@ -92,22 +112,23 @@ class AsdfFile : public std::enable_shared_from_this { /** * Private constructor for creating the `Ndarray` wrapper from a raw asdf_value_t * */ - explicit Ndarray(std::shared_ptr file, asdf_value_t *ptr); + explicit Ndarray(const AsdfFile& file, asdf_value_t *ptr); /** * Private constructor for creating the `Ndarray` wrapper from a raw asdf_ndarray_t * */ - explicit Ndarray(std::shared_ptr file, asdf_ndarray_t *ptr) - : m_file(file), m_ndarray_ptr(ptr) {} + explicit Ndarray(asdf_ndarray_t *ptr) + : m_ndarray_ptr(ptr) {} - std::shared_ptr m_file; asdf_ndarray_t* m_ndarray_ptr; + ImageTile::ImageType m_image_type; + mutable std::unique_ptr> m_shape; }; /** * Return the N-th ndarray from the top-level of the ASDF tree iterating the top-level * keys in order. */ - Ndarray getNdarray(int index); + std::unique_ptr getNdarray(int index); /** * Return any ndarray from any path in the ASDF tree @@ -116,7 +137,7 @@ class AsdfFile : public std::enable_shared_from_this { * If the given path exists but is not an ndarray, an AsdfValueTypeMismatchException is * thrown. */ - Ndarray getNdarray(const std::string& path); + std::unique_ptr getNdarray(const std::string& path); /* TODO: More general methods for reading metadata from the ASDF tree; for the first version * not needed though. */ diff --git a/SEFramework/SEFramework/ASDF/AsdfImageSource.h b/SEFramework/SEFramework/ASDF/AsdfImageSource.h new file mode 100644 index 000000000..155dc4334 --- /dev/null +++ b/SEFramework/SEFramework/ASDF/AsdfImageSource.h @@ -0,0 +1,135 @@ +/** + * 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 + */ +/* + * AsdfImageSource.h + * + * Created on: Sep 03, 2025 + * Author: embray + */ + +#ifndef _SEFRAMEWORK_IMAGE_ASDFIMAGESOURCE_H_ +#define _SEFRAMEWORK_IMAGE_ASDFIMAGESOURCE_H_ + +#include +#include + +#include "SEFramework/ASDF/AsdfFile.h" +#include "SEFramework/Image/ImageSource.h" + + +namespace SourceXtractor { + +using Euclid::FilePool::FileManager; +using Euclid::FilePool::FileHandler; + +/** + * Read images from ASDF files + * + * Similar to FitsImageSource but simpler currently as it only allows reading image tiles from + * ASDF, not writing, nor does it support metadata lookup (yet). These features will come later. + */ +class AsdfImageSource : public ImageSource, public std::enable_shared_from_this { +public: + + /** + * Exception raised when instantiating an AsdfImageSource with an ndarray path that is + * not a supported image type (not 2- or 3-D, and not a supported ImageTile::ImageType + */ + class InvalidImageNdarrayException : public Elements::Exception {}; + + /** + * Constructor + * @param filename + * Path to the ASDF file + * @param image_index + * An integer index of the the ndarray to load the image from. When using a numerical + * index, only ndarrays in the top-level of the ASDF tree are considered in the enumeration, + * with ndarrays containing non-image data (i.e. not 2 or 3 dimensions, or containing a + * non-scalar datatype) skipped over. By default the first image ndarray found in the + * file is loaded. + */ + explicit AsdfImageSource(const std::string& filename, int image_index = 0, + std::shared_ptr manager = FileManager::getDefault()) : + AsdfImageSource(filename, image_index, std::nullopt, std::move(manager)) {} + + /** + * Constructor + * @param filename + * Path to the ASDF file + * @param image_path + * JSON Path into the tree of the image ndarray to load. This can be any arbitrarily + * nested ndarray. ndarrays in the root of the tree can be addressed just like "key" + * or "/key" though the leading slash may be omitted. If the given path is not for a + * valid image ndarray an AsdfFile::AsdfValueTypeMismatchException is thrown. + */ + explicit AsdfImageSource(const std::string& filename, const std::string& image_path, + std::shared_ptr manager = FileManager::getDefault()) : + AsdfImageSource(filename, 0, image_path, std::move(manager)) {} + + std::shared_ptr getImageTile(int x, int y, int width, int height) const override; + + std::string getRepr() const override { + return m_filename; + } + + /** + * Returns the width of the image in pixels + */ + int getWidth() const override { + return m_ndarray->shape().at(m_ndarray->ndim() - 2); + } + + /** + * Returns the height of the image in pixels + */ + int getHeight() const override { + return m_ndarray->shape().at(m_ndarray->ndim() - 1); + } + + /** + * Returns the height of the image in pixels + */ + int getDepth() const { + return m_ndarray->ndim() == 3 ? m_ndarray->shape().at(0) : 0; + } + + ImageTile::ImageType getType() const override { + return m_ndarray->getImageType(); + } + + void setLayer(int layer); + + void saveTile(ImageTile& tile) override; + +private: + AsdfImageSource(const std::string& filename, int image_index, + std::optional image_path, + std::shared_ptr manager); + + std::string m_filename; + std::shared_ptr m_file_manager; + std::shared_ptr m_handler; + + std::string m_image_path; + + std::shared_ptr m_ndarray; + int m_current_layer; +}; +} + +#endif /* _SEFRAMEWORK_IMAGE_ASDFIMAGESOURCE_H_ */ diff --git a/SEFramework/SEFramework/ASDF/AsdfReader.h b/SEFramework/SEFramework/ASDF/AsdfReader.h new file mode 100644 index 000000000..9ad25b65a --- /dev/null +++ b/SEFramework/SEFramework/ASDF/AsdfReader.h @@ -0,0 +1,56 @@ +/** 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 SEFramework/ASDF/AsdfReader.h + * @date 09/12/25 + * @author embray + */ + +#ifndef _SEFRAMEWORK_ASDF_ASDFREADER_H +#define _SEFRAMEWORK_ASDF_ASDFREADER_H + +#include "SEFramework/Image/BufferedImage.h" +#include "SEFramework/Image/ImageFileReader.h" +#include "SEFramework/ASDF/AsdfImageSource.h" + +namespace SourceXtractor { + +/** + * @class AsdfReader + * @brief + * + */ +class AsdfReader : public ImageFileReader { + +public: + using ImageFileReader::ImageFileReader; + + static bool test(std::istream& stream); + + using ImageFileReader::get; + /** + * Get the N-th supported image ndarray from the ASDF file + */ + std::shared_ptr get(int image_index) override; + std::shared_ptr get(const std::string& image_path) override; +}; /* End of AsdfReader class */ + +} /* namespace SourceXtractor */ + + +#endif /* _SEFRAMEWORK_ASDF_ASDFREADER_H */ + diff --git a/SEFramework/SEFramework/FITS/FitsReader.h b/SEFramework/SEFramework/FITS/FitsReader.h index 6f0128867..5d0afe775 100644 --- a/SEFramework/SEFramework/FITS/FitsReader.h +++ b/SEFramework/SEFramework/FITS/FitsReader.h @@ -1,4 +1,4 @@ -/** Copyright © 2019 Université de Genève, LMU Munich - Faculty of Physics, IAP-CNRS/Sorbonne Université +/** 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 @@ -15,13 +15,13 @@ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ /** - * @file SEFramework/Image/FitsReader.h + * @file SEFramework/FITS/FitsReader.h * @date 09/01/25 * @author embray */ -#ifndef _SEFRAMEWORK_IMAGE_FITSREADER_H -#define _SEFRAMEWORK_IMAGE_FITSREADER_H +#ifndef _SEFRAMEWORK_FITS_FITSREADER_H +#define _SEFRAMEWORK_FITS_FITSREADER_H #include "SEFramework/Image/BufferedImage.h" #include "SEFramework/Image/ImageFileReader.h" @@ -71,4 +71,5 @@ class FitsReader : public ImageFileReader { } /* namespace SourceXtractor */ -#endif +#endif /* _SEFRAMEWORK_FITS_FITSREADER_H */ + diff --git a/SEFramework/src/lib/ASDF/AsdfFile.cpp b/SEFramework/src/lib/ASDF/AsdfFile.cpp index 47e6cac00..f9aa8d965 100644 --- a/SEFramework/src/lib/ASDF/AsdfFile.cpp +++ b/SEFramework/src/lib/ASDF/AsdfFile.cpp @@ -52,10 +52,6 @@ AsdfFile::AsdfFile(const boost::filesystem::path& path, bool writeable) AsdfFile::~AsdfFile() {} -asdf_file_t* AsdfFile::getAsdfFilePtr() { - return m_asdf_ptr.get(); -} - void AsdfFile::open() { asdf_file_t* ptr = asdf_open_file(m_path.native().c_str(), "r"); @@ -78,7 +74,7 @@ void AsdfFile::open() { } // TODO: Support negative indexing as well -AsdfFile::Ndarray AsdfFile::getNdarray(int index) { +std::unique_ptr AsdfFile::getNdarray(int index) { AsdfValuePtr root = getValue("/"); if (!asdf_value_is_mapping(root.get())) { @@ -94,7 +90,7 @@ AsdfFile::Ndarray AsdfFile::getNdarray(int index) { asdf_value_t* value = asdf_mapping_item_value(item); if (asdf_value_is_ndarray(value)) { if (count == index) { - return Ndarray(shared_from_this(), value); + return std::unique_ptr(new Ndarray(*this, value)); } else if (count > index) { break; } @@ -106,9 +102,9 @@ AsdfFile::Ndarray AsdfFile::getNdarray(int index) { << m_path.native(); } -AsdfFile::Ndarray AsdfFile::getNdarray(const std::string &path) { +std::unique_ptr AsdfFile::getNdarray(const std::string &path) { AsdfValuePtr value = getValue(path); - return Ndarray(shared_from_this(), value.get()); + return std::unique_ptr(new Ndarray(*this, value.get())); } @@ -123,8 +119,40 @@ AsdfFile::AsdfValuePtr AsdfFile::getValue(const std::string& path) { } -AsdfFile::Ndarray::Ndarray(std::shared_ptr file, asdf_value_t *value) - : Ndarray(file, (asdf_ndarray_t*)nullptr) { +static ImageTile::ImageType convertImageType(asdf_datatype_t datatype) { + ImageTile::ImageType image_type; + + switch (datatype) { + case ASDF_DATATYPE_FLOAT32: + image_type = ImageTile::FloatImage; + break; + case ASDF_DATATYPE_FLOAT64: + image_type = ImageTile::DoubleImage; + break; + case ASDF_DATATYPE_INT32: + image_type = ImageTile::IntImage; + break; + case ASDF_DATATYPE_UINT32: + image_type = ImageTile::UIntImage; + break; + case ASDF_DATATYPE_INT64: + image_type = ImageTile::LongLongImage; + break; + default: + // TODO: Support more datatypes supported by ASDF + // Currently SE++ (and ImageTile::ImageType) is constrainted by the basic BITPIX datatypes + // supported by FITS. There's no strong need for that other than the fact that it currently + // only supports FITS. Nevertheless for now this will cover most common cases. + throw AsdfFile::AsdfUnsupportedDatatypeException() << "Unsupported ASDF ndarray datatype: " + << asdf_ndarray_datatype_to_string(datatype); + } + + return image_type; +} + + +AsdfFile::Ndarray::Ndarray(const AsdfFile& file, asdf_value_t *value) + : Ndarray((asdf_ndarray_t*)nullptr) { asdf_ndarray_t *ndarray_ptr = nullptr; asdf_value_err_t err = asdf_value_as_ndarray(value, &ndarray_ptr); switch (err) { @@ -134,15 +162,16 @@ AsdfFile::Ndarray::Ndarray(std::shared_ptr file, asdf_value_t *value) case ASDF_VALUE_ERR_TYPE_MISMATCH: { const char* path = asdf_value_path(value); throw AsdfValueTypeMismatchException() << "Value at " << path << " is not an ndarray: " - << file->m_path.native(); + << file.m_path.native(); } default: { - const char *error_message = asdf_error(file->getAsdfFilePtr()); + const char *error_message = asdf_error(file.getAsdfFilePtr()); throw AsdfValueNotFoundException() << "An error occurred reading the ASDF file " - << file->m_path.native() << ": " << error_message; + << file.m_path.native() << ": " << error_message; } } m_ndarray_ptr = ndarray_ptr; + m_image_type = convertImageType(ndarray_ptr->datatype); } @@ -175,4 +204,21 @@ void AsdfFile::Ndarray::fillImageTile(const std::shared_ptr image_til } } + +bool AsdfFile::Ndarray::isSupportedImage() const { + uint64_t ndim = m_ndarray_ptr->ndim; + + if (ndim < 2 || ndim > 3) { + return false; + } + + try { + getImageType(); + } catch (const AsdfFile::AsdfUnsupportedDatatypeException&) { + return false; + } + + return true; +} + } // namespace SourceXtractor diff --git a/SEFramework/src/lib/ASDF/AsdfImageSource.cpp b/SEFramework/src/lib/ASDF/AsdfImageSource.cpp new file mode 100644 index 000000000..24ded31bc --- /dev/null +++ b/SEFramework/src/lib/ASDF/AsdfImageSource.cpp @@ -0,0 +1,118 @@ +/** + * 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 + */ +/* + * AsdfImageSource.cpp + * + * Created on: Sep 03, 2025 + * Author: embray + */ + +#include + +#include "SEFramework/ASDF/AsdfFile.h" +#include "SEFramework/ASDF/AsdfImageSource.h" + + +namespace SourceXtractor { +AsdfImageSource::AsdfImageSource(const std::string& filename, int image_index, + std::optional image_path, + std::shared_ptr manager) + : m_filename(filename) + , m_file_manager(std::move(manager)) + , m_handler(m_file_manager->getFileHandler(filename)) { + + auto acc = m_handler->getAccessor(); + auto& file = acc->m_fd; + + if (image_path) { + m_ndarray = std::move(file.getNdarray(image_path.value())); + } else { + // Find the 'th ndarray that actually looks like a usable image + int found = 0; + int ndarray_index = 0; + while (found <= image_index) { + try { + auto ndarray = file.getNdarray(ndarray_index); + + if (ndarray->isSupportedImage()) { + if (found == image_index) { + m_ndarray = std::move(ndarray); + break; + } else { + found++; + } + } + } catch (const AsdfFile::AsdfUnsupportedDatatypeException&) { + // continue + } catch (const AsdfFile::AsdfValueNotFoundException&) { + break; + } + ndarray_index++; + } + + if (!m_ndarray) { + throw AsdfImageSource::InvalidImageNdarrayException(); + } + } + + uint64_t ndim = m_ndarray->ndim(); + if (ndim < 2 || ndim > 3) { + if (image_path) { + throw InvalidImageNdarrayException() << "Can't find 2D image or data cube at " + << image_path.value() << "in ASDF file: " << filename; + } else { + throw InvalidImageNdarrayException() << "Can't find 2D image or data cube in ASDF file: " + << filename; + } + } + + try { + m_ndarray->getImageType(); + } catch (const AsdfFile::AsdfUnsupportedDatatypeException&) { + if (image_path) { + throw InvalidImageNdarrayException() << "Unsupported datatype for ndarray at " + << image_path.value() << "in ASDF file: " << filename; + } else { + throw InvalidImageNdarrayException() << "Unsupported datatype in ASDF file: " + << filename; + } + } +} + + +void AsdfImageSource::setLayer(int layer) { + uint64_t depth = m_ndarray->ndim() >= 3 ? m_ndarray->shape().at(0) : 0; + if (layer < 0 || layer >= static_cast(depth)) { + throw Elements::Exception() << "Trying to access an inexistent data cube layer (" << layer << ") in " << m_filename; + } + m_current_layer = layer; +} + + +std::shared_ptr AsdfImageSource::getImageTile(int x, int y, int width, int height) const { + auto tile = ImageTile::create(m_ndarray->getImageType(), x, y, width, height, + std::const_pointer_cast(shared_from_this())); + m_ndarray->fillImageTile(tile, m_current_layer); + return tile; +} + + +void AsdfImageSource::saveTile(ImageTile& /* tile */) { + throw Elements::Exception() << "Saving to ASDF not yet supported"; +}; +} diff --git a/SEFramework/src/lib/ASDF/AsdfReader.cpp b/SEFramework/src/lib/ASDF/AsdfReader.cpp new file mode 100644 index 000000000..9667304ed --- /dev/null +++ b/SEFramework/src/lib/ASDF/AsdfReader.cpp @@ -0,0 +1,65 @@ +/** 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 SEFramework/ASDF/AsdfReader.cpp + * @date 09/12/25 + * @author embray + */ + +#include + +#include + +#include "SEFramework/ASDF/AsdfFile.h" +#include "SEFramework/ASDF/AsdfReader.h" +#include "SEFramework/Image/ImageFileReader.h" +#include "SEFramework/Image/ImageSource.h" + + +namespace SourceXtractor { + + +bool AsdfReader::test(std::istream& stream) { + // Should be more than enough to determine the #ASDF header comment + std::string buf(256, '\0'); + if (!stream.read(&buf[0], 256)) { + return false; + } + + buf.erase(buf.find_first_of("\r\n")); + boost::regex asdf_header_regex(R"(^#ASDF \d+\.\d+\.\d+$)"); + return boost::regex_match(buf, asdf_header_regex); +} + + +std::shared_ptr AsdfReader::get(int image_index) { + return std::make_shared(m_filename, image_index); +} + + +std::shared_ptr AsdfReader::get(const std::string& extname) { + return std::make_shared(m_filename, extname); +} + + +static ImageFileReader::StaticFileType asdf_reader{ + "asdf", + {"asdf"} +}; + +} + diff --git a/SEFramework/src/lib/FITS/FitsReader.cpp b/SEFramework/src/lib/FITS/FitsReader.cpp index 3486fec0c..e50673916 100644 --- a/SEFramework/src/lib/FITS/FitsReader.cpp +++ b/SEFramework/src/lib/FITS/FitsReader.cpp @@ -15,7 +15,7 @@ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ /** - * @file SEFramework/Image/FitsReader.cpp + * @file SEFramework/FITS/FitsReader.cpp * @date 09/01/25 * @author embray */ diff --git a/SEFramework/tests/src/ASDF/AsdfFile_test.cpp b/SEFramework/tests/src/ASDF/AsdfFile_test.cpp index 99fa83a12..6fa200b2d 100644 --- a/SEFramework/tests/src/ASDF/AsdfFile_test.cpp +++ b/SEFramework/tests/src/ASDF/AsdfFile_test.cpp @@ -53,10 +53,10 @@ BOOST_AUTO_TEST_CASE( missing_file ) { //----------------------------------------------------------------------------- BOOST_AUTO_TEST_CASE( get_ndarray_by_index ) { - auto asdf_file = std::make_shared(Elements::getAuxiliaryPath("basic.asdf").native()); - auto ndarray = asdf_file->getNdarray(0); - BOOST_CHECK_EQUAL(ndarray.ndim(), 1); - auto shape = ndarray.shape(); + AsdfFile asdf_file(Elements::getAuxiliaryPath("basic.asdf").native()); + auto ndarray = asdf_file.getNdarray(0); + BOOST_CHECK_EQUAL(ndarray->ndim(), 1); + auto shape = ndarray->shape(); std::vector expected_shape{8}; BOOST_CHECK_EQUAL_COLLECTIONS( shape.begin(), shape.end(), @@ -67,10 +67,10 @@ BOOST_AUTO_TEST_CASE( get_ndarray_by_index ) { //----------------------------------------------------------------------------- BOOST_AUTO_TEST_CASE( get_ndarray_by_path ) { - auto asdf_file = std::make_shared(Elements::getAuxiliaryPath("basic.asdf").native()); - auto ndarray = asdf_file->getNdarray("data"); - BOOST_CHECK_EQUAL(ndarray.ndim(), 1); - auto shape = ndarray.shape(); + AsdfFile asdf_file(Elements::getAuxiliaryPath("basic.asdf").native()); + auto ndarray = asdf_file.getNdarray("data"); + BOOST_CHECK_EQUAL(ndarray->ndim(), 1); + auto shape = ndarray->shape(); std::vector expected_shape{8}; BOOST_CHECK_EQUAL_COLLECTIONS( shape.begin(), shape.end(), From 69420f4057f4b3d6c9ca9bd1534969e4f1673c32 Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Mon, 15 Sep 2025 14:32:00 +0200 Subject: [PATCH 11/29] support specifying the desired ImageType for the ImageSources returned by the ImageFileReader.get() and iterator interfaces Make this work again for FITS--for ASDF we need more work on the libasdf side. --- .../SEFramework/ASDF/AsdfImageSource.h | 9 ++++- SEFramework/SEFramework/ASDF/AsdfReader.h | 7 +++- SEFramework/SEFramework/FITS/FitsReader.h | 10 +++-- .../SEFramework/Image/ImageFileReader.h | 38 +++++++++++++++---- SEFramework/src/lib/ASDF/AsdfImageSource.cpp | 3 ++ SEFramework/src/lib/ASDF/AsdfReader.cpp | 9 +++-- SEFramework/src/lib/FITS/FitsReader.cpp | 13 ++++--- SEFramework/src/lib/Image/ImageFileReader.cpp | 4 +- .../Configuration/DetectionImageConfig.cpp | 2 +- 9 files changed, 67 insertions(+), 28 deletions(-) diff --git a/SEFramework/SEFramework/ASDF/AsdfImageSource.h b/SEFramework/SEFramework/ASDF/AsdfImageSource.h index 155dc4334..0a605fce0 100644 --- a/SEFramework/SEFramework/ASDF/AsdfImageSource.h +++ b/SEFramework/SEFramework/ASDF/AsdfImageSource.h @@ -30,6 +30,7 @@ #include "SEFramework/ASDF/AsdfFile.h" #include "SEFramework/Image/ImageSource.h" +#include "SEFramework/Image/ImageTile.h" namespace SourceXtractor { @@ -64,8 +65,9 @@ class AsdfImageSource : public ImageSource, public std::enable_shared_from_this< * file is loaded. */ explicit AsdfImageSource(const std::string& filename, int image_index = 0, + ImageTile::ImageType image_type = ImageTile::AutoType, std::shared_ptr manager = FileManager::getDefault()) : - AsdfImageSource(filename, image_index, std::nullopt, std::move(manager)) {} + AsdfImageSource(filename, image_index, std::nullopt, image_type, std::move(manager)) {} /** * Constructor @@ -78,8 +80,9 @@ class AsdfImageSource : public ImageSource, public std::enable_shared_from_this< * valid image ndarray an AsdfFile::AsdfValueTypeMismatchException is thrown. */ explicit AsdfImageSource(const std::string& filename, const std::string& image_path, + ImageTile::ImageType image_type = ImageTile::AutoType, std::shared_ptr manager = FileManager::getDefault()) : - AsdfImageSource(filename, 0, image_path, std::move(manager)) {} + AsdfImageSource(filename, 0, image_path, image_type, std::move(manager)) {} std::shared_ptr getImageTile(int x, int y, int width, int height) const override; @@ -119,9 +122,11 @@ class AsdfImageSource : public ImageSource, public std::enable_shared_from_this< private: AsdfImageSource(const std::string& filename, int image_index, std::optional image_path, + ImageTile::ImageType image_type, std::shared_ptr manager); std::string m_filename; + ImageTile::ImageType m_image_type; std::shared_ptr m_file_manager; std::shared_ptr m_handler; diff --git a/SEFramework/SEFramework/ASDF/AsdfReader.h b/SEFramework/SEFramework/ASDF/AsdfReader.h index 9ad25b65a..4671566dd 100644 --- a/SEFramework/SEFramework/ASDF/AsdfReader.h +++ b/SEFramework/SEFramework/ASDF/AsdfReader.h @@ -25,6 +25,7 @@ #include "SEFramework/Image/BufferedImage.h" #include "SEFramework/Image/ImageFileReader.h" +#include "SEFramework/Image/ImageTile.h" #include "SEFramework/ASDF/AsdfImageSource.h" namespace SourceXtractor { @@ -45,8 +46,10 @@ class AsdfReader : public ImageFileReader { /** * Get the N-th supported image ndarray from the ASDF file */ - std::shared_ptr get(int image_index) override; - std::shared_ptr get(const std::string& image_path) override; + std::shared_ptr get(int image_index, + ImageTile::ImageType image_type = ImageTile::AutoType) override; + std::shared_ptr get(const std::string& image_path, + ImageTile::ImageType image_type = ImageTile::AutoType) override; }; /* End of AsdfReader class */ } /* namespace SourceXtractor */ diff --git a/SEFramework/SEFramework/FITS/FitsReader.h b/SEFramework/SEFramework/FITS/FitsReader.h index 5d0afe775..aa4f130e8 100644 --- a/SEFramework/SEFramework/FITS/FitsReader.h +++ b/SEFramework/SEFramework/FITS/FitsReader.h @@ -25,6 +25,7 @@ #include "SEFramework/Image/BufferedImage.h" #include "SEFramework/Image/ImageFileReader.h" +#include "SEFramework/Image/ImageTile.h" #include "SEFramework/FITS/FitsImageSource.h" namespace SourceXtractor { @@ -51,10 +52,13 @@ class FitsReader : public ImageFileReader { * * Use FitsReader.getHdu to get by FITS HDU number */ - std::shared_ptr get(int image_index) override; - std::shared_ptr get(const std::string& image_path) override; + std::shared_ptr get(int image_index, + ImageTile::ImageType image_type = ImageTile::AutoType) override; + std::shared_ptr get(const std::string& image_path, + ImageTile::ImageType image_type = ImageTile::AutoType) override; - std::shared_ptr getHdu(int hdu_num); + std::shared_ptr getHdu(int hdu_num, + ImageTile::ImageType image_type = ImageTile::AutoType); // TODO: Integrate this interface into the ImageFileReader base class template diff --git a/SEFramework/SEFramework/Image/ImageFileReader.h b/SEFramework/SEFramework/Image/ImageFileReader.h index abce18f95..c90d77d9e 100644 --- a/SEFramework/SEFramework/Image/ImageFileReader.h +++ b/SEFramework/SEFramework/Image/ImageFileReader.h @@ -26,6 +26,7 @@ #include #include "SEFramework/Image/ImageSource.h" +#include "SEFramework/Image/ImageTile.h" namespace SourceXtractor { @@ -62,9 +63,15 @@ class ImageFileReader { virtual ~ImageFileReader() = default; - virtual std::shared_ptr get(); - virtual std::shared_ptr get(int image_index) = 0; - virtual std::shared_ptr get(const std::string& image_path) = 0; + /** + * Equivalent to calling get(m_image_index/m_image_path, ImageTile::AutoType) + */ + std::shared_ptr get(); + + virtual std::shared_ptr get( + int image_index, ImageTile::ImageType image_type = ImageTile::AutoType) = 0; + virtual std::shared_ptr get( + const std::string& image_path, ImageTile::ImageType image_type = ImageTile::AutoType) = 0; /* ImageFileReader iterator interface */ class Iterator { @@ -75,10 +82,11 @@ class ImageFileReader { using pointer = std::shared_ptr*; using reference = std::shared_ptr&; - Iterator() noexcept : m_reader(nullptr), m_index(0) {} + Iterator() noexcept : m_reader(nullptr), m_index(0), m_image_type(ImageTile::AutoType) {} - Iterator(ImageFileReader* reader, int index) - : m_reader(reader), m_index(index) { + Iterator(ImageFileReader* reader, int index, + ImageTile::ImageType image_type = ImageTile::AutoType) + : m_reader(reader), m_index(index), m_image_type(image_type) { advance(); } @@ -98,6 +106,15 @@ class ImageFileReader { return !(*this == other); } + /* To support ImageFileReader->iter(ImageTile::ImageType) */ + Iterator begin() { + return *this; + } + + Iterator end() { + return Iterator(); + } + private: void advance() { // If m_reader is null the iterator is considered "done" @@ -107,7 +124,7 @@ class ImageFileReader { } try { - m_current = m_reader->get(m_index); + m_current = m_reader->get(m_index, m_image_type); } catch (...) { m_current = nullptr; } @@ -116,10 +133,15 @@ class ImageFileReader { ImageFileReader* m_reader; int m_index; std::shared_ptr m_current; + ImageTile::ImageType m_image_type; }; + Iterator iter(ImageTile::ImageType image_type = ImageTile::AutoType) { + return Iterator(this, 0, image_type); + } + Iterator begin() { - return Iterator(this, 0); + return Iterator(this, 0, ImageTile::AutoType); } Iterator end() { diff --git a/SEFramework/src/lib/ASDF/AsdfImageSource.cpp b/SEFramework/src/lib/ASDF/AsdfImageSource.cpp index 24ded31bc..49818017b 100644 --- a/SEFramework/src/lib/ASDF/AsdfImageSource.cpp +++ b/SEFramework/src/lib/ASDF/AsdfImageSource.cpp @@ -31,8 +31,10 @@ namespace SourceXtractor { AsdfImageSource::AsdfImageSource(const std::string& filename, int image_index, std::optional image_path, + ImageTile::ImageType image_type, std::shared_ptr manager) : m_filename(filename) + , m_image_type(image_type) , m_file_manager(std::move(manager)) , m_handler(m_file_manager->getFileHandler(filename)) { @@ -105,6 +107,7 @@ void AsdfImageSource::setLayer(int layer) { std::shared_ptr AsdfImageSource::getImageTile(int x, int y, int width, int height) const { + // TODO: (#5) support image data type conversion auto tile = ImageTile::create(m_ndarray->getImageType(), x, y, width, height, std::const_pointer_cast(shared_from_this())); m_ndarray->fillImageTile(tile, m_current_layer); diff --git a/SEFramework/src/lib/ASDF/AsdfReader.cpp b/SEFramework/src/lib/ASDF/AsdfReader.cpp index 9667304ed..3f55369c0 100644 --- a/SEFramework/src/lib/ASDF/AsdfReader.cpp +++ b/SEFramework/src/lib/ASDF/AsdfReader.cpp @@ -46,13 +46,14 @@ bool AsdfReader::test(std::istream& stream) { } -std::shared_ptr AsdfReader::get(int image_index) { - return std::make_shared(m_filename, image_index); +std::shared_ptr AsdfReader::get(int image_index, ImageTile::ImageType image_type) { + return std::make_shared(m_filename, image_index, image_type); } -std::shared_ptr AsdfReader::get(const std::string& extname) { - return std::make_shared(m_filename, extname); +std::shared_ptr AsdfReader::get(const std::string& extname, + ImageTile::ImageType image_type) { + return std::make_shared(m_filename, extname, image_type); } diff --git a/SEFramework/src/lib/FITS/FitsReader.cpp b/SEFramework/src/lib/FITS/FitsReader.cpp index e50673916..2845c6961 100644 --- a/SEFramework/src/lib/FITS/FitsReader.cpp +++ b/SEFramework/src/lib/FITS/FitsReader.cpp @@ -40,7 +40,7 @@ bool FitsReader::test(std::istream& stream) { } -std::shared_ptr FitsReader::get(int image_index) { +std::shared_ptr FitsReader::get(int image_index, ImageTile::ImageType image_type) { // TODO: It turns out none of this is necessary: It's already implemented in // the FitsFile class, which actually loops over all the HDUs when opening the file // and gets the HDU numbers of images (can be accessed with FitsFile::getImageHdus) @@ -48,7 +48,7 @@ std::shared_ptr FitsReader::get(int image_index) { auto it = m_image_hdu_map.find(image_index); if (it != m_image_hdu_map.end()) { - return std::make_shared(m_filename, it->second); + return std::make_shared(m_filename, it->second, image_type); } // Try loading HDUs from the file until we find the next image HDU @@ -81,13 +81,14 @@ std::shared_ptr FitsReader::get(int image_index) { } -std::shared_ptr FitsReader::get(const std::string& extname) { - return std::make_shared(m_filename, extname); +std::shared_ptr FitsReader::get(const std::string& extname, + ImageTile::ImageType image_type) { + return std::make_shared(m_filename, extname, image_type); } -std::shared_ptr FitsReader::getHdu(int hdu_num) { - return std::make_shared(m_filename, hdu_num); +std::shared_ptr FitsReader::getHdu(int hdu_num, ImageTile::ImageType image_type) { + return std::make_shared(m_filename, hdu_num, image_type); } diff --git a/SEFramework/src/lib/Image/ImageFileReader.cpp b/SEFramework/src/lib/Image/ImageFileReader.cpp index 338bd1037..c1c3f905e 100644 --- a/SEFramework/src/lib/Image/ImageFileReader.cpp +++ b/SEFramework/src/lib/Image/ImageFileReader.cpp @@ -92,9 +92,9 @@ ImageFileReader::ImageFileReader(const std::string& filename, int image_index) std::shared_ptr ImageFileReader::get() { if (m_image_index >= 0) { - return get(m_image_index); + return get(m_image_index, ImageTile::AutoType); } else if (m_image_path) { - return get(*m_image_path); + return get(*m_image_path, ImageTile::AutoType); } return get(0); diff --git a/SEImplementation/src/lib/Configuration/DetectionImageConfig.cpp b/SEImplementation/src/lib/Configuration/DetectionImageConfig.cpp index 8f0d6ed9f..9fcafabc4 100644 --- a/SEImplementation/src/lib/Configuration/DetectionImageConfig.cpp +++ b/SEImplementation/src/lib/Configuration/DetectionImageConfig.cpp @@ -98,7 +98,7 @@ void DetectionImageConfig::initialize(const UserValues& args) { auto image_reader = ImageFileReader::create(m_detection_image_path); - for (const auto& img_source: *image_reader) { + for (const auto& img_source: image_reader->iter(ImageTile::FloatImage)) { DetectionImageExtension extension = DetectionImageExtension::create(img_source, args); m_extensions.emplace_back(std::move(extension)); } From 559048723d9c22f9b08191bd948280b92577ae81 Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Fri, 3 Oct 2025 15:53:32 +0200 Subject: [PATCH 12/29] Support image type conversion in AsdfImageSource, building on top of latest updates to libasdf's development API --- SEFramework/src/lib/ASDF/AsdfFile.cpp | 36 ++++++++++++++++++++++++--- 1 file changed, 32 insertions(+), 4 deletions(-) diff --git a/SEFramework/src/lib/ASDF/AsdfFile.cpp b/SEFramework/src/lib/ASDF/AsdfFile.cpp index f9aa8d965..c001534b4 100644 --- a/SEFramework/src/lib/ASDF/AsdfFile.cpp +++ b/SEFramework/src/lib/ASDF/AsdfFile.cpp @@ -119,10 +119,10 @@ AsdfFile::AsdfValuePtr AsdfFile::getValue(const std::string& path) { } -static ImageTile::ImageType convertImageType(asdf_datatype_t datatype) { +static ImageTile::ImageType convertDatatypeToImageType(const asdf_datatype_t *datatype) { ImageTile::ImageType image_type; - switch (datatype) { + switch (datatype->type) { case ASDF_DATATYPE_FLOAT32: image_type = ImageTile::FloatImage; break; @@ -144,13 +144,40 @@ static ImageTile::ImageType convertImageType(asdf_datatype_t datatype) { // supported by FITS. There's no strong need for that other than the fact that it currently // only supports FITS. Nevertheless for now this will cover most common cases. throw AsdfFile::AsdfUnsupportedDatatypeException() << "Unsupported ASDF ndarray datatype: " - << asdf_ndarray_datatype_to_string(datatype); + << asdf_ndarray_datatype_to_string(datatype->type); } return image_type; } +static asdf_scalar_datatype_t convertImageTypeToDatatype(ImageTile::ImageType image_type) { + // This is the same default used for FITS + asdf_scalar_datatype_t datatype = ASDF_DATATYPE_FLOAT32; + + switch (image_type) { + default: + case ImageTile::FloatImage: + datatype = ASDF_DATATYPE_FLOAT32; + break; + case ImageTile::DoubleImage: + datatype = ASDF_DATATYPE_FLOAT64; + break; + case ImageTile::IntImage: + datatype = ASDF_DATATYPE_INT32; + break; + case ImageTile::UIntImage: + datatype = ASDF_DATATYPE_UINT32; + break; + case ImageTile::LongLongImage: + datatype = ASDF_DATATYPE_INT64; + break; + } + + return datatype; +} + + AsdfFile::Ndarray::Ndarray(const AsdfFile& file, asdf_value_t *value) : Ndarray((asdf_ndarray_t*)nullptr) { asdf_ndarray_t *ndarray_ptr = nullptr; @@ -171,7 +198,7 @@ AsdfFile::Ndarray::Ndarray(const AsdfFile& file, asdf_value_t *value) } } m_ndarray_ptr = ndarray_ptr; - m_image_type = convertImageType(ndarray_ptr->datatype); + m_image_type = convertDatatypeToImageType(&ndarray_ptr->datatype); } @@ -185,6 +212,7 @@ void AsdfFile::Ndarray::fillImageTile(const std::shared_ptr image_til image_tile->getWidth(), image_tile->getHeight(), &plane_origin, + convertImageTypeToDatatype(image_tile->getType()), &data ); From d0fe03f227be037eda59b25d69401c56fa807a53 Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Fri, 3 Oct 2025 16:03:43 +0200 Subject: [PATCH 13/29] Support reading ExternalFlagConfig from ASDF files (I'm not even exactly sure what this is for yet, but it nicely demonstrates how relatively easy it is now to abstract out the file type) --- .../ExternalFlag/ExternalFlagConfig.cpp | 55 ++++--------------- 1 file changed, 12 insertions(+), 43 deletions(-) diff --git a/SEImplementation/src/lib/Plugin/ExternalFlag/ExternalFlagConfig.cpp b/SEImplementation/src/lib/Plugin/ExternalFlag/ExternalFlagConfig.cpp index a0202efc3..9e3155b1b 100644 --- a/SEImplementation/src/lib/Plugin/ExternalFlag/ExternalFlagConfig.cpp +++ b/SEImplementation/src/lib/Plugin/ExternalFlag/ExternalFlagConfig.cpp @@ -23,15 +23,10 @@ #include #include -#include -using boost::regex; -using boost::regex_match; -using boost::smatch; - #include "Configuration/ProgramOptionsHelper.h" -#include "SEFramework/FITS/FitsReader.h" - +#include "SEFramework/Image/BufferedImage.h" +#include "SEFramework/Image/ImageFileReader.h" #include "SEImplementation/Plugin/ExternalFlag/ExternalFlagConfig.h" namespace po = boost::program_options; @@ -41,7 +36,7 @@ using poh = Euclid::Configuration::ProgramOptionsHelper; namespace SourceXtractor { namespace { - + const std::string FLAG_IMAGE {"flag-image"}; const std::string FLAG_TYPE {"flag-type"}; @@ -67,7 +62,7 @@ auto ExternalFlagConfig::getProgramOptions() -> std::map()); } - + // Check that the type is a valid option if (available_types.count(type) == 0) { throw Elements::Exception() << "Invalid option " << poh::wildcard(FLAG_TYPE, name) @@ -90,38 +85,15 @@ void ExternalFlagConfig::preInitialize(const UserValues& args) { void ExternalFlagConfig::initialize(const UserValues& args) { for (auto& name : poh::findWildcardNames({FLAG_IMAGE, FLAG_TYPE}, args)) { - + auto& filename = args.at(poh::wildcard(FLAG_IMAGE, name)).as(); std::vector> flag_images; - boost::regex hdu_regex(".*\\[[0-9]*\\]$"); - - for (int i=0;; i++) { - std::shared_ptr fits_image_source; - if (boost::regex_match(filename, hdu_regex)) { - if (i==0) { - fits_image_source = std::make_shared(filename, 0, ImageTile::LongLongImage); - } else { - break; - } - } else { - try { - fits_image_source = std::make_shared(filename, i+1, ImageTile::LongLongImage); - } catch (...) { - if (i==0) { - // Skip past primary HDU if it doesn't have an image - continue; - } else { - if (flag_images.size() == 0) { - throw; - } - break; - } - } - } - - flag_images.emplace_back(BufferedImage::create(fits_image_source)); + auto image_reader = ImageFileReader::create(filename); + + for (const auto& img_source: image_reader->iter(ImageTile::LongLongImage)) { + flag_images.emplace_back(BufferedImage::create(img_source)); } - + std::string type_str; if (args.count(poh::wildcard(FLAG_TYPE, name)) == 0) { type_str = "OR"; @@ -129,7 +101,7 @@ void ExternalFlagConfig::initialize(const UserValues& args) { type_str = boost::to_upper_copy(args.at(poh::wildcard(FLAG_TYPE, name)).as()); } Type type = available_types.at(type_str); - + m_flag_info_list.emplace_back(name, FlagInfo{std::move(flag_images), type}); } } @@ -139,6 +111,3 @@ auto ExternalFlagConfig::getFlagInfoList() const -> const std::vector Date: Fri, 3 Oct 2025 16:28:52 +0200 Subject: [PATCH 14/29] Basic support for reading ASDF images in MeasurementImageConfig Needed to add setLayer in the ImageSource base class for convenience; by default this can just be a no-op where it's not applicable. --- .../SEFramework/ASDF/AsdfImageSource.h | 2 +- .../SEFramework/FITS/FitsImageSource.h | 2 +- SEFramework/SEFramework/Image/ImageSource.h | 3 ++ .../Configuration/MeasurementImageConfig.cpp | 32 ++++++++++++------- 4 files changed, 26 insertions(+), 13 deletions(-) diff --git a/SEFramework/SEFramework/ASDF/AsdfImageSource.h b/SEFramework/SEFramework/ASDF/AsdfImageSource.h index 0a605fce0..4c82f5b48 100644 --- a/SEFramework/SEFramework/ASDF/AsdfImageSource.h +++ b/SEFramework/SEFramework/ASDF/AsdfImageSource.h @@ -115,7 +115,7 @@ class AsdfImageSource : public ImageSource, public std::enable_shared_from_this< return m_ndarray->getImageType(); } - void setLayer(int layer); + void setLayer(int layer) override; void saveTile(ImageTile& tile) override; diff --git a/SEFramework/SEFramework/FITS/FitsImageSource.h b/SEFramework/SEFramework/FITS/FitsImageSource.h index 042d567e6..da186d359 100644 --- a/SEFramework/SEFramework/FITS/FitsImageSource.h +++ b/SEFramework/SEFramework/FITS/FitsImageSource.h @@ -107,7 +107,7 @@ class FitsImageSource : public ImageSource, public std::enable_shared_from_this< return m_depth; } - void setLayer(int layer); + void setLayer(int layer) override; std::shared_ptr getImageTile(int x, int y, int width, int height) const override; diff --git a/SEFramework/SEFramework/Image/ImageSource.h b/SEFramework/SEFramework/Image/ImageSource.h index 227191ece..b74929a1c 100644 --- a/SEFramework/SEFramework/Image/ImageSource.h +++ b/SEFramework/SEFramework/Image/ImageSource.h @@ -72,6 +72,9 @@ class ImageSource { /// Returns the height of the image in pixels virtual int getHeight() const = 0; + /// Sets the current layer to use if the image is in a data cube + virtual void setLayer(int /* unused */) { } + virtual ImageTile::ImageType getType() const = 0; /** diff --git a/SEImplementation/src/lib/Configuration/MeasurementImageConfig.cpp b/SEImplementation/src/lib/Configuration/MeasurementImageConfig.cpp index ba8140805..ad661ba8c 100644 --- a/SEImplementation/src/lib/Configuration/MeasurementImageConfig.cpp +++ b/SEImplementation/src/lib/Configuration/MeasurementImageConfig.cpp @@ -31,9 +31,11 @@ #include +#include #include +#include +#include #include -#include #include #include @@ -81,8 +83,8 @@ void validateImagePaths(const PyMeasurementImage& image) { } std::shared_ptr createMeasurementImage( - std::shared_ptr fits_image_source, double flux_scale) { - std::shared_ptr image = BufferedImage::create(fits_image_source); + std::shared_ptr image_source, double flux_scale) { + std::shared_ptr image = BufferedImage::create(image_source); if (flux_scale != 1.) { image = MultiplyImage::create(image, flux_scale); } @@ -117,8 +119,9 @@ std::shared_ptr createWeightMap(const PyMeasurementImage& py_image) throw Elements::Exception() << "Please give an appropriate weight type for image: " << py_image.weight_file; } - auto weight_image_source = - std::make_shared(py_image.weight_file, py_image.weight_hdu+1, ImageTile::FloatImage); + auto reader = ImageFileReader::create(py_image.weight_file); + // TODO: Extend py_image to also support ASDF ndarray paths + auto weight_image_source = reader->get(py_image.weight_hdu, ImageTile::FloatImage); std::shared_ptr weight_map = BufferedImage::create(weight_image_source); if (py_image.is_data_cube) { weight_image_source->setLayer(py_image.weight_layer); @@ -151,7 +154,7 @@ WeightImage::PixelType extractWeightThreshold(const PyMeasurementImage& py_image } else { threshold = std::numeric_limits::max(); } - break; + break; } return threshold; } @@ -185,15 +188,22 @@ void MeasurementImageConfig::initialize(const UserValues&) { info.m_image_layer = py_image.image_layer; info.m_weight_layer = py_image.weight_layer; - auto fits_image_source = - std::make_shared(py_image.file, py_image.image_hdu+1, ImageTile::FloatImage); + auto reader = ImageFileReader::create(py_image.file); + auto image_source = reader->get(py_image.image_hdu, ImageTile::FloatImage); if (py_image.is_data_cube) { - fits_image_source->setLayer(py_image.image_layer); + image_source->setLayer(py_image.image_layer); } - info.m_measurement_image = createMeasurementImage(fits_image_source, py_image.flux_scale); - info.m_coordinate_system = std::make_shared(*fits_image_source); + info.m_measurement_image = createMeasurementImage(image_source, py_image.flux_scale); + + // TODO: Instantiating a WCS from ASDF files is not yet supported, use the dummy + // WCS for now. + if (auto fits_image_source = std::dynamic_pointer_cast(image_source)) { + info.m_coordinate_system = std::make_shared(*fits_image_source); + } else { + info.m_coordinate_system = std::make_shared(WCS::identity(2)); + } info.m_gain = py_image.gain / flux_scale; info.m_saturation_level = py_image.saturation * flux_scale; From aabcf19943f41ec3bc1175de75043fb793dc73ca Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Fri, 3 Oct 2025 16:34:01 +0200 Subject: [PATCH 15/29] Basic support for reading ASDF files in WeightImageConfig --- .../lib/Configuration/WeightImageConfig.cpp | 42 ++++--------------- 1 file changed, 7 insertions(+), 35 deletions(-) diff --git a/SEImplementation/src/lib/Configuration/WeightImageConfig.cpp b/SEImplementation/src/lib/Configuration/WeightImageConfig.cpp index 593a4f21d..0885e801d 100644 --- a/SEImplementation/src/lib/Configuration/WeightImageConfig.cpp +++ b/SEImplementation/src/lib/Configuration/WeightImageConfig.cpp @@ -23,18 +23,14 @@ #include #include -#include -using boost::regex; -using boost::regex_match; -using boost::smatch; #include "Configuration/ConfigManager.h" +#include "SEFramework/Image/BufferedImage.h" +#include "SEFramework/Image/ImageFileReader.h" #include "SEFramework/Image/ImageSource.h" #include "SEFramework/Image/ProcessingImageSource.h" #include "SEFramework/Image/ProcessedImage.h" -#include "SEFramework/FITS/FitsReader.h" -#include "SEFramework/FITS/FitsImageSource.h" #include "SEImplementation/Configuration/DetectionImageConfig.h" @@ -137,38 +133,14 @@ void WeightImageConfig::initialize(const UserValues& args) { throw Elements::Exception() << "Setting absolute weight but providing *no* weight image does not make sense."; if (weight_image_filename != "") { - boost::regex hdu_regex(".*\\[[0-9]*\\]$"); - - for (int i=0;; i++) { - std::shared_ptr fits_image_source; - if (boost::regex_match(weight_image_filename, hdu_regex)) { - if (i==0) { - fits_image_source = std::make_shared(weight_image_filename, 0, ImageTile::FloatImage); - } else { - break; - } - } else { - try { - fits_image_source = std::make_shared(weight_image_filename, i+1, ImageTile::FloatImage); - } catch (...) { - if (i==0) { - // Skip past primary HDU if it doesn't have an image - continue; - } else { - if (m_weight_images.size() == 0) { - throw; - } - break; - } - } - } - - std::shared_ptr weight_image = BufferedImage::create(fits_image_source); + auto reader = ImageFileReader::create(weight_image_filename); + for (const auto& img_source: reader->iter(ImageTile::FloatImage)) { + std::shared_ptr weight_image = BufferedImage::create( + img_source); weight_image = convertWeightMap(weight_image, m_weight_type, m_weight_scaling); - - // we should have a corresponding detection image auto flux_scale = getDependency().getOriginalFluxScale(m_weight_images.size()); + // we should have a corresponding detection image WeightImage::PixelType scaled_weight_threshold = m_weight_threshold; if (flux_scale != 1. && m_absolute_weight) { weight_image = MultiplyImage::create(weight_image, flux_scale * flux_scale); From d500ccce0ee1baeed44a395b73b347f8463a7967 Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Fri, 3 Oct 2025 16:46:41 +0200 Subject: [PATCH 16/29] Replace the ASDF test files with ones that contain the same data as the FITS tests--just their mirror images roughly Header metadata isn't copied over yet though. --- SEFramework/auxdir/basic.asdf | Bin 824 -> 0 bytes SEFramework/auxdir/with_primary.asdf | Bin 0 -> 1145 bytes SEFramework/tests/src/ASDF/AsdfFile_test.cpp | 31 ++++++++++++------- 3 files changed, 20 insertions(+), 11 deletions(-) delete mode 100644 SEFramework/auxdir/basic.asdf create mode 100644 SEFramework/auxdir/with_primary.asdf diff --git a/SEFramework/auxdir/basic.asdf b/SEFramework/auxdir/basic.asdf deleted file mode 100644 index b39b0f5a7a44fb60a25baf1bb9552f2bd27be268..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 824 zcma)5&2H2%5N@GRDX+j~s;cy|adxG$LOGF=QdLchME5|2s>qvoHDNl!6R_vzzc9;oH)yZgw$N-$KT9+{>IsjF5-7EVna6ccCDXYOrrBR`VdocOv$~= z=qw?BjeSB>rFGL%PBqud&>b5!!yB#z zL1-e}@&*f~Lzol-S|+i$7NuGUZ3a-NN-X$Xgs@*AF2W#~ODvix%T(13nkijzZ2kR# z2T=2h(01?uuBnGqbqEKR&7D_BgNm3COngRB+e2fcp86`ARswm>k#{T!vR0C2B{!xU z!d4|Lt#~bG!XP_y&QU+)pZZ42uA1)Cs?}YBDLoHb30F-`8xsk;2TP)#k@hnFPu4#= zqy~ANbDggVV%IaeYdb!^(OSI|)GE{yup=~Iy3GR^)#!{?VR>5|E*BKZ8pj8enJ!T% zoeNrgDG^ISIa6>8`SkkAV~ly<4wKW3;l|(J8yBwZtIOme{qo}bkLdp8*W}mF9=_?E rC(b@~cFWmk&OUc`+u8OHfWzed+fz6>kB>gCJ%4)Ro)2TXfaUia>;UQY diff --git a/SEFramework/auxdir/with_primary.asdf b/SEFramework/auxdir/with_primary.asdf new file mode 100644 index 0000000000000000000000000000000000000000..14622e356f114eaaf5b5092a45ac27f3af6c664d GIT binary patch literal 1145 zcmbVM!EVz)5G^-k@dLnRtEwu9coU~dTKiDp)B#0JT9o?GQYvc`Z){6(Hrm|;OqC;7 z;J`2N72H6aIdJ0xxFB)l%&eURlMsQ}gLh^;Z{Ez?)o$6R&I4rN5-w?5%JxoO``EEN z4g_lu+&;J4EwBx3+qEB}O(aQR2`PFJrvBI>!tWO~UDwfOk8xU5z7E-dRHiqGx}5OD zn*A26FP{=lAq5!G6&cHraSOR2MQJAvozMXrQ7#H7WHB9)fLdrLlyYPhi$NsAaToVk zoDuaIkBLm}oq~qYki?L!hzn=}H6k{&PzAa%uumtu`}YEKkC=U0F{fG>3CW-{fYE1? z4pp*H-h7lr)QhPkK9MBr8d7o)uQwn2<3sn>x4SQ@-|xS^_WsqAJ3s!2c4z3tE}q{`zTW)wA^v#q^3A(1S8BH9 wWh@5}GRnKFL4?YyL8Dqr9}iS-og5ya=CRW_%e@@%JZL)FH}~LcR?T$opTVa_00000 literal 0 HcmV?d00001 diff --git a/SEFramework/tests/src/ASDF/AsdfFile_test.cpp b/SEFramework/tests/src/ASDF/AsdfFile_test.cpp index 6fa200b2d..e319e41cc 100644 --- a/SEFramework/tests/src/ASDF/AsdfFile_test.cpp +++ b/SEFramework/tests/src/ASDF/AsdfFile_test.cpp @@ -29,6 +29,15 @@ using namespace SourceXtractor; + +struct AsdfFileFixture { + std::string primary_path; + + AsdfFileFixture() { + primary_path = Elements::getAuxiliaryPath("with_primary.asdf").native(); + } +}; + //----------------------------------------------------------------------------- BOOST_AUTO_TEST_SUITE (AsdfFile_test) @@ -36,9 +45,9 @@ BOOST_AUTO_TEST_SUITE (AsdfFile_test) //----------------------------------------------------------------------------- -BOOST_AUTO_TEST_CASE( open_file ) { +BOOST_FIXTURE_TEST_CASE( open_file, AsdfFileFixture ) { BOOST_CHECK_NO_THROW({ - AsdfFile asdf_file(Elements::getAuxiliaryPath("basic.asdf").native()); + AsdfFile asdf_file(primary_path); }); } @@ -52,12 +61,12 @@ BOOST_AUTO_TEST_CASE( missing_file ) { //----------------------------------------------------------------------------- -BOOST_AUTO_TEST_CASE( get_ndarray_by_index ) { - AsdfFile asdf_file(Elements::getAuxiliaryPath("basic.asdf").native()); +BOOST_FIXTURE_TEST_CASE( get_ndarray_by_index, AsdfFileFixture ) { + AsdfFile asdf_file(primary_path); auto ndarray = asdf_file.getNdarray(0); - BOOST_CHECK_EQUAL(ndarray->ndim(), 1); + BOOST_CHECK_EQUAL(ndarray->ndim(), 2); auto shape = ndarray->shape(); - std::vector expected_shape{8}; + std::vector expected_shape{1, 1}; BOOST_CHECK_EQUAL_COLLECTIONS( shape.begin(), shape.end(), expected_shape.begin(), expected_shape.end() @@ -66,12 +75,12 @@ BOOST_AUTO_TEST_CASE( get_ndarray_by_index ) { //----------------------------------------------------------------------------- -BOOST_AUTO_TEST_CASE( get_ndarray_by_path ) { - AsdfFile asdf_file(Elements::getAuxiliaryPath("basic.asdf").native()); - auto ndarray = asdf_file.getNdarray("data"); - BOOST_CHECK_EQUAL(ndarray->ndim(), 1); +BOOST_FIXTURE_TEST_CASE( get_ndarray_by_path, AsdfFileFixture ) { + AsdfFile asdf_file(primary_path); + auto ndarray = asdf_file.getNdarray("PRIMARY"); + BOOST_CHECK_EQUAL(ndarray->ndim(), 2); auto shape = ndarray->shape(); - std::vector expected_shape{8}; + std::vector expected_shape{1, 1}; BOOST_CHECK_EQUAL_COLLECTIONS( shape.begin(), shape.end(), expected_shape.begin(), expected_shape.end() From 813d64be87f801bd9092e52a7e2ab80a1ea01678 Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Fri, 3 Oct 2025 18:29:46 +0200 Subject: [PATCH 17/29] Add tests for AsdfImageSource approximately mirroring the applicable tests for FitsImageSource_test --- SEFramework/CMakeLists.txt | 3 + SEFramework/auxdir/multiple_hdu.asdf | Bin 0 -> 1459 bytes .../tests/src/ASDF/AsdfImageSource_test.cpp | 108 ++++++++++++++++++ 3 files changed, 111 insertions(+) create mode 100644 SEFramework/auxdir/multiple_hdu.asdf create mode 100644 SEFramework/tests/src/ASDF/AsdfImageSource_test.cpp diff --git a/SEFramework/CMakeLists.txt b/SEFramework/CMakeLists.txt index a4ef0f751..a7c8f39f4 100644 --- a/SEFramework/CMakeLists.txt +++ b/SEFramework/CMakeLists.txt @@ -206,6 +206,9 @@ if(ASDF_FOUND) elements_add_unit_test(AsdfFile_test tests/src/ASDF/AsdfFile_test.cpp LINK_LIBRARIES SEFramework TYPE Boost) + elements_add_unit_test(AsdfImageSource_test tests/src/ASDF/AsdfImageSource_test.cpp + LINK_LIBRARIES SEFramework + TYPE Boost) endif() #=============================================================================== # Declare the Python programs here diff --git a/SEFramework/auxdir/multiple_hdu.asdf b/SEFramework/auxdir/multiple_hdu.asdf new file mode 100644 index 0000000000000000000000000000000000000000..8cca55c7413b4cc6dcdd586f043084d9a0b891e0 GIT binary patch literal 1459 zcmbVM&ubGw6prXY7ykhjVM-}g%x8mMZkk$3(^lIEtu-u@>?FI8-3c?3)U{N4 z@gib9c<}1Me?ai$P4VnOJ?I}Gh~iNSUVO9p(QVQYaW3=b?fbs>%{Q|%Zf@jnAYDsp zNo8Ed)<)4>&713aNKQj?V#{19K&&ehMRO62A%SZK7rf$X#BLgx+tq}kswx_*P(~7R ztwJ>*WoXsiGQ%t|I-fbMiXF@dlt4$ci<`os%s@qlpioJkv`C#c2;*_&P>(cljTmUs z5u#xv5;a#i&9YXZUb|71G7pO|o{TF9`PhSM**K0`z~NHgKvSSb*KE>C&t9u#XK=Pk zrm}mA<8nbk8NkRrL3~-v<*K`tUe@x6z%~}Rtr|)q5UWzhoVS~5J&0y`*mtXh3vHzn zEOeXn;hW5DqhWucd+kkWl+_ZUr9Dy}u0e!67b&9K!|K{P%;qCsEH1EbV`fJTlJ@9( zA>+RjMesumURp617rKgMPlkie)j){TCab_%>CjeMD}n~Bs&$IR^c3XEfgqIG1nx%N z6{1cc#~mz-w)HsDOG-$e8pWR;&0ltF1LQZ?fqs2!6uUfzBlWH?W+4%gew;d4Qi^7- zFw&D$loeuikTTkK8^h1&#!(7Pk+-6K*h~5UMQ#Fc=pIrfMbk9pW3I3qOU6cwPnREm zIf%VF^L_u^^Ph*g2PdM=c1U{jtw)*W_@ChRKFfZ7{o-&ZEywA`+==L72i=?b{Y$~O z3m@Kj@2@|3_U7xZbfkai{RxEhRC+LpP--x#XQoe7($nAFPwVBLeLh&bGhe*;D=>1= + +#include + +#include "SEFramework/ASDF/AsdfImageSource.h" +#include "SEFramework/Image/ImageAccessor.h" + +using namespace SourceXtractor; + +struct AsdfImageSourceFixture { + std::string mhdu_path, primary_path; + + AsdfImageSourceFixture() { + mhdu_path = Elements::getAuxiliaryPath("multiple_hdu.asdf").native(); + primary_path = Elements::getAuxiliaryPath("with_primary.asdf").native(); + } +}; + +//----------------------------------------------------------------------------- + +BOOST_AUTO_TEST_SUITE(AsdfImageSource_test) + +//----------------------------------------------------------------------------- + +BOOST_FIXTURE_TEST_CASE(primary_test, AsdfImageSourceFixture) { + auto img_src = std::make_shared(primary_path, 0, ImageTile::FloatImage); + BOOST_CHECK_EQUAL(img_src->getHeight(), 1); + BOOST_CHECK_EQUAL(img_src->getWidth(), 1); + BOOST_CHECK_CLOSE(img_src->getImageTile(0, 0, 1, 1)->getValue(0, 0), 1024.44f, 1e-8); +} + +//----------------------------------------------------------------------------- + +// NOTE: This claims to test a compressed image? At least based on the name +// But the FITS file this test was based on does not contain a compressed image either... +BOOST_FIXTURE_TEST_CASE(compressed_image_test, AsdfImageSourceFixture) { + auto img_src = std::make_shared(mhdu_path, 0, ImageTile::FloatImage); + BOOST_CHECK_EQUAL(img_src->getHeight(), 1); + BOOST_CHECK_EQUAL(img_src->getWidth(), 1); + BOOST_CHECK_CLOSE(img_src->getImageTile(0, 0, 1, 1)->getValue(0, 0), 256.2f, 1e-8); +} + +//----------------------------------------------------------------------------- + +BOOST_FIXTURE_TEST_CASE(second_image_test, AsdfImageSourceFixture) { + auto img_src = std::make_shared(mhdu_path, 1, ImageTile::FloatImage); + BOOST_CHECK_EQUAL(img_src->getHeight(), 1); + BOOST_CHECK_EQUAL(img_src->getWidth(), 1); + BOOST_CHECK_CLOSE(img_src->getImageTile(0, 0, 1, 1)->getValue(0, 0), 1024.44f, 1e-8); +} + +//----------------------------------------------------------------------------- + +BOOST_FIXTURE_TEST_CASE(second_image_by_name_test, AsdfImageSourceFixture) { + auto img_src = std::make_shared(mhdu_path, "IMAGE2", ImageTile::FloatImage); + BOOST_CHECK_EQUAL(img_src->getHeight(), 1); + BOOST_CHECK_EQUAL(img_src->getWidth(), 1); + BOOST_CHECK_CLOSE(img_src->getImageTile(0, 0, 1, 1)->getValue(0, 0), 1024.44f, 1e-8); +} + +//----------------------------------------------------------------------------- + +BOOST_FIXTURE_TEST_CASE(ndarray_is_table_test, AsdfImageSourceFixture) { + BOOST_CHECK_THROW(std::make_shared(mhdu_path, 3), Elements::Exception); + BOOST_CHECK_THROW(std::make_shared(mhdu_path, "TABLE"), Elements::Exception); +} + +//----------------------------------------------------------------------------- + +BOOST_FIXTURE_TEST_CASE(bad_ndarray_test, AsdfImageSourceFixture) { + BOOST_CHECK_THROW(std::make_shared(mhdu_path, 5), Elements::Exception); +} + +//----------------------------------------------------------------------------- + +BOOST_FIXTURE_TEST_CASE(empty_ndarray_test, AsdfImageSourceFixture) { + BOOST_CHECK_THROW(std::make_shared(mhdu_path, 2), Elements::Exception); + BOOST_CHECK_THROW(std::make_shared(mhdu_path, "PRIMARY"), Elements::Exception); +} + +//----------------------------------------------------------------------------- + +BOOST_AUTO_TEST_SUITE_END() + +//----------------------------------------------------------------------------- + From 88ac0eecc557dd01eac4fb96e98a59a6ca22853e Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Mon, 6 Oct 2025 10:21:10 +0200 Subject: [PATCH 18/29] Replace the FitsReader::readFile convenience method with ImageFileReader::readImage on the base class--better name and more general Also fixes handling of image_type in FitsReader::get --- SEFramework/SEFramework/FITS/FitsReader.h | 8 ---- .../SEFramework/Image/ImageFileReader.h | 43 +++++++++++++++++++ SEFramework/src/lib/FITS/FitsReader.cpp | 2 +- .../tests/src/FITS/FitsReader_test.cpp | 4 +- .../lib/Configuration/SegmentationConfig.cpp | 2 +- 5 files changed, 47 insertions(+), 12 deletions(-) diff --git a/SEFramework/SEFramework/FITS/FitsReader.h b/SEFramework/SEFramework/FITS/FitsReader.h index aa4f130e8..b9a099048 100644 --- a/SEFramework/SEFramework/FITS/FitsReader.h +++ b/SEFramework/SEFramework/FITS/FitsReader.h @@ -23,7 +23,6 @@ #ifndef _SEFRAMEWORK_FITS_FITSREADER_H #define _SEFRAMEWORK_FITS_FITSREADER_H -#include "SEFramework/Image/BufferedImage.h" #include "SEFramework/Image/ImageFileReader.h" #include "SEFramework/Image/ImageTile.h" #include "SEFramework/FITS/FitsImageSource.h" @@ -60,13 +59,6 @@ class FitsReader : public ImageFileReader { std::shared_ptr getHdu(int hdu_num, ImageTile::ImageType image_type = ImageTile::AutoType); - // TODO: Integrate this interface into the ImageFileReader base class - template - static std::shared_ptr> readFile(const std::string& filename) { - auto image_source = std::make_shared(filename, 0, ImageTile::getTypeValue(T())); - return BufferedImage::create(image_source); - } - private: std::map m_image_hdu_map; diff --git a/SEFramework/SEFramework/Image/ImageFileReader.h b/SEFramework/SEFramework/Image/ImageFileReader.h index c90d77d9e..411b4864b 100644 --- a/SEFramework/SEFramework/Image/ImageFileReader.h +++ b/SEFramework/SEFramework/Image/ImageFileReader.h @@ -25,6 +25,7 @@ #include +#include "SEFramework/Image/BufferedImage.h" #include "SEFramework/Image/ImageSource.h" #include "SEFramework/Image/ImageTile.h" @@ -170,6 +171,48 @@ class ImageFileReader { */ static std::string getDefault(); + /** + * Convenience shortcut for returning a buffered image of a given datatype + * the file + * + * @param filename + * The filename of the file + * @param image_index (optional) + * The index of the image to return from the file if it contains multiple + * images (such as a multi-extension FITS file) + * @return + * A new instance of an Image. An Elements::Exception is thrown if the file type + * cannot be determined or does not contain a valid image. + */ + template + static std::shared_ptr> readImage(const std::string& filename, int image_index = 0) { + auto reader = create(filename); + auto image_source = reader->get(image_index, ImageTile::getTypeValue(T())); + return BufferedImage::create(image_source); + } + + /** + * Convenience shortcut for returning a buffered image of a given datatype + * the file + * + * @param filename + * The filename of the file + * @param image_path + * The name or path within the file for the image if there are multiple + * images in the file (e.g. the EXTNAME of a FITS HDU or a path within + * an ASDF file) + * @return + * A new instance of an Image. An Elements::Exception is thrown if the file type + * cannot be determined or does not contain a valid image. + */ + template + static std::shared_ptr> readImage(const std::string& filename, + const std::string& image_path) { + auto reader = create(filename); + auto image_source = reader->get(image_path, ImageTile::getTypeValue(T())); + return BufferedImage::create(image_source); + } + /** * Create an instance of an `ImageFileReader` from the filename (which may contain an optional * extension suffix) diff --git a/SEFramework/src/lib/FITS/FitsReader.cpp b/SEFramework/src/lib/FITS/FitsReader.cpp index 2845c6961..f608f0f05 100644 --- a/SEFramework/src/lib/FITS/FitsReader.cpp +++ b/SEFramework/src/lib/FITS/FitsReader.cpp @@ -63,7 +63,7 @@ std::shared_ptr FitsReader::get(int image_index, ImageTile::ImageTy while (known_index < image_index) { try { - auto image_source = getHdu(++hdu_num); + auto image_source = getHdu(++hdu_num, image_type); m_image_hdu_map[++known_index] = hdu_num; if (known_index == image_index) { diff --git a/SEFramework/tests/src/FITS/FitsReader_test.cpp b/SEFramework/tests/src/FITS/FitsReader_test.cpp index 56c0e5c68..c90e5dd32 100644 --- a/SEFramework/tests/src/FITS/FitsReader_test.cpp +++ b/SEFramework/tests/src/FITS/FitsReader_test.cpp @@ -51,7 +51,7 @@ BOOST_AUTO_TEST_SUITE (FitsReader_test) //----------------------------------------------------------------------------- BOOST_FIXTURE_TEST_CASE( read_file, FitsReaderFixture ) { - auto img = FitsReader::readFile(m_tmp_fits.path().native()); + auto img = FitsReader::readImage(m_tmp_fits.path().native()); BOOST_CHECK_EQUAL(img->getWidth(), 1); BOOST_CHECK_EQUAL(img->getHeight(), 1); BOOST_CHECK_EQUAL(img->getChunk(0, 0, 1, 1)->getValue(0, 0), 42); @@ -103,7 +103,7 @@ BOOST_AUTO_TEST_CASE ( open_with_extname ) { //----------------------------------------------------------------------------- BOOST_AUTO_TEST_CASE( missing_file ) { - BOOST_CHECK_THROW(FitsReader::readFile("/not/existing/path"), Elements::Exception); + BOOST_CHECK_THROW(FitsReader::readImage("/not/existing/path"), Elements::Exception); } //----------------------------------------------------------------------------- diff --git a/SEImplementation/src/lib/Configuration/SegmentationConfig.cpp b/SEImplementation/src/lib/Configuration/SegmentationConfig.cpp index 4eb6a4338..8dcc582a5 100644 --- a/SEImplementation/src/lib/Configuration/SegmentationConfig.cpp +++ b/SEImplementation/src/lib/Configuration/SegmentationConfig.cpp @@ -154,7 +154,7 @@ std::shared_ptr SegmentationConfig::loadFilter std::shared_ptr SegmentationConfig::loadFITSFilter(const std::string& filename) const { // read in the FITS file - auto convolution_kernel = FitsReader::readFile(filename); + auto convolution_kernel = FitsReader::readImage(filename); // give some feedback on the filter segConfigLogger.info() << "Loaded segmentation filter: " << filename << " height: " << convolution_kernel->getHeight() << " width: " << convolution_kernel->getWidth(); From cad4a1fe9d1a7ed791e4c0fcdea03867978e9f71 Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Mon, 6 Oct 2025 11:02:05 +0200 Subject: [PATCH 19/29] Add tests for AsdfReader mirroring FitsReader_test This caught (and now fixes) some bugs in how the ImageTile::ImageType was being passed through to AsdfImageSource::getImageTile --- SEFramework/CMakeLists.txt | 3 + .../SEFramework/ASDF/AsdfImageSource.h | 15 ++- SEFramework/src/lib/ASDF/AsdfImageSource.cpp | 2 +- SEFramework/tests/src/ASDF/1px.asdf.h | 86 ++++++++++++++ .../tests/src/ASDF/AsdfReader_test.cpp | 110 ++++++++++++++++++ 5 files changed, 213 insertions(+), 3 deletions(-) create mode 100644 SEFramework/tests/src/ASDF/1px.asdf.h create mode 100644 SEFramework/tests/src/ASDF/AsdfReader_test.cpp diff --git a/SEFramework/CMakeLists.txt b/SEFramework/CMakeLists.txt index a7c8f39f4..164581f2d 100644 --- a/SEFramework/CMakeLists.txt +++ b/SEFramework/CMakeLists.txt @@ -209,6 +209,9 @@ if(ASDF_FOUND) elements_add_unit_test(AsdfImageSource_test tests/src/ASDF/AsdfImageSource_test.cpp LINK_LIBRARIES SEFramework TYPE Boost) + elements_add_unit_test(AsdfReader_test tests/src/ASDF/AsdfReader_test.cpp + LINK_LIBRARIES SEFramework + TYPE Boost) endif() #=============================================================================== # Declare the Python programs here diff --git a/SEFramework/SEFramework/ASDF/AsdfImageSource.h b/SEFramework/SEFramework/ASDF/AsdfImageSource.h index 4c82f5b48..a482ef125 100644 --- a/SEFramework/SEFramework/ASDF/AsdfImageSource.h +++ b/SEFramework/SEFramework/ASDF/AsdfImageSource.h @@ -1,5 +1,5 @@ /** - * Copyright © 2019-2025 Université de Genève, LMU Munich - Faculty of Physics, IAP-CNRS/Sorbonne Université + * Copyright © 2019 Université de Genève, LMU Munich - Faculty of Physics, IAP-CNRS/Sorbonne Université * * This library is free software; you can redistribute it and/or modify it under * the terms of the GNU Lesser General Public License as published by the Free @@ -90,6 +90,10 @@ class AsdfImageSource : public ImageSource, public std::enable_shared_from_this< return m_filename; } + int getNdim() const { + return m_ndarray->ndim(); + } + /** * Returns the width of the image in pixels */ @@ -111,8 +115,15 @@ class AsdfImageSource : public ImageSource, public std::enable_shared_from_this< return m_ndarray->ndim() == 3 ? m_ndarray->shape().at(0) : 0; } + /** + * Returns the data type of the pixel values + * + * NOTE: Does not necessarily return the data type of the raw array, but + * rather the type specified when constructing the AsdfImageSource, to which + * the original data is converted. + */ ImageTile::ImageType getType() const override { - return m_ndarray->getImageType(); + return m_image_type; } void setLayer(int layer) override; diff --git a/SEFramework/src/lib/ASDF/AsdfImageSource.cpp b/SEFramework/src/lib/ASDF/AsdfImageSource.cpp index 49818017b..175e45c4e 100644 --- a/SEFramework/src/lib/ASDF/AsdfImageSource.cpp +++ b/SEFramework/src/lib/ASDF/AsdfImageSource.cpp @@ -108,7 +108,7 @@ void AsdfImageSource::setLayer(int layer) { std::shared_ptr AsdfImageSource::getImageTile(int x, int y, int width, int height) const { // TODO: (#5) support image data type conversion - auto tile = ImageTile::create(m_ndarray->getImageType(), x, y, width, height, + auto tile = ImageTile::create(m_image_type, x, y, width, height, std::const_pointer_cast(shared_from_this())); m_ndarray->fillImageTile(tile, m_current_layer); return tile; diff --git a/SEFramework/tests/src/ASDF/1px.asdf.h b/SEFramework/tests/src/ASDF/1px.asdf.h new file mode 100644 index 000000000..771ebc00d --- /dev/null +++ b/SEFramework/tests/src/ASDF/1px.asdf.h @@ -0,0 +1,86 @@ +/** Copyright © 2019 Université de Genève, LMU Munich - Faculty of Physics, IAP-CNRS/Sorbonne Université + * + * This library is free software; you can redistribute it and/or modify it under + * the terms of the GNU Lesser General Public License as published by the Free + * Software Foundation; either version 3.0 of the License, or (at your option) + * any later version. + * + * This library is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + */ +static const unsigned char image_asdf[] = { + 0x23, 0x41, 0x53, 0x44, 0x46, 0x20, 0x31, 0x2e, 0x30, 0x2e, 0x30, 0x0a, + 0x23, 0x41, 0x53, 0x44, 0x46, 0x5f, 0x53, 0x54, 0x41, 0x4e, 0x44, 0x41, + 0x52, 0x44, 0x20, 0x31, 0x2e, 0x36, 0x2e, 0x30, 0x0a, 0x25, 0x59, 0x41, + 0x4d, 0x4c, 0x20, 0x31, 0x2e, 0x31, 0x0a, 0x25, 0x54, 0x41, 0x47, 0x20, + 0x21, 0x20, 0x74, 0x61, 0x67, 0x3a, 0x73, 0x74, 0x73, 0x63, 0x69, 0x2e, + 0x65, 0x64, 0x75, 0x3a, 0x61, 0x73, 0x64, 0x66, 0x2f, 0x0a, 0x2d, 0x2d, + 0x2d, 0x20, 0x21, 0x63, 0x6f, 0x72, 0x65, 0x2f, 0x61, 0x73, 0x64, 0x66, + 0x2d, 0x31, 0x2e, 0x31, 0x2e, 0x30, 0x0a, 0x61, 0x73, 0x64, 0x66, 0x5f, + 0x6c, 0x69, 0x62, 0x72, 0x61, 0x72, 0x79, 0x3a, 0x20, 0x21, 0x63, 0x6f, + 0x72, 0x65, 0x2f, 0x73, 0x6f, 0x66, 0x74, 0x77, 0x61, 0x72, 0x65, 0x2d, + 0x31, 0x2e, 0x30, 0x2e, 0x30, 0x20, 0x7b, 0x61, 0x75, 0x74, 0x68, 0x6f, + 0x72, 0x3a, 0x20, 0x54, 0x68, 0x65, 0x20, 0x41, 0x53, 0x44, 0x46, 0x20, + 0x44, 0x65, 0x76, 0x65, 0x6c, 0x6f, 0x70, 0x65, 0x72, 0x73, 0x2c, 0x20, + 0x68, 0x6f, 0x6d, 0x65, 0x70, 0x61, 0x67, 0x65, 0x3a, 0x20, 0x27, 0x68, + 0x74, 0x74, 0x70, 0x3a, 0x2f, 0x2f, 0x67, 0x69, 0x74, 0x68, 0x75, 0x62, + 0x2e, 0x63, 0x6f, 0x6d, 0x2f, 0x61, 0x73, 0x64, 0x66, 0x2d, 0x66, 0x6f, + 0x72, 0x6d, 0x61, 0x74, 0x2f, 0x61, 0x73, 0x64, 0x66, 0x27, 0x2c, 0x0a, + 0x20, 0x20, 0x6e, 0x61, 0x6d, 0x65, 0x3a, 0x20, 0x61, 0x73, 0x64, 0x66, + 0x2c, 0x20, 0x76, 0x65, 0x72, 0x73, 0x69, 0x6f, 0x6e, 0x3a, 0x20, 0x34, + 0x2e, 0x31, 0x2e, 0x31, 0x2e, 0x64, 0x65, 0x76, 0x33, 0x38, 0x2b, 0x67, + 0x39, 0x37, 0x61, 0x39, 0x66, 0x65, 0x34, 0x39, 0x7d, 0x0a, 0x68, 0x69, + 0x73, 0x74, 0x6f, 0x72, 0x79, 0x3a, 0x0a, 0x20, 0x20, 0x65, 0x78, 0x74, + 0x65, 0x6e, 0x73, 0x69, 0x6f, 0x6e, 0x73, 0x3a, 0x0a, 0x20, 0x20, 0x2d, + 0x20, 0x21, 0x63, 0x6f, 0x72, 0x65, 0x2f, 0x65, 0x78, 0x74, 0x65, 0x6e, + 0x73, 0x69, 0x6f, 0x6e, 0x5f, 0x6d, 0x65, 0x74, 0x61, 0x64, 0x61, 0x74, + 0x61, 0x2d, 0x31, 0x2e, 0x30, 0x2e, 0x30, 0x0a, 0x20, 0x20, 0x20, 0x20, + 0x65, 0x78, 0x74, 0x65, 0x6e, 0x73, 0x69, 0x6f, 0x6e, 0x5f, 0x63, 0x6c, + 0x61, 0x73, 0x73, 0x3a, 0x20, 0x61, 0x73, 0x64, 0x66, 0x2e, 0x65, 0x78, + 0x74, 0x65, 0x6e, 0x73, 0x69, 0x6f, 0x6e, 0x2e, 0x5f, 0x6d, 0x61, 0x6e, + 0x69, 0x66, 0x65, 0x73, 0x74, 0x2e, 0x4d, 0x61, 0x6e, 0x69, 0x66, 0x65, + 0x73, 0x74, 0x45, 0x78, 0x74, 0x65, 0x6e, 0x73, 0x69, 0x6f, 0x6e, 0x0a, + 0x20, 0x20, 0x20, 0x20, 0x65, 0x78, 0x74, 0x65, 0x6e, 0x73, 0x69, 0x6f, + 0x6e, 0x5f, 0x75, 0x72, 0x69, 0x3a, 0x20, 0x61, 0x73, 0x64, 0x66, 0x3a, + 0x2f, 0x2f, 0x61, 0x73, 0x64, 0x66, 0x2d, 0x66, 0x6f, 0x72, 0x6d, 0x61, + 0x74, 0x2e, 0x6f, 0x72, 0x67, 0x2f, 0x63, 0x6f, 0x72, 0x65, 0x2f, 0x65, + 0x78, 0x74, 0x65, 0x6e, 0x73, 0x69, 0x6f, 0x6e, 0x73, 0x2f, 0x63, 0x6f, + 0x72, 0x65, 0x2d, 0x31, 0x2e, 0x36, 0x2e, 0x30, 0x0a, 0x20, 0x20, 0x20, + 0x20, 0x6d, 0x61, 0x6e, 0x69, 0x66, 0x65, 0x73, 0x74, 0x5f, 0x73, 0x6f, + 0x66, 0x74, 0x77, 0x61, 0x72, 0x65, 0x3a, 0x20, 0x21, 0x63, 0x6f, 0x72, + 0x65, 0x2f, 0x73, 0x6f, 0x66, 0x74, 0x77, 0x61, 0x72, 0x65, 0x2d, 0x31, + 0x2e, 0x30, 0x2e, 0x30, 0x20, 0x7b, 0x6e, 0x61, 0x6d, 0x65, 0x3a, 0x20, + 0x61, 0x73, 0x64, 0x66, 0x5f, 0x73, 0x74, 0x61, 0x6e, 0x64, 0x61, 0x72, + 0x64, 0x2c, 0x20, 0x76, 0x65, 0x72, 0x73, 0x69, 0x6f, 0x6e, 0x3a, 0x20, + 0x31, 0x2e, 0x31, 0x2e, 0x31, 0x7d, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x73, + 0x6f, 0x66, 0x74, 0x77, 0x61, 0x72, 0x65, 0x3a, 0x20, 0x21, 0x63, 0x6f, + 0x72, 0x65, 0x2f, 0x73, 0x6f, 0x66, 0x74, 0x77, 0x61, 0x72, 0x65, 0x2d, + 0x31, 0x2e, 0x30, 0x2e, 0x30, 0x20, 0x7b, 0x6e, 0x61, 0x6d, 0x65, 0x3a, + 0x20, 0x61, 0x73, 0x64, 0x66, 0x2c, 0x20, 0x76, 0x65, 0x72, 0x73, 0x69, + 0x6f, 0x6e, 0x3a, 0x20, 0x34, 0x2e, 0x31, 0x2e, 0x31, 0x2e, 0x64, 0x65, + 0x76, 0x33, 0x38, 0x2b, 0x67, 0x39, 0x37, 0x61, 0x39, 0x66, 0x65, 0x34, + 0x39, 0x7d, 0x0a, 0x50, 0x52, 0x49, 0x4d, 0x41, 0x52, 0x59, 0x3a, 0x20, + 0x21, 0x63, 0x6f, 0x72, 0x65, 0x2f, 0x6e, 0x64, 0x61, 0x72, 0x72, 0x61, + 0x79, 0x2d, 0x31, 0x2e, 0x31, 0x2e, 0x30, 0x0a, 0x20, 0x20, 0x73, 0x6f, + 0x75, 0x72, 0x63, 0x65, 0x3a, 0x20, 0x30, 0x0a, 0x20, 0x20, 0x64, 0x61, + 0x74, 0x61, 0x74, 0x79, 0x70, 0x65, 0x3a, 0x20, 0x66, 0x6c, 0x6f, 0x61, + 0x74, 0x36, 0x34, 0x0a, 0x20, 0x20, 0x62, 0x79, 0x74, 0x65, 0x6f, 0x72, + 0x64, 0x65, 0x72, 0x3a, 0x20, 0x62, 0x69, 0x67, 0x0a, 0x20, 0x20, 0x73, + 0x68, 0x61, 0x70, 0x65, 0x3a, 0x20, 0x5b, 0x31, 0x2c, 0x20, 0x31, 0x5d, + 0x0a, 0x2e, 0x2e, 0x2e, 0x0a, 0xd3, 0x42, 0x4c, 0x4b, 0x00, 0x30, 0x00, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, + 0x00, 0x00, 0x08, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x08, 0x00, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x08, 0x15, 0x57, 0x0f, 0x33, 0xdd, + 0x31, 0xaf, 0xec, 0x2c, 0x53, 0x93, 0x50, 0x93, 0xcb, 0x36, 0x1d, 0x40, + 0x45, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x23, 0x41, 0x53, 0x44, 0x46, + 0x20, 0x42, 0x4c, 0x4f, 0x43, 0x4b, 0x20, 0x49, 0x4e, 0x44, 0x45, 0x58, + 0x0a, 0x25, 0x59, 0x41, 0x4d, 0x4c, 0x20, 0x31, 0x2e, 0x31, 0x0a, 0x2d, + 0x2d, 0x2d, 0x0a, 0x2d, 0x20, 0x37, 0x30, 0x31, 0x0a, 0x2e, 0x2e, 0x2e +}; +static const size_t image_asdf_len = 804; diff --git a/SEFramework/tests/src/ASDF/AsdfReader_test.cpp b/SEFramework/tests/src/ASDF/AsdfReader_test.cpp new file mode 100644 index 000000000..98aee4b97 --- /dev/null +++ b/SEFramework/tests/src/ASDF/AsdfReader_test.cpp @@ -0,0 +1,110 @@ +/** Copyright © 2019 Université de Genève, LMU Munich - Faculty of Physics, IAP-CNRS/Sorbonne Université + * + * This library is free software; you can redistribute it and/or modify it under + * the terms of the GNU Lesser General Public License as published by the Free + * Software Foundation; either version 3.0 of the License, or (at your option) + * any later version. + * + * This library is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + */ +/** + * @file tests/src/ASDF/AsdfReader_test.cpp + * @date 10/06/25 + * @author embray + */ + +#include +#include + +#include +#include +#include + +#include "SEFramework/ASDF/AsdfImageSource.h" +#include "SEFramework/ASDF/AsdfReader.h" +#include "SEFramework/Image/ImageFileReader.h" + +#include "1px.asdf.h" + +using namespace SourceXtractor; + +struct AsdfReaderFixture { + Elements::TempFile m_tmp_asdf; + + AsdfReaderFixture() { + std::ofstream out{m_tmp_asdf.path().c_str()}; + out.write(reinterpret_cast(image_asdf), image_asdf_len); + } +}; + +//----------------------------------------------------------------------------- + +BOOST_AUTO_TEST_SUITE (AsdfReader_test) + + +//----------------------------------------------------------------------------- + +BOOST_FIXTURE_TEST_CASE( read_file, AsdfReaderFixture ) { + auto img = AsdfReader::readImage(m_tmp_asdf.path().native()); + BOOST_CHECK_EQUAL(img->getWidth(), 1); + BOOST_CHECK_EQUAL(img->getHeight(), 1); + BOOST_CHECK_EQUAL(img->getChunk(0, 0, 1, 1)->getValue(0, 0), 42); +} + +//----------------------------------------------------------------------------- + +BOOST_FIXTURE_TEST_CASE ( image_source, AsdfReaderFixture ) { + auto img = AsdfImageSource(m_tmp_asdf.path().native()); + BOOST_CHECK_EQUAL(img.getNdim(), 2); +} + +//----------------------------------------------------------------------------- + +BOOST_AUTO_TEST_CASE ( detect_file_type ) { + auto reader = ImageFileReader::create(Elements::getAuxiliaryPath("with_primary.asdf").native()); + auto* asdf_reader = dynamic_cast(reader.get()); + BOOST_CHECK(asdf_reader != nullptr); +} + +//----------------------------------------------------------------------------- + +BOOST_AUTO_TEST_CASE ( iterate ) { + auto reader = ImageFileReader::create(Elements::getAuxiliaryPath("with_primary.asdf").native()); + // Should just return one image, from the primary HDU (the next HDU in this file is a table) + int image_count = 0; + for (const auto& img_source: *reader) { + auto asdf_source = std::dynamic_pointer_cast(img_source); + BOOST_CHECK(asdf_source != nullptr); + image_count++; + } + BOOST_CHECK_EQUAL(image_count, 1); +} + +//----------------------------------------------------------------------------- + +BOOST_AUTO_TEST_CASE ( open_with_extname ) { + auto reader = ImageFileReader::create(Elements::getAuxiliaryPath("multiple_hdu.asdf").native()); + auto* asdf_reader = dynamic_cast(reader.get()); + BOOST_CHECK(asdf_reader != nullptr); + auto img = reader->get("IMAGE2"); + BOOST_CHECK_EQUAL(img->getWidth(), 1); + BOOST_CHECK_EQUAL(img->getHeight(), 1); +} + +//----------------------------------------------------------------------------- + +BOOST_AUTO_TEST_CASE( missing_file ) { + BOOST_CHECK_THROW(AsdfReader::readImage("/not/existing/path"), Elements::Exception); +} + +//----------------------------------------------------------------------------- + + +BOOST_AUTO_TEST_SUITE_END () From 94e132902ca9d48a9f7d02e39e28a189575e3caa Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Mon, 27 Oct 2025 16:53:20 +0100 Subject: [PATCH 20/29] refactor: allow WCS() to be instantiated from any ImageSource (in theory) which will make it easier to avoid having to write different special cases for FITS and ASDF files specifically (even though they're the only applicable cases for now) For the ASDF case specifically support extracting a FitsWCS. Right now just provides a wrapper around the relevant libasdf GWCS extension types--still need to then plug that into WCSLIB and add tests --- SEFramework/CMakeLists.txt | 5 +- SEFramework/SEFramework/ASDF/AsdfFile.h | 63 ++++++++++ .../SEFramework/CoordinateSystem/WCS.h | 10 +- SEFramework/auxdir/wcs_header.asdf | Bin 0 -> 2773 bytes SEFramework/src/lib/ASDF/AsdfFile.cpp | 111 +++++++++++++++++- SEFramework/src/lib/CoordinateSystem/WCS.cpp | 40 ++++++- SEFramework/tests/src/ASDF/AsdfFile_test.cpp | 43 +++++++ .../Configuration/DetectionImageConfig.h | 7 ++ .../Configuration/DetectionImageConfig.cpp | 11 +- .../Configuration/MeasurementImageConfig.cpp | 10 +- 10 files changed, 276 insertions(+), 24 deletions(-) create mode 100644 SEFramework/auxdir/wcs_header.asdf diff --git a/SEFramework/CMakeLists.txt b/SEFramework/CMakeLists.txt index 164581f2d..0d3172eee 100644 --- a/SEFramework/CMakeLists.txt +++ b/SEFramework/CMakeLists.txt @@ -31,6 +31,8 @@ elements_depends_on_subdirs(ModelFitting) find_package(GMock) find_package(CCfits) find_package(BoostDLL) +# TODO: For ASDF support we also need optional GWCS support to work with WCS +# from ASDF files; this needs to be checked as well find_package(ASDF) find_package(FFTW COMPONENTS single double) find_package(WCSLIB REQUIRED) @@ -40,9 +42,10 @@ find_package(WCSLIB REQUIRED) # ASDF support is optional--don't compile it if libasdf was not found #=============================================================================== if(NOT ASDF_FOUND) - message("libasdf not found, compiling without ASDF support!") + message("libasdf not found, compiling without ASDF support") else() file(GLOB ASDF_SRC src/lib/ASDF/*.cpp) + add_definitions(-DWITH_ASDF) endif() diff --git a/SEFramework/SEFramework/ASDF/AsdfFile.h b/SEFramework/SEFramework/ASDF/AsdfFile.h index b7594922c..df8651993 100644 --- a/SEFramework/SEFramework/ASDF/AsdfFile.h +++ b/SEFramework/SEFramework/ASDF/AsdfFile.h @@ -29,6 +29,7 @@ #include #include +#include #include @@ -124,6 +125,58 @@ class AsdfFile { mutable std::unique_ptr> m_shape; }; + /** + * Wrapper around asdf_gwcs_fits_t, which represents a FITS- (WCSLIB) + * compatible WCS for use with SourceXtractor::WCS + */ + class FitsWCS { + public: + friend class AsdfFile; + + ~FitsWCS() { + asdf_gwcs_fits_destroy(m_gwcs_fits_ptr); + } + + std::array crpix() const noexcept { + return {m_gwcs_fits_ptr->crpix[0], m_gwcs_fits_ptr->crpix[1]}; + } + + std::array crval() const noexcept { + return {m_gwcs_fits_ptr->crval[0], m_gwcs_fits_ptr->crval[1]}; + } + + std::array cdelt() const noexcept { + return {m_gwcs_fits_ptr->cdelt[0], m_gwcs_fits_ptr->cdelt[1]}; + } + + std::array, 2> pc() const noexcept { + return {{ + {m_gwcs_fits_ptr->pc[0][0], m_gwcs_fits_ptr->pc[0][1]}, + {m_gwcs_fits_ptr->pc[1][0], m_gwcs_fits_ptr->pc[1][1]} + }}; + } + + std::array ctype() const noexcept { + return {{ + m_gwcs_fits_ptr->ctype[0] ? std::string_view(m_gwcs_fits_ptr->ctype[0]) : std::string_view(), + m_gwcs_fits_ptr->ctype[1] ? std::string_view(m_gwcs_fits_ptr->ctype[1]) : std::string_view() + }}; + } + + private: + /** + * Private constructor for creating the `Ndarray` wrapper from a raw asdf_value_t * + */ + explicit FitsWCS(const AsdfFile& file, asdf_value_t *ptr); + /** + * Private constructor for creating the `Ndarray` wrapper from a raw asdf_ndarray_t * + */ + explicit FitsWCS(asdf_gwcs_fits_t *ptr) + : m_gwcs_fits_ptr(ptr) {} + + asdf_gwcs_fits_t* m_gwcs_fits_ptr; + }; + /** * Return the N-th ndarray from the top-level of the ASDF tree iterating the top-level * keys in order. @@ -141,6 +194,16 @@ class AsdfFile { /* TODO: More general methods for reading metadata from the ASDF tree; for the first version * not needed though. */ + /** + * Return FITS-compatible WCS metadata (wrapped in `AsdfFile::FitsWCS`) + * + * If called without any arguments it will look for the first applicable GWCS + * object in the ASDF metadata tree. Called with a path argument it will + * look for one specifically at that path. In either case if no matching + * GWCS is found will return a null-ish pointer. + */ + std::unique_ptr getFitsWCS(); + std::unique_ptr getFitsWCS(const std::string& path); private: struct AsdfValueDestroy { diff --git a/SEFramework/SEFramework/CoordinateSystem/WCS.h b/SEFramework/SEFramework/CoordinateSystem/WCS.h index 5f0524e42..352e8d461 100644 --- a/SEFramework/SEFramework/CoordinateSystem/WCS.h +++ b/SEFramework/SEFramework/CoordinateSystem/WCS.h @@ -29,6 +29,10 @@ #include "SEFramework/CoordinateSystem/CoordinateSystem.h" #include "SEFramework/FITS/FitsImageSource.h" +#include "SEFramework/Image/ImageSource.h" +#ifdef WITH_ASDF +#include "SEFramework/ASDF/AsdfImageSource.h" +#endif struct wcsprm; @@ -37,6 +41,10 @@ namespace SourceXtractor { class WCS : public CoordinateSystem { public: explicit WCS(const FitsImageSource& fits_image_source); +#ifdef WITH_ASDF + explicit WCS(const AsdfImageSource& asdf_image_source); +#endif + explicit WCS(const ImageSource& image_source); explicit WCS(const WCS& original); virtual ~WCS(); @@ -55,7 +63,7 @@ class WCS : public CoordinateSystem { explicit WCS(std::unique_ptr> wcs) : m_wcs(std::move(wcs)) {} - void init(char* headers, int number_of_records); + void initFits(char* headers, int number_of_records); std::unique_ptr> m_wcs; }; diff --git a/SEFramework/auxdir/wcs_header.asdf b/SEFramework/auxdir/wcs_header.asdf new file mode 100644 index 0000000000000000000000000000000000000000..8a122faccc5956b77c2d2643f2a5975a9bc6a657 GIT binary patch literal 2773 zcmb_eJ&YSg6g~)0m;)gp5apbr#0BeJ+h_aE798iqxhTpdpxhOpQ?wfIj_rZn@**1@Q;u2;AM8G53!HF0|tc5lQ1gS_m zWSaz0Oc=M&AR3Yw_el#C2SUUx$LUisNV>Kc4Yi|Q#D-X?bX>OFOcL+k(MFwGfxCYi~ zFz@5ya)Co1Gk**q2~TBLS)AR_nEn)iC1NHZ>5}0z>kNc=lrRrURRgb5PNXrs^@0eC z#VYW+sUQ*a3A9E)g$M}n_y9|xt&4?A9ss`Z%A=jw2d)D@@KLkNE1;mWt}d11MtxLU zP8&^N#&kpi*Vw9;%Vi76|BKRtAhBtIN!H~)m|xW95$ZA3ZxbU?dsQFAt416lJR;nc zVL|*`BMYU{t9;Z(mhcutRN3F4+>Hk*rydSm36oSUdI>b|ur(@L$gI<%(b0MPL=X=; zJoAu*RDkA!q>lW9$m9`JcvxsMUt1Imhg^2$6wrHAZ~)$=L)@ofUmKPMkWV*CLLazp zHyVL{Zp(?8Jmw;G#w9>kjt;Gg5bcLBR8-A#RwxMD#(}DnVS?qYQrWIssMc&Y?Ycoi z?!)YvvX|k#mbvqM5&%ScjvZ`R?a8$5euM{_w;>x^Ka*NGvywG()r9FWsKRWXE4lqc z@(fVU=2;#~{TzyC8peq@cnqRMWZWPMEfip3Ad0=&qFg+jC@qv12m17FI7sR$8EaqY z%wZH#4>}g+!M<0@j4;%+%HK)zt%(^Uw4ZZ4pcWDXxNKzBa1H(68RU@wZXgPSLro(H zb=oXC2={-SM{&e$@~+LWg|ekjJX1P7j;P0Y1_dKM!bs>59TpT;C5pL#k&w!Z4?Zi@ z8ZJCIq%_2Wa7X399cf692{U4{LF{5780{t?+d@0LnX_qg&nNu}QxCX0?0inOP?!Wk z!M5$f?Gx)~=gM=BKF?*ZIeqGG{rMM~Uw`_{M}Ob@+26yT&mXz>(v$aB44uu&H_zVK ze9^gb?Z~T-uY7j7y(zt1|Khtu5t)>d~FBJ?F5dZ)H literal 0 HcmV?d00001 diff --git a/SEFramework/src/lib/ASDF/AsdfFile.cpp b/SEFramework/src/lib/ASDF/AsdfFile.cpp index c001534b4..3c70870f6 100644 --- a/SEFramework/src/lib/ASDF/AsdfFile.cpp +++ b/SEFramework/src/lib/ASDF/AsdfFile.cpp @@ -102,14 +102,69 @@ std::unique_ptr AsdfFile::getNdarray(int index) { << m_path.native(); } + std::unique_ptr AsdfFile::getNdarray(const std::string &path) { AsdfValuePtr value = getValue(path); return std::unique_ptr(new Ndarray(*this, value.get())); } +bool fitsWcsValuePredicate(asdf_value_t* value) { + asdf_gwcs_t* gwcs = nullptr; + asdf_value_err_t err = asdf_value_as_gwcs(value, &gwcs); + + if (ASDF_VALUE_OK != err || !gwcs) { + return false; + } + + return asdf_gwcs_is_fits((asdf_file_t*)asdf_value_file(value), gwcs); +} + + +std::unique_ptr AsdfFile::getFitsWCS() { + AsdfValuePtr root = getValue("/"); + + if (!root) { + throw AsdfValueNotFoundException() << "Could not load the root of the ASDF tree in file: " + << m_path.native(); + } + + // Find the first applicable GWCS, if any + asdf_find_item_t* item = asdf_value_find(root.get(), fitsWcsValuePredicate); + + if (!item) { + logger.warn() << "No FITS-compatible WCS could be found in the ASDF file: " + << m_path.native(); + return std::unique_ptr{}; + } + + asdf_value_t* value = asdf_find_item_value(item); + return std::unique_ptr(new FitsWCS(*this, value)); +} + + +std::unique_ptr AsdfFile::getFitsWCS(const std::string& path) { + asdf_value_t* value = asdf_get_value(m_asdf_ptr.get(), path.c_str()); + asdf_gwcs_t* gwcs = nullptr; + + if (!value) { + throw AsdfValueNotFoundException() << "No value at given WCS path " << path + << " in ASDF file: " << m_path.native(); + } + + asdf_value_err_t err = asdf_value_as_gwcs(value, &gwcs); + + if (ASDF_VALUE_OK != err || !asdf_gwcs_is_fits(m_asdf_ptr.get(), gwcs)) { + throw AsdfValueNotFoundException() << "Value at given WCS path " << path + << " is not a FITS-compatible GWCS in ASDF file: " << m_path.native(); + } + + return std::unique_ptr(new FitsWCS(*this, value)); +} + + AsdfFile::AsdfValuePtr AsdfFile::getValue(const std::string& path) { - asdf_value_t *v = asdf_get_value(m_asdf_ptr.get(), path.c_str()); + asdf_value_t* v = asdf_get_value(m_asdf_ptr.get(), path.c_str()); if (!v) { throw AsdfValueNotFoundException() << "No value at path " << path << " in file: " << m_path.native(); @@ -119,7 +174,7 @@ AsdfFile::AsdfValuePtr AsdfFile::getValue(const std::string& path) { } -static ImageTile::ImageType convertDatatypeToImageType(const asdf_datatype_t *datatype) { +static ImageTile::ImageType convertDatatypeToImageType(const asdf_datatype_t* datatype) { ImageTile::ImageType image_type; switch (datatype->type) { @@ -178,9 +233,9 @@ static asdf_scalar_datatype_t convertImageTypeToDatatype(ImageTile::ImageType im } -AsdfFile::Ndarray::Ndarray(const AsdfFile& file, asdf_value_t *value) +AsdfFile::Ndarray::Ndarray(const AsdfFile& file, asdf_value_t* value) : Ndarray((asdf_ndarray_t*)nullptr) { - asdf_ndarray_t *ndarray_ptr = nullptr; + asdf_ndarray_t* ndarray_ptr = nullptr; asdf_value_err_t err = asdf_value_as_ndarray(value, &ndarray_ptr); switch (err) { case ASDF_VALUE_OK: @@ -192,7 +247,7 @@ AsdfFile::Ndarray::Ndarray(const AsdfFile& file, asdf_value_t *value) << file.m_path.native(); } default: { - const char *error_message = asdf_error(file.getAsdfFilePtr()); + const char* error_message = asdf_error(file.getAsdfFilePtr()); throw AsdfValueNotFoundException() << "An error occurred reading the ASDF file " << file.m_path.native() << ": " << error_message; } @@ -249,4 +304,50 @@ bool AsdfFile::Ndarray::isSupportedImage() const { return true; } + +/** + * NOTE: The asdf_value_t* should be for the full GWCS object, not just the + * fitswcs_imaging part + * + * The full GWCS is needed in order to properly read the FITS WCS out of it. + */ +AsdfFile::FitsWCS::FitsWCS(const AsdfFile& file, asdf_value_t* value) + : FitsWCS((asdf_gwcs_fits_t*)nullptr) { + asdf_gwcs_t* gwcs_ptr = nullptr; + asdf_gwcs_fits_t* gwcs_fits_ptr = nullptr; + asdf_value_err_t err = asdf_value_as_gwcs(value, &gwcs_ptr); + switch (err) { + case ASDF_VALUE_OK: + // Value exists and is an ndarray: OK + break; + case ASDF_VALUE_ERR_TYPE_MISMATCH: { + const char* path = asdf_value_path(value); + throw AsdfValueTypeMismatchException() << "Value at " << path << " is not a GWCS: " + << file.m_path.native(); + } + default: { + const char* error_message = asdf_error(file.getAsdfFilePtr()); + throw AsdfValueNotFoundException() << "An error occurred reading the ASDF file " + << file.m_path.native() << ": " << error_message; + } + } + + if (!asdf_gwcs_is_fits(file.getAsdfFilePtr(), gwcs_ptr)) { + const char* path = asdf_value_path(value); + throw AsdfValueTypeMismatchException() << "Value at " << path << " does not contain " + "a FITS-compatible WCS: " << file.m_path.native(); + } + + // This structure is already checked by asdf_gwcs_is_fits so we should expect all + // these values to be valid now... + const asdf_gwcs_step_t* step0 = &gwcs_ptr->steps[0]; + gwcs_fits_ptr = (asdf_gwcs_fits_t*)step0->transform; + assert(gwcs_fits_ptr); + + m_gwcs_fits_ptr = gwcs_fits_ptr; + // We don't need the asdf_value_t anymore at this point and can release it. + asdf_value_destroy(value); +} + + } // namespace SourceXtractor diff --git a/SEFramework/src/lib/CoordinateSystem/WCS.cpp b/SEFramework/src/lib/CoordinateSystem/WCS.cpp index fd5997701..b07298655 100644 --- a/SEFramework/src/lib/CoordinateSystem/WCS.cpp +++ b/SEFramework/src/lib/CoordinateSystem/WCS.cpp @@ -33,6 +33,11 @@ #include #include +#ifdef WITH_ASDF +#include +#include +#endif + #include "ElementsKernel/Exception.h" #include "ElementsKernel/Logging.h" @@ -156,7 +161,7 @@ WCS::WCS(const FitsImageSource& fits_image_source) : m_wcs(nullptr, nullptr) { int number_of_records = 0; auto fits_headers = fits_image_source.getFitsHeaders(number_of_records); - init(&(*fits_headers)[0], number_of_records); + initFits(&(*fits_headers)[0], number_of_records); } WCS::WCS(const WCS& original) : m_wcs(nullptr, nullptr) { @@ -164,6 +169,10 @@ WCS::WCS(const WCS& original) : m_wcs(nullptr, nullptr) { //FIXME Horrible hack: I couldn't figure out how to properly do a deep copy wcsprm so instead // of making a copy, I use the ascii headers output from the original to recreate a new one + // (embray): Major sympathies here. Looking through the wcslib headers I found there is a + // wcscopy() which is just a wrapper around wcssub() which should do it. I'll give that a try + // later. + int number_of_records; char *raw_header; @@ -171,13 +180,13 @@ WCS::WCS(const WCS& original) : m_wcs(nullptr, nullptr) { throw Elements::Exception() << "Failed to get the FITS headers for the WCS coordinate system when copying WCS"; } - init(raw_header, number_of_records); + initFits(raw_header, number_of_records); free(raw_header); } -void WCS::init(char* headers, int number_of_records) { +void WCS::initFits(char* headers, int number_of_records) { wcserr_enable(1); int nreject = 0, nwcs = 0, nreject_strict = 0; @@ -216,6 +225,31 @@ void WCS::init(char* headers, int number_of_records) { } +#ifdef HAVE_ASDF +/** WCS initializer from an ASDF file + * + * Currently this makes a brash assumption: if there is any compatible GWCS + * object in the file it "must" be the right one. This assumption can be wrong + * but in practice most ASDF files have one data array, one WCS. + * + * Later we will figure out how to work in some config option(s) to explicitly + * provide a path to the correct WCS to use if there is any ambiguity. + */ +WCS::WCS(const AsdfImageSource& fits_image_source) : m_wcs(nullptr, nullptr) { +} +#endif + + +/** + * Initializer for a generic ImageSource + * + * This just creates a dummy identity WCS and logs a warning + */ +WCS::WCS(const ImageSource &) : WCS(identity(2)) { + logger.warn() << "No WCS info on generic image source; creating an identity WCS"; +} + + WCS::~WCS() { } diff --git a/SEFramework/tests/src/ASDF/AsdfFile_test.cpp b/SEFramework/tests/src/ASDF/AsdfFile_test.cpp index e319e41cc..a9c38f617 100644 --- a/SEFramework/tests/src/ASDF/AsdfFile_test.cpp +++ b/SEFramework/tests/src/ASDF/AsdfFile_test.cpp @@ -32,9 +32,11 @@ using namespace SourceXtractor; struct AsdfFileFixture { std::string primary_path; + std::string wcs_path; AsdfFileFixture() { primary_path = Elements::getAuxiliaryPath("with_primary.asdf").native(); + wcs_path = Elements::getAuxiliaryPath("wcs_header.asdf").native(); } }; @@ -87,6 +89,47 @@ BOOST_FIXTURE_TEST_CASE( get_ndarray_by_path, AsdfFileFixture ) { ); } + +BOOST_FIXTURE_TEST_CASE( get_fits_wcs_auto, AsdfFileFixture ) { + AsdfFile asdf_file(wcs_path); + auto wcs = asdf_file.getFitsWCS(); + + BOOST_CHECK(wcs != nullptr); + + auto crpix = wcs->crpix(); + std::array expected_crpix{12099.5, -88700.5}; + for (int idx = 0; idx < 2; idx++) { + BOOST_CHECK_CLOSE(crpix[idx], expected_crpix[idx], 1e-6); + } + + auto crval = wcs->crval(); + std::array expected_crval{270., 64.60237301}; + for (int idx = 0; idx < 2; idx++) { + BOOST_CHECK_CLOSE(crval[idx], expected_crval[idx], 1e-6); + } + + auto cdelt = wcs->cdelt(); + std::array expected_cdelt{1.52777778e-05, 1.52777778e-05}; + for (int idx = 0; idx < 2; idx++) { + BOOST_CHECK_CLOSE(cdelt[idx], expected_cdelt[idx], 1e-6); + } + + auto pc = wcs->pc(); + std::array, 2> expected_pc{{{1., 0.}, {-0., 1.}}}; + for (int idx = 0; idx < 2; idx++) { + for (int jdx = 0; jdx < 2; jdx++) { + BOOST_CHECK_CLOSE(pc[idx][jdx], expected_pc[idx][jdx], 1e-6); + } + } + + auto ctype = wcs->ctype(); + std::array expected_ctype{{"RA---TAN", "DEC--TAN"}}; + BOOST_CHECK_EQUAL_COLLECTIONS( + ctype.begin(), ctype.end(), + expected_ctype.begin(), expected_ctype.end() + ); +} + //----------------------------------------------------------------------------- diff --git a/SEImplementation/SEImplementation/Configuration/DetectionImageConfig.h b/SEImplementation/SEImplementation/Configuration/DetectionImageConfig.h index 90efe9cfa..7d2a5270b 100644 --- a/SEImplementation/SEImplementation/Configuration/DetectionImageConfig.h +++ b/SEImplementation/SEImplementation/Configuration/DetectionImageConfig.h @@ -28,6 +28,9 @@ #include "SEFramework/Image/Image.h" #include "SEFramework/Image/ImageSource.h" #include "SEFramework/CoordinateSystem/CoordinateSystem.h" +#ifdef WITH_ASDF +#include "SEFramework/ASDF/AsdfImageSource.h" +#endif namespace SourceXtractor { @@ -98,6 +101,10 @@ class DetectionImageConfig : public Euclid::Configuration::Configuration { double saturation, double flux_scale, int interpolation_gap); DetectionImageExtension(std::shared_ptr fits_image_source, double gain, double saturation, double flux_scale, int interpolation_gap); +#ifdef WITH_ASDF + DetectionImageExtension(std::shared_ptr fits_image_source, double gain, + double saturation, double flux_scale, int interpolation_gap); +#endif static DetectionImageExtension create(std::shared_ptr, const UserValues& args); diff --git a/SEImplementation/src/lib/Configuration/DetectionImageConfig.cpp b/SEImplementation/src/lib/Configuration/DetectionImageConfig.cpp index 9fcafabc4..e9b68e3e5 100644 --- a/SEImplementation/src/lib/Configuration/DetectionImageConfig.cpp +++ b/SEImplementation/src/lib/Configuration/DetectionImageConfig.cpp @@ -22,17 +22,17 @@ #include "Configuration/ConfigManager.h" #include +#include using boost::regex; using boost::regex_match; using boost::smatch; +#include "SEFramework/CoordinateSystem/WCS.h" #include "SEFramework/Image/BufferedImage.h" #include "SEFramework/Image/ImageFileReader.h" #include "SEFramework/Image/ProcessedImage.h" #include "SEFramework/FITS/FitsImageSource.h" -#include "SEFramework/CoordinateSystem/WCS.h" - #include "SEImplementation/Configuration/DetectionImageConfig.h" using namespace Euclid::Configuration; @@ -79,8 +79,9 @@ void DetectionImageConfig::initialize(const UserValues& args) { if (args.find(REFERENCE_IMAGE) != args.end()) { DetectionImageExtension extension; - auto reference_image_source = std::make_shared( - args.find(REFERENCE_IMAGE)->second.as(), 0, ImageTile::FloatImage); + auto image_reader = ImageFileReader::create( + args.find(REFERENCE_IMAGE)->second.as()); + auto reference_image_source = image_reader->get(0); extension.m_coordinate_system = std::make_shared(*reference_image_source); m_extensions.emplace_back(std::move(extension)); @@ -109,7 +110,7 @@ DetectionImageConfig::DetectionImageExtension::DetectionImageExtension( std::shared_ptr image_source, double gain, double saturation, double flux_scale, int interpolation_gap) { init(image_source, gain, saturation, flux_scale, interpolation_gap); - m_coordinate_system = std::make_shared(WCS::identity(2)); + m_coordinate_system = std::make_shared(*image_source); rescale(); } diff --git a/SEImplementation/src/lib/Configuration/MeasurementImageConfig.cpp b/SEImplementation/src/lib/Configuration/MeasurementImageConfig.cpp index ad661ba8c..0f7c79519 100644 --- a/SEImplementation/src/lib/Configuration/MeasurementImageConfig.cpp +++ b/SEImplementation/src/lib/Configuration/MeasurementImageConfig.cpp @@ -31,7 +31,6 @@ #include -#include #include #include #include @@ -196,14 +195,7 @@ void MeasurementImageConfig::initialize(const UserValues&) { } info.m_measurement_image = createMeasurementImage(image_source, py_image.flux_scale); - - // TODO: Instantiating a WCS from ASDF files is not yet supported, use the dummy - // WCS for now. - if (auto fits_image_source = std::dynamic_pointer_cast(image_source)) { - info.m_coordinate_system = std::make_shared(*fits_image_source); - } else { - info.m_coordinate_system = std::make_shared(WCS::identity(2)); - } + info.m_coordinate_system = std::make_shared(*image_source); info.m_gain = py_image.gain / flux_scale; info.m_saturation_level = py_image.saturation * flux_scale; From 52cda33832a93fb06e06c5c26cdf64c0f5a8e01f Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Tue, 28 Oct 2025 15:14:46 +0100 Subject: [PATCH 21/29] Finish implementing reading of WCS data from the ASDF file (where supported) and initializing the wcsprm, etc. Adds some tests based on the once for FITS but as noted it's a different WCS with different results. I tested the expected against wcslib and GWCS as well. --- .../SEFramework/ASDF/AsdfImageSource.h | 3 + SEFramework/auxdir/wcs_header.asdf | Bin 2773 -> 2777 bytes SEFramework/src/lib/ASDF/AsdfImageSource.cpp | 7 ++ SEFramework/src/lib/CoordinateSystem/WCS.cpp | 108 ++++++++++++++---- .../tests/src/CoordinateSystem/WCS_test.cpp | 91 ++++++++++++++- 5 files changed, 181 insertions(+), 28 deletions(-) diff --git a/SEFramework/SEFramework/ASDF/AsdfImageSource.h b/SEFramework/SEFramework/ASDF/AsdfImageSource.h index a482ef125..3d0e92ae9 100644 --- a/SEFramework/SEFramework/ASDF/AsdfImageSource.h +++ b/SEFramework/SEFramework/ASDF/AsdfImageSource.h @@ -130,6 +130,9 @@ class AsdfImageSource : public ImageSource, public std::enable_shared_from_this< void saveTile(ImageTile& tile) override; + /** Get optional FITS WCS */ + std::unique_ptr getFitsWCS() const; + private: AsdfImageSource(const std::string& filename, int image_index, std::optional image_path, diff --git a/SEFramework/auxdir/wcs_header.asdf b/SEFramework/auxdir/wcs_header.asdf index 8a122faccc5956b77c2d2643f2a5975a9bc6a657..05bbdff400c24f39085c1a472ff7dae42a5792df 100644 GIT binary patch delta 84 zcmcaAdQ)`6W+qNU9R AsdfImageSource::getImageTile(int x, int y, int width void AsdfImageSource::saveTile(ImageTile& /* tile */) { throw Elements::Exception() << "Saving to ASDF not yet supported"; }; + + +std::unique_ptr AsdfImageSource::getFitsWCS() const { + auto acc = m_handler->getAccessor(); + auto& file = acc->m_fd; + return file.getFitsWCS(); +} } diff --git a/SEFramework/src/lib/CoordinateSystem/WCS.cpp b/SEFramework/src/lib/CoordinateSystem/WCS.cpp index b07298655..f59cb3430 100644 --- a/SEFramework/src/lib/CoordinateSystem/WCS.cpp +++ b/SEFramework/src/lib/CoordinateSystem/WCS.cpp @@ -47,6 +47,30 @@ static auto logger = Elements::Logging::getLogger("WCS"); decltype(&wcssub) safe_wcssub = &wcssub; + +/** + * wcslib < 5.18 is not fully safe thread, as some functions (like discpy, called by lincpy) + * rely on global variables for determining the allocation sizes. For those versions, this is called + * instead, wrapping the call with a mutex. + */ +static int wrapped_wcssub(int alloc, const struct wcsprm* wcssrc, int* nsub, int axes[], struct wcsprm* wcsdst) { + static std::mutex cpy_mutex; + std::lock_guard lock(cpy_mutex); + + return wcssub(alloc, wcssrc, nsub, axes, wcsdst); +} + + +static void installSafeWcssub(void) { + int wcsver[3]; + wcslib_version(wcsver); + if (wcsver[0] < 5 || (wcsver[0] == 5 && wcsver[1] < 18)) { + logger.info() << "wcslib " << wcsver[0] << "." << wcsver[1] + << " is not fully thread safe, using wrapped lincpy call!"; + safe_wcssub = &wrapped_wcssub; + } +} + /** * Translate the return code from wcspih to an elements exception */ @@ -145,18 +169,6 @@ static void wcsReportWarnings(const char *err_buffer) { } } -/** - * wcslib < 5.18 is not fully safe thread, as some functions (like discpy, called by lincpy) - * rely on global variables for determining the allocation sizes. For those versions, this is called - * instead, wrapping the call with a mutex. - */ -static int wrapped_wcssub(int alloc, const struct wcsprm* wcssrc, int* nsub, int axes[], struct wcsprm* wcsdst) { - static std::mutex cpy_mutex; - std::lock_guard lock(cpy_mutex); - - return wcssub(alloc, wcssrc, nsub, axes, wcsdst); -} - WCS::WCS(const FitsImageSource& fits_image_source) : m_wcs(nullptr, nullptr) { int number_of_records = 0; auto fits_headers = fits_image_source.getFitsHeaders(number_of_records); @@ -215,17 +227,11 @@ void WCS::initFits(char* headers, int number_of_records) { wcsvfree(&nwcs_copy, &ptr); }); - int wcsver[3]; - wcslib_version(wcsver); - if (wcsver[0] < 5 || (wcsver[0] == 5 && wcsver[1] < 18)) { - logger.info() << "wcslib " << wcsver[0] << "." << wcsver[1] - << " is not fully thread safe, using wrapped lincpy call!"; - safe_wcssub = &wrapped_wcssub; - } + installSafeWcssub(); } -#ifdef HAVE_ASDF +#ifdef WITH_ASDF /** WCS initializer from an ASDF file * * Currently this makes a brash assumption: if there is any compatible GWCS @@ -235,9 +241,67 @@ void WCS::initFits(char* headers, int number_of_records) { * Later we will figure out how to work in some config option(s) to explicitly * provide a path to the correct WCS to use if there is any ambiguity. */ -WCS::WCS(const AsdfImageSource& fits_image_source) : m_wcs(nullptr, nullptr) { +WCS::WCS(const AsdfImageSource& asdf_image_source) : m_wcs(nullptr, nullptr) { + // First get whether we even have a WCS in the image + auto fits_wcs = asdf_image_source.getFitsWCS(); + + if (!fits_wcs) { + auto tmp = WCS::identity(2); + m_wcs = std::move(tmp.m_wcs); + return; + } + + wcserr_enable(1); + + int nwcs = 0; + wcsprm *wcs; + wcs = (wcsprm *)malloc(sizeof(*wcs)); + + if (!wcs) { + throw Elements::Exception() << "failed to allocate memory for wcslib"; + } + + // Write warnings to a buffer + wcsprintf_set(nullptr); + + // Initialize the wcsprm with memory allocated for 2 dimensions + wcs->flag = -1; + wcsini(1, 2, wcs); + + // Populate from the AsdfFile::FitsWCS + for (int idx = 0; idx < 2; idx++) { + // WARNING: The GWCS fitwcs_imaging schema (and by extension the libasdf + // GWCS extension) use 0-indexed values for crpix: + // https://github.com/asdf-format/asdf-wcs-schemas/blob/main/resources/schemas/stsci.edu/gwcs/fitswcs_imaging-1.0.0.yaml + wcs->crpix[idx] = fits_wcs->crpix()[idx] + 1.0; + wcs->crval[idx] = fits_wcs->crval()[idx]; + wcs->cdelt[idx] = fits_wcs->cdelt()[idx]; + + const auto ctype = fits_wcs->ctype()[idx]; + if (!ctype.empty()) { + std::strncpy(wcs->ctype[idx], ctype.data(), 9); + } else { + wcs->ctype[idx][0] = '\0'; + } + + for (int jdx = 0; jdx < 2; jdx++) { + wcs->pc[idx * wcs->naxis + jdx] = fits_wcs->pc()[idx][jdx]; + } + } + + int ret = wcsset(wcs); + wcsRaiseOnParseError(ret); + wcsReportWarnings(wcsprintf_buf()); + + m_wcs = decltype(m_wcs)(wcs, [nwcs](wcsprm* ptr) { + int nwcs_copy = nwcs; + wcsfree(ptr); + wcsvfree(&nwcs_copy, &ptr); + }); + + installSafeWcssub(); } -#endif +#endif /* WITH_ASDF */ /** diff --git a/SEFramework/tests/src/CoordinateSystem/WCS_test.cpp b/SEFramework/tests/src/CoordinateSystem/WCS_test.cpp index 25b680a66..4ff84ff53 100644 --- a/SEFramework/tests/src/CoordinateSystem/WCS_test.cpp +++ b/SEFramework/tests/src/CoordinateSystem/WCS_test.cpp @@ -20,6 +20,10 @@ #include #include +#ifdef WITH_ASDF +#include "SEFramework/ASDF/AsdfImageSource.h" +#endif + using namespace SourceXtractor; // The wcs_header.fits file contains the headers extracted from @@ -28,12 +32,21 @@ using namespace SourceXtractor; // Other images with other projections do not trigger the error struct WCSFixture { std::string m_fits_path; - std::shared_ptr m_wcs; + std::shared_ptr m_wcs_fits; +#ifdef WITH_ASDF + std::string m_asdf_path; + std::shared_ptr m_wcs_asdf; +#endif WCSFixture() { m_fits_path = Elements::getAuxiliaryPath("wcs_header.fits").native(); FitsImageSource fits_source(m_fits_path, 0); - m_wcs = std::make_shared(fits_source); + m_wcs_fits = std::make_shared(fits_source); +#ifdef WITH_ASDF + m_asdf_path = Elements::getAuxiliaryPath("wcs_header.asdf").native(); + AsdfImageSource asdf_source(m_asdf_path); + m_wcs_asdf = std::make_shared(asdf_source); +#endif } }; @@ -50,7 +63,7 @@ BOOST_FIXTURE_TEST_CASE(ImageToWorld_test, WCSFixture) { for (size_t i = 0; i < img_coords.size(); ++i) { auto img = img_coords[i]; auto true_world = world_coords[i]; - auto world = m_wcs->imageToWorld(img); + auto world = m_wcs_fits->imageToWorld(img); BOOST_CHECK_CLOSE(world.m_alpha, true_world.m_alpha, 1e-4); BOOST_CHECK_CLOSE(world.m_delta, true_world.m_delta, 1e-4); } @@ -65,7 +78,7 @@ BOOST_FIXTURE_TEST_CASE(WorldToImage_test, WCSFixture) { for (size_t i = 0; i < img_coords.size(); ++i) { auto true_img = img_coords[i]; auto world = world_coords[i]; - auto img = m_wcs->worldToImage(world); + auto img = m_wcs_fits->worldToImage(world); BOOST_CHECK_CLOSE(img.m_x, true_img.m_x, 1e-4); BOOST_CHECK_CLOSE(img.m_y, true_img.m_y, 1e-4); } @@ -74,11 +87,11 @@ BOOST_FIXTURE_TEST_CASE(WorldToImage_test, WCSFixture) { //----------------------------------------------------------------------------- BOOST_FIXTURE_TEST_CASE(ImageOutOfBounds_test, WCSFixture) { - auto world = m_wcs->imageToWorld(ImageCoordinate(-10, -5)); + auto world = m_wcs_fits->imageToWorld(ImageCoordinate(-10, -5)); BOOST_CHECK_CLOSE(world.m_alpha, 231.36376564, 1e-4); BOOST_CHECK_CLOSE(world.m_delta, 30.74723277, 1e-4); - world = m_wcs->imageToWorld(ImageCoordinate(2100, 2100)); + world = m_wcs_fits->imageToWorld(ImageCoordinate(2100, 2100)); BOOST_CHECK_CLOSE(world.m_alpha, 231.62793621, 1e-4); BOOST_CHECK_CLOSE(world.m_delta, 30.8452462, 1e-4); } @@ -94,6 +107,72 @@ BOOST_FIXTURE_TEST_CASE(ImageOutOfBounds_test, WCSFixture) { // //----------------------------------------------------------------------------- + +#ifdef WITH_ASDF +//----------------------------------------------------------------------------- +// tests for WCS from ASDF; this mirrors the FITS tests somewhat but the +// expected results are not the same, as the existing format for embedding a +// FITS WCS in ASDF GWCS (currently used primarily by JWST and RST L3 +// products) is quite a bit simplified and does not include SIP polynomials, +// etc. +// +// Expected values were checked against wcslib manually. +BOOST_FIXTURE_TEST_CASE(ASDF_ImageToWorld_test, WCSFixture) { + std::vector img_coords{{0, 0}, {10, 8}, {55.5, 980.5}}; + std::vector world_coords{ + {269.5464167447208, 65.95659889599523}, + {269.5467894638456, 65.95672214987262}, + {269.548235426113, 65.97157594650089} + }; + + for (size_t i = 0; i < img_coords.size(); ++i) { + auto img = img_coords[i]; + auto true_world = world_coords[i]; + auto world = m_wcs_asdf->imageToWorld(img); + BOOST_CHECK_CLOSE(world.m_alpha, true_world.m_alpha, 1e-4); + BOOST_CHECK_CLOSE(world.m_delta, true_world.m_delta, 1e-4); + } +} + +//----------------------------------------------------------------------------- + +BOOST_FIXTURE_TEST_CASE(ASDF_WorldToImage_test, WCSFixture) { + std::vector img_coords{{0., 0.}, {10., 8.}, {55.5, 980.5}}; + std::vector world_coords{ + {269.5464167447208, 65.95659889599523}, + {269.5467894638456, 65.95672214987262}, + {269.548235426113, 65.97157594650089} + }; + + for (size_t i = 0; i < img_coords.size(); ++i) { + auto true_img = img_coords[i]; + auto world = world_coords[i]; + auto img = m_wcs_asdf->worldToImage(world); + if (i == 0) { + // Value should be close to zero in the first case; BOOST_CHECK_CLOSE + // fails on this + BOOST_CHECK_SMALL(img.m_x, 1e-8); + BOOST_CHECK_SMALL(img.m_y, 1e-8); + } else { + BOOST_CHECK_CLOSE(img.m_x, true_img.m_x, 1e-4); + BOOST_CHECK_CLOSE(img.m_y, true_img.m_y, 1e-4); + } + } +} + +//----------------------------------------------------------------------------- + +BOOST_FIXTURE_TEST_CASE(ASDF_ImageOutOfBounds_test, WCSFixture) { + auto world = m_wcs_asdf->imageToWorld(ImageCoordinate(-10, -5)); + BOOST_CHECK_CLOSE(world.m_alpha, 269.54604322, 1e-4); + BOOST_CHECK_CLOSE(world.m_delta, 65.95652145, 1e-4); + + world = m_wcs_asdf->imageToWorld(ImageCoordinate(2100, 2100)); + BOOST_CHECK_CLOSE(world.m_alpha, 269.62467272, 1e-4); + BOOST_CHECK_CLOSE(world.m_delta, 65.98887495, 1e-4); +} +#endif /* WITH_ASDF */ + BOOST_AUTO_TEST_SUITE_END() //----------------------------------------------------------------------------- From 8f6ee10d0c0543bb045b1f0cf5eb39f1fbbe3f68 Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Tue, 28 Oct 2025 15:36:11 +0100 Subject: [PATCH 22/29] fix: clean up this bit of cargo cult; in this case nwcs = 1 always --- SEFramework/src/lib/CoordinateSystem/WCS.cpp | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/SEFramework/src/lib/CoordinateSystem/WCS.cpp b/SEFramework/src/lib/CoordinateSystem/WCS.cpp index f59cb3430..040be2039 100644 --- a/SEFramework/src/lib/CoordinateSystem/WCS.cpp +++ b/SEFramework/src/lib/CoordinateSystem/WCS.cpp @@ -253,7 +253,6 @@ WCS::WCS(const AsdfImageSource& asdf_image_source) : m_wcs(nullptr, nullptr) { wcserr_enable(1); - int nwcs = 0; wcsprm *wcs; wcs = (wcsprm *)malloc(sizeof(*wcs)); @@ -293,10 +292,8 @@ WCS::WCS(const AsdfImageSource& asdf_image_source) : m_wcs(nullptr, nullptr) { wcsRaiseOnParseError(ret); wcsReportWarnings(wcsprintf_buf()); - m_wcs = decltype(m_wcs)(wcs, [nwcs](wcsprm* ptr) { - int nwcs_copy = nwcs; + m_wcs = decltype(m_wcs)(wcs, [](wcsprm* ptr) { wcsfree(ptr); - wcsvfree(&nwcs_copy, &ptr); }); installSafeWcssub(); From 01bb894902cd379df7d2107cf110bcf9ebb0794d Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Tue, 28 Oct 2025 16:01:26 +0100 Subject: [PATCH 23/29] add ability to get the path of an ndarray (and later an AsdfImageSource) this will be useful when matching WCS's to images --- SEFramework/SEFramework/ASDF/AsdfFile.h | 4 ++++ SEFramework/src/lib/ASDF/AsdfFile.cpp | 3 ++- SEFramework/tests/src/ASDF/AsdfFile_test.cpp | 1 + 3 files changed, 7 insertions(+), 1 deletion(-) diff --git a/SEFramework/SEFramework/ASDF/AsdfFile.h b/SEFramework/SEFramework/ASDF/AsdfFile.h index df8651993..928013b53 100644 --- a/SEFramework/SEFramework/ASDF/AsdfFile.h +++ b/SEFramework/SEFramework/ASDF/AsdfFile.h @@ -109,6 +109,8 @@ class AsdfFile { */ void fillImageTile(const std::shared_ptr image_tile, int layer); + const std::string& getPath() const { return m_path; } + private: /** * Private constructor for creating the `Ndarray` wrapper from a raw asdf_value_t * @@ -121,6 +123,8 @@ class AsdfFile { : m_ndarray_ptr(ptr) {} asdf_ndarray_t* m_ndarray_ptr; + // Store the JSON Path to the ndarray + std::string m_path; ImageTile::ImageType m_image_type; mutable std::unique_ptr> m_shape; }; diff --git a/SEFramework/src/lib/ASDF/AsdfFile.cpp b/SEFramework/src/lib/ASDF/AsdfFile.cpp index 3c70870f6..f2f89d17c 100644 --- a/SEFramework/src/lib/ASDF/AsdfFile.cpp +++ b/SEFramework/src/lib/ASDF/AsdfFile.cpp @@ -237,12 +237,12 @@ AsdfFile::Ndarray::Ndarray(const AsdfFile& file, asdf_value_t* value) : Ndarray((asdf_ndarray_t*)nullptr) { asdf_ndarray_t* ndarray_ptr = nullptr; asdf_value_err_t err = asdf_value_as_ndarray(value, &ndarray_ptr); + const char* path = asdf_value_path(value); switch (err) { case ASDF_VALUE_OK: // Value exists and is an ndarray: OK break; case ASDF_VALUE_ERR_TYPE_MISMATCH: { - const char* path = asdf_value_path(value); throw AsdfValueTypeMismatchException() << "Value at " << path << " is not an ndarray: " << file.m_path.native(); } @@ -253,6 +253,7 @@ AsdfFile::Ndarray::Ndarray(const AsdfFile& file, asdf_value_t* value) } } m_ndarray_ptr = ndarray_ptr; + m_path = path; m_image_type = convertDatatypeToImageType(&ndarray_ptr->datatype); } diff --git a/SEFramework/tests/src/ASDF/AsdfFile_test.cpp b/SEFramework/tests/src/ASDF/AsdfFile_test.cpp index a9c38f617..b973f930a 100644 --- a/SEFramework/tests/src/ASDF/AsdfFile_test.cpp +++ b/SEFramework/tests/src/ASDF/AsdfFile_test.cpp @@ -87,6 +87,7 @@ BOOST_FIXTURE_TEST_CASE( get_ndarray_by_path, AsdfFileFixture ) { shape.begin(), shape.end(), expected_shape.begin(), expected_shape.end() ); + BOOST_CHECK_EQUAL(ndarray->getPath(), "/PRIMARY"); } From 5b449a6d8173d2e2e3dbc930face9138ca2f6818 Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Fri, 31 Oct 2025 12:07:49 +0100 Subject: [PATCH 24/29] add extra options, when compiled with ASDF support, for setting the WCS path for detection images and reference images this just sets up the wiring--still need to implement parsing and interpretation of the WCS path (and figurout where to document the format) --- .../SEFramework/ASDF/AsdfImageSource.h | 3 + .../SEFramework/CoordinateSystem/WCS.h | 10 ++- SEFramework/src/lib/ASDF/AsdfImageSource.cpp | 7 ++ SEFramework/src/lib/CoordinateSystem/WCS.cpp | 37 ++++++---- SEImplementation/CMakeLists.txt | 7 ++ .../Configuration/DetectionImageConfig.h | 3 +- .../Configuration/DetectionImageConfig.cpp | 69 ++++++++++++++++++- 7 files changed, 118 insertions(+), 18 deletions(-) diff --git a/SEFramework/SEFramework/ASDF/AsdfImageSource.h b/SEFramework/SEFramework/ASDF/AsdfImageSource.h index 3d0e92ae9..f08243d60 100644 --- a/SEFramework/SEFramework/ASDF/AsdfImageSource.h +++ b/SEFramework/SEFramework/ASDF/AsdfImageSource.h @@ -133,6 +133,9 @@ class AsdfImageSource : public ImageSource, public std::enable_shared_from_this< /** Get optional FITS WCS */ std::unique_ptr getFitsWCS() const; + /** Get optional FITS WCS with optional path lookup map */ + std::unique_ptr getFitsWCS(std::optional wcs_path) const; + private: AsdfImageSource(const std::string& filename, int image_index, std::optional image_path, diff --git a/SEFramework/SEFramework/CoordinateSystem/WCS.h b/SEFramework/SEFramework/CoordinateSystem/WCS.h index 352e8d461..ec85e0304 100644 --- a/SEFramework/SEFramework/CoordinateSystem/WCS.h +++ b/SEFramework/SEFramework/CoordinateSystem/WCS.h @@ -31,6 +31,7 @@ #include "SEFramework/FITS/FitsImageSource.h" #include "SEFramework/Image/ImageSource.h" #ifdef WITH_ASDF +#include "SEFramework/ASDF/AsdfFile.h" #include "SEFramework/ASDF/AsdfImageSource.h" #endif @@ -41,11 +42,12 @@ namespace SourceXtractor { class WCS : public CoordinateSystem { public: explicit WCS(const FitsImageSource& fits_image_source); + explicit WCS(const ImageSource& image_source); + explicit WCS(const WCS& original); #ifdef WITH_ASDF explicit WCS(const AsdfImageSource& asdf_image_source); + explicit WCS(const AsdfImageSource& asdf_image_source, std::optional wcs_path); #endif - explicit WCS(const ImageSource& image_source); - explicit WCS(const WCS& original); virtual ~WCS(); @@ -65,6 +67,10 @@ class WCS : public CoordinateSystem { void initFits(char* headers, int number_of_records); +#ifdef WITH_ASDF + void initAsdf(std::unique_ptr fits_wcs); +#endif + std::unique_ptr> m_wcs; }; diff --git a/SEFramework/src/lib/ASDF/AsdfImageSource.cpp b/SEFramework/src/lib/ASDF/AsdfImageSource.cpp index bff746843..71d5450b4 100644 --- a/SEFramework/src/lib/ASDF/AsdfImageSource.cpp +++ b/SEFramework/src/lib/ASDF/AsdfImageSource.cpp @@ -125,4 +125,11 @@ std::unique_ptr AsdfImageSource::getFitsWCS() const { auto& file = acc->m_fd; return file.getFitsWCS(); } + +// TODO: Implement handling of wcs_path_map +std::unique_ptr AsdfImageSource::getFitsWCS(std::optional wcs_path) const { + auto acc = m_handler->getAccessor(); + auto& file = acc->m_fd; + return file.getFitsWCS(); +} } diff --git a/SEFramework/src/lib/CoordinateSystem/WCS.cpp b/SEFramework/src/lib/CoordinateSystem/WCS.cpp index 040be2039..4f4027b93 100644 --- a/SEFramework/src/lib/CoordinateSystem/WCS.cpp +++ b/SEFramework/src/lib/CoordinateSystem/WCS.cpp @@ -232,19 +232,7 @@ void WCS::initFits(char* headers, int number_of_records) { #ifdef WITH_ASDF -/** WCS initializer from an ASDF file - * - * Currently this makes a brash assumption: if there is any compatible GWCS - * object in the file it "must" be the right one. This assumption can be wrong - * but in practice most ASDF files have one data array, one WCS. - * - * Later we will figure out how to work in some config option(s) to explicitly - * provide a path to the correct WCS to use if there is any ambiguity. - */ -WCS::WCS(const AsdfImageSource& asdf_image_source) : m_wcs(nullptr, nullptr) { - // First get whether we even have a WCS in the image - auto fits_wcs = asdf_image_source.getFitsWCS(); - +void WCS::initAsdf(std::unique_ptr fits_wcs) { if (!fits_wcs) { auto tmp = WCS::identity(2); m_wcs = std::move(tmp.m_wcs); @@ -298,6 +286,29 @@ WCS::WCS(const AsdfImageSource& asdf_image_source) : m_wcs(nullptr, nullptr) { installSafeWcssub(); } + + +/** WCS initializer from an ASDF file + * + * Currently this makes a brash assumption: if there is any compatible GWCS + * object in the file it "must" be the right one. This assumption can be wrong + * but in practice most ASDF files have one data array, one WCS. + * + * Later we will figure out how to work in some config option(s) to explicitly + * provide a path to the correct WCS to use if there is any ambiguity. + */ +WCS::WCS(const AsdfImageSource& asdf_image_source) : m_wcs(nullptr, nullptr) { + // First get whether we even have a WCS in the image + auto fits_wcs = asdf_image_source.getFitsWCS(); + initAsdf(std::move(fits_wcs)); +} + + +WCS::WCS(const AsdfImageSource& asdf_image_source, std::optional wcs_path) : m_wcs(nullptr, nullptr) { + // First get whether we even have a WCS in the image + auto fits_wcs = asdf_image_source.getFitsWCS(wcs_path); + initAsdf(std::move(fits_wcs)); +} #endif /* WITH_ASDF */ diff --git a/SEImplementation/CMakeLists.txt b/SEImplementation/CMakeLists.txt index 1b52a93d3..100b2f612 100644 --- a/SEImplementation/CMakeLists.txt +++ b/SEImplementation/CMakeLists.txt @@ -37,6 +37,7 @@ find_package(BoostPython ${PYTHON_EXPLICIT_VERSION}) find_package(Boost 1.53 REQUIRED) find_package(OnnxRuntime) find_package(Log4CPP REQUIRED) +find_package(ASDF) if (${Boost_VERSION} LESS "106700") message(WARNING " @@ -48,6 +49,12 @@ if (${PYTHON_VERSION_MAJOR} VERSION_LESS 3) message(DEPRECATION "Python 2 support is deprecated and will be removed (found ${PYTHON_VERSION_STRING})") endif () +if(NOT ASDF_FOUND) + message("libasdf not found, compiling without ASDF support") +else() + add_definitions(-DWITH_ASDF) +endif() + #=============================================================================== # Declare the library dependencies here # Example: diff --git a/SEImplementation/SEImplementation/Configuration/DetectionImageConfig.h b/SEImplementation/SEImplementation/Configuration/DetectionImageConfig.h index 7d2a5270b..92896a120 100644 --- a/SEImplementation/SEImplementation/Configuration/DetectionImageConfig.h +++ b/SEImplementation/SEImplementation/Configuration/DetectionImageConfig.h @@ -103,7 +103,8 @@ class DetectionImageConfig : public Euclid::Configuration::Configuration { double saturation, double flux_scale, int interpolation_gap); #ifdef WITH_ASDF DetectionImageExtension(std::shared_ptr fits_image_source, double gain, - double saturation, double flux_scale, int interpolation_gap); + double saturation, double flux_scale, int interpolation_gap, + std::optional wcs_path); #endif static DetectionImageExtension create(std::shared_ptr, const UserValues& args); diff --git a/SEImplementation/src/lib/Configuration/DetectionImageConfig.cpp b/SEImplementation/src/lib/Configuration/DetectionImageConfig.cpp index e9b68e3e5..2d47ba100 100644 --- a/SEImplementation/src/lib/Configuration/DetectionImageConfig.cpp +++ b/SEImplementation/src/lib/Configuration/DetectionImageConfig.cpp @@ -33,6 +33,10 @@ using boost::smatch; #include "SEFramework/Image/ProcessedImage.h" #include "SEFramework/FITS/FitsImageSource.h" +#ifdef WITH_ASDF +#include "SEFramework/ASDF/AsdfImageSource.h" +#endif + #include "SEImplementation/Configuration/DetectionImageConfig.h" using namespace Euclid::Configuration; @@ -48,15 +52,32 @@ static const std::string DETECTION_IMAGE_SATURATION { "detection-image-saturatio static const std::string DETECTION_IMAGE_INTERPOLATION { "detection-image-interpolation" }; static const std::string DETECTION_IMAGE_INTERPOLATION_GAP { "detection-image-interpolation-gap" }; +#ifdef WITH_ASDF +static const std::string DETECTION_IMAGE_ASDF_WCS_PATH { "detection-image-asdf-wcs-path" }; +static const std::string REFERENCE_IMAGE_ASDF_WCS_PATH { "reference-image-asdf-wcs-path" }; +#endif + DetectionImageConfig::DetectionImageConfig(long manager_id) : Configuration(manager_id) {} std::map DetectionImageConfig::getProgramOptions() { return { {"Detection image", { {DETECTION_IMAGE.c_str(), po::value(), - "Path to a fits format image to be used as detection image."}, + // NOTE: In principle ImageFileReader can be used to programmatically + // build a string listing the supported file types, but as currently + // there are only two and one of them only with optional support we + // can use an ifdef for now. +#ifdef WITH_ASDF + "Path to a FITS or ASDF format image to be used as detection image."}, +#else + "Path to a FITS format image to be used as detection image."}, +#endif {REFERENCE_IMAGE.c_str(), po::value(), - "Path to a fits format image to be used as coordinates reference only."}, +#ifdef WITH_ASDF + "Path to a FITS or ASDF format image to be used as coordinates reference only."}, +#else + "Path to a FITS or ASDF format image to be used as coordinates reference only."}, +#endif {DETECTION_IMAGE_GAIN.c_str(), po::value(), "Detection image gain in e-/ADU (0 = infinite gain)"}, {DETECTION_IMAGE_FLUX_SCALE.c_str(), po::value(), @@ -67,6 +88,13 @@ std::map DetectionImageConfig "Interpolate bad pixels in detection image"}, {DETECTION_IMAGE_INTERPOLATION_GAP.c_str(), po::value()->default_value(5), "Maximum number if pixels to interpolate over"} +#ifdef WITH_ASDF + , + {DETECTION_IMAGE_ASDF_WCS_PATH.c_str(), po::value(), + "When reading from an ASDF file, the JSON path to the WCS to use for the detection image"}, + {REFERENCE_IMAGE_ASDF_WCS_PATH.c_str(), po::value(), + "When reading from an ASDF file, the JSON path to the WCS to use for the reference image"} +#endif }}}; } @@ -82,7 +110,18 @@ void DetectionImageConfig::initialize(const UserValues& args) { auto image_reader = ImageFileReader::create( args.find(REFERENCE_IMAGE)->second.as()); auto reference_image_source = image_reader->get(0); +#ifndef WITH_ASDF extension.m_coordinate_system = std::make_shared(*reference_image_source); +#else + if (auto asdf_src = std::dynamic_pointer_cast(reference_image_source)) { + std::optional wcs_path; + if (auto it = args.find(REFERENCE_IMAGE_ASDF_WCS_PATH); it != args.end()) + wcs_path = it->second.as(); + extension.m_coordinate_system = std::make_shared(*asdf_src, wcs_path); + } else { + extension.m_coordinate_system = std::make_shared(*reference_image_source); + } +#endif m_extensions.emplace_back(std::move(extension)); m_is_reference_image = true; @@ -162,6 +201,23 @@ DetectionImageConfig::DetectionImageExtension::DetectionImageExtension( } +#ifdef WITH_ASDF +/** + * Special case for loading from ASDF when enabled, to use the detection-image-asdf-wcs-path + * option if given + */ +DetectionImageConfig::DetectionImageExtension::DetectionImageExtension( + std::shared_ptr asdf_image_source, double gain, double saturation, + double flux_scale, int interpolation_gap, std::optional wcs_path) { + init(asdf_image_source, gain, saturation, flux_scale, interpolation_gap); + m_coordinate_system = std::make_shared(*asdf_image_source, wcs_path); + rescale(); +} + + +#endif /* WITH_ASDF */ + + void DetectionImageConfig::DetectionImageExtension::init( std::shared_ptr image_source, double gain, double saturation, double flux_scale, int interpolation_gap) { @@ -195,6 +251,15 @@ DetectionImageConfig::DetectionImageExtension DetectionImageConfig::DetectionIma if (auto fits_src = std::dynamic_pointer_cast(image_source)) { return DetectionImageExtension(fits_src, gain, saturation, flux_scale, interpolation_gap); } +#ifdef WITH_ASDF + else if (auto asdf_src = std::dynamic_pointer_cast(image_source)) { + std::optional wcs_path; + if (auto it = args.find(DETECTION_IMAGE_ASDF_WCS_PATH); it != args.end()) + wcs_path = it->second.as(); + return DetectionImageExtension(asdf_src, gain, saturation, flux_scale, interpolation_gap, + wcs_path); + } +#endif return DetectionImageExtension(image_source, gain, saturation, flux_scale, interpolation_gap); } From cecb22f43db62b0fe96bbe99fb628acf5a3f05ab Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Fri, 31 Oct 2025 13:23:36 +0100 Subject: [PATCH 25/29] Implement parsing of ASDF WCS path mappings and add tests --- .../SEFramework/ASDF/AsdfImageSource.h | 4 + SEFramework/auxdir/wcs_header.asdf | Bin 2777 -> 6044 bytes SEFramework/src/lib/ASDF/AsdfImageSource.cpp | 69 +++++++++++++++++- .../tests/src/ASDF/AsdfImageSource_test.cpp | 53 +++++++++++++- 4 files changed, 124 insertions(+), 2 deletions(-) diff --git a/SEFramework/SEFramework/ASDF/AsdfImageSource.h b/SEFramework/SEFramework/ASDF/AsdfImageSource.h index f08243d60..88725f9f9 100644 --- a/SEFramework/SEFramework/ASDF/AsdfImageSource.h +++ b/SEFramework/SEFramework/ASDF/AsdfImageSource.h @@ -142,6 +142,10 @@ class AsdfImageSource : public ImageSource, public std::enable_shared_from_this< ImageTile::ImageType image_type, std::shared_ptr manager); + using WcsPathMap = std::unordered_map; + + WcsPathMap parseWcsPath(const std::string& wcs_path) const; + std::string m_filename; ImageTile::ImageType m_image_type; std::shared_ptr m_file_manager; diff --git a/SEFramework/auxdir/wcs_header.asdf b/SEFramework/auxdir/wcs_header.asdf index 05bbdff400c24f39085c1a472ff7dae42a5792df..9f5bad3b89fa37965e74f890223a8d1310b14adc 100644 GIT binary patch delta 679 zcmaKoze~eF6vw%1o01G-5m6~ftXM?b%cbUGqk=f7xHvcnP7>Q_Gc=JTR*MK-#UF!j z)YZjFEclP;;@{%p`REv+>z0W=&Xeu*6)C?87EtExo?)@5oDUYT10a4KM}#nvqepTbkF1qqG6Mdh^S7s6jxF;GB!2h a(p4}rF@jK*V9L}COqrQjaOvskaRC5y6e24C diff --git a/SEFramework/src/lib/ASDF/AsdfImageSource.cpp b/SEFramework/src/lib/ASDF/AsdfImageSource.cpp index 71d5450b4..d36afde13 100644 --- a/SEFramework/src/lib/ASDF/AsdfImageSource.cpp +++ b/SEFramework/src/lib/ASDF/AsdfImageSource.cpp @@ -21,14 +21,21 @@ * Created on: Sep 03, 2025 * Author: embray */ +#include +#include #include +#include "ElementsKernel/Logging.h" + #include "SEFramework/ASDF/AsdfFile.h" #include "SEFramework/ASDF/AsdfImageSource.h" namespace SourceXtractor { + +static auto logger = Elements::Logging::getLogger("ASDF"); + AsdfImageSource::AsdfImageSource(const std::string& filename, int image_index, std::optional image_path, ImageTile::ImageType image_type, @@ -120,6 +127,51 @@ void AsdfImageSource::saveTile(ImageTile& /* tile */) { }; +/** + * Prefix an ASDF path with '/' if it doesn't already have, to ensure all paths are + * normalized; also trims whitespace + */ +inline std::string asdfNormalizePath(std::string s) { + s = boost::trim_copy(s); + if (!s.empty() && s.front() != '/') { + s.insert(s.begin(), '/'); + } + return s; +} + + +AsdfImageSource::WcsPathMap AsdfImageSource::parseWcsPath(const std::string& wcs_path) const { + WcsPathMap path_map; + + size_t nsep = std::count(wcs_path.begin(), wcs_path.end(), ':'); + + // Simple case, a single WCS which we indicate with '*' + if (nsep == 0) { + path_map["*"] = asdfNormalizePath(wcs_path); + return path_map; + } + + std::vector tokens; + boost::split(tokens, wcs_path, boost::is_any_of(":")); + for (auto& token : tokens) { + boost::trim(token); + } + + if (tokens.size() % 2 != 0) { + throw Elements::Exception() << "Invalid WCS path given: " << wcs_path << "; odd number of " + "path elements"; + } + + for (size_t idx = 0; idx < tokens.size(); idx += 2) { + std::string key = asdfNormalizePath(tokens[idx]); + std::string value = asdfNormalizePath(tokens[idx + 1]); + path_map[key] = value; + } + + return path_map; +} + + std::unique_ptr AsdfImageSource::getFitsWCS() const { auto acc = m_handler->getAccessor(); auto& file = acc->m_fd; @@ -130,6 +182,21 @@ std::unique_ptr AsdfImageSource::getFitsWCS() const { std::unique_ptr AsdfImageSource::getFitsWCS(std::optional wcs_path) const { auto acc = m_handler->getAccessor(); auto& file = acc->m_fd; - return file.getFitsWCS(); + + if (wcs_path == std::nullopt) { + return file.getFitsWCS(); + } + + WcsPathMap path_map = parseWcsPath(*wcs_path); + const std::string& ndarray_path = m_ndarray->getPath(); + + if (path_map.find(ndarray_path) != path_map.end()) { + return file.getFitsWCS(path_map.find(ndarray_path)->second); + } else if (path_map.find("*") != path_map.end()) { + return file.getFitsWCS(path_map.find("*")->second); + } else { + logger.warn() << "No WCS path found for ndarray at " << ndarray_path << "; trying any WCS"; + return file.getFitsWCS(); + } } } diff --git a/SEFramework/tests/src/ASDF/AsdfImageSource_test.cpp b/SEFramework/tests/src/ASDF/AsdfImageSource_test.cpp index 140dd45fc..e63fe0abb 100644 --- a/SEFramework/tests/src/ASDF/AsdfImageSource_test.cpp +++ b/SEFramework/tests/src/ASDF/AsdfImageSource_test.cpp @@ -30,11 +30,12 @@ using namespace SourceXtractor; struct AsdfImageSourceFixture { - std::string mhdu_path, primary_path; + std::string mhdu_path, primary_path, wcs_header; AsdfImageSourceFixture() { mhdu_path = Elements::getAuxiliaryPath("multiple_hdu.asdf").native(); primary_path = Elements::getAuxiliaryPath("with_primary.asdf").native(); + wcs_header = Elements::getAuxiliaryPath("wcs_header.asdf").native(); } }; @@ -100,6 +101,56 @@ BOOST_FIXTURE_TEST_CASE(empty_ndarray_test, AsdfImageSourceFixture) { BOOST_CHECK_THROW(std::make_shared(mhdu_path, "PRIMARY"), Elements::Exception); } + +BOOST_FIXTURE_TEST_CASE(get_fitswcs_auto, AsdfImageSourceFixture) { + auto img_src = std::make_shared(wcs_header, "data", ImageTile::FloatImage); + auto fits_wcs = img_src->getFitsWCS(); + BOOST_CHECK(fits_wcs != nullptr); + // The full WCS is checked in AsdfFile_test, but here make sure we just grab the expected + // one. There are two WCS in the file each with different expected ctypes. The first one, + // "wcs", should be TAN, the other is AIR + auto ctype = fits_wcs->ctype(); + std::array expected_ctype{{"RA---TAN", "DEC--TAN"}}; + BOOST_CHECK_EQUAL_COLLECTIONS( + ctype.begin(), ctype.end(), + expected_ctype.begin(), expected_ctype.end() + ); +} + + +BOOST_FIXTURE_TEST_CASE(get_fitswcs_simple_path, AsdfImageSourceFixture) { + auto img_src = std::make_shared(wcs_header, "data", ImageTile::FloatImage); + auto fits_wcs = img_src->getFitsWCS("wcs2"); + BOOST_CHECK(fits_wcs != nullptr); + // The full WCS is checked in AsdfFile_test, but here make sure we just grab the expected + // one. There are two WCS in the file each with different expected ctypes. The first one, + // "wcs", should be TAN, the other is AIR + auto ctype = fits_wcs->ctype(); + std::array expected_ctype{{"RA---AIR", "DEC--AIR"}}; + BOOST_CHECK_EQUAL_COLLECTIONS( + ctype.begin(), ctype.end(), + expected_ctype.begin(), expected_ctype.end() + ); +} + + +BOOST_FIXTURE_TEST_CASE(get_fitswcs_mapped_path, AsdfImageSourceFixture) { + auto img_src = std::make_shared(wcs_header, "data", ImageTile::FloatImage); + auto fits_wcs = img_src->getFitsWCS("data:wcs2"); + BOOST_CHECK(fits_wcs != nullptr); + // The full WCS is checked in AsdfFile_test, but here make sure we just grab the expected + // one. There are two WCS in the file each with different expected ctypes. The first one, + // "wcs", should be TAN, the other is AIR + auto ctype = fits_wcs->ctype(); + std::array expected_ctype{{"RA---AIR", "DEC--AIR"}}; + BOOST_CHECK_EQUAL_COLLECTIONS( + ctype.begin(), ctype.end(), + expected_ctype.begin(), expected_ctype.end() + ); + + BOOST_CHECK_THROW(img_src->getFitsWCS("data:does-not-exist"), Elements::Exception); +} + //----------------------------------------------------------------------------- BOOST_AUTO_TEST_SUITE_END() From e4077f4a58efacdc19b59278a80932c5372e5e7f Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Fri, 31 Oct 2025 14:03:10 +0100 Subject: [PATCH 26/29] build: ensure libasdf was compiled with GWCS support kinda hacky but gets the job done for now; later may change this when the libasdf gwcs module gets moved out to a separate library/plugin --- cmake/modules/FindASDF.cmake | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/cmake/modules/FindASDF.cmake b/cmake/modules/FindASDF.cmake index 28351906a..764339c63 100644 --- a/cmake/modules/FindASDF.cmake +++ b/cmake/modules/FindASDF.cmake @@ -24,9 +24,17 @@ if(NOT ASDF_FOUND) HINTS ENV ASDF_ROOT_DIR ASDF_INSTALL_DIR PATH_SUFFIXES lib) - # You can add other deps here if libasdf needs them at link time - set(ASDF_LIBRARIES ${ASDF_LIBRARY} ${FYAML_LIBRARY}) - set(ASDF_INCLUDE_DIRS ${ASDF_INCLUDE_DIR}) + # Check for gwcs support as well + find_path(ASDF_GWCS_INCLUDE_DIR gwcs.h + HINTS ASDF_INCLUDE_DIR + PATH_SUFFIXES asdf/gwcs) + if (ASDF_GWCS_INCLUDE_DIR STREQUAL "ASDF_GWCS_INCLUDE_DIR-NOTFOUND") + message(WARNING "libasdf found, but no GWCS support available.") + else() + # You can add other deps here if libasdf needs them at link time + set(ASDF_LIBRARIES ${ASDF_LIBRARY} ${FYAML_LIBRARY}) + set(ASDF_INCLUDE_DIRS ${ASDF_INCLUDE_DIR}) + endif() include(FindPackageHandleStandardArgs) find_package_handle_standard_args(ASDF @@ -37,5 +45,4 @@ if(NOT ASDF_FOUND) list(REMOVE_DUPLICATES ASDF_LIBRARIES) list(REMOVE_DUPLICATES ASDF_INCLUDE_DIRS) - endif() From 67ef443c045dcf41d3f28f92052716a5f7d7d370 Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Mon, 3 Nov 2025 15:45:02 +0100 Subject: [PATCH 27/29] fix: a handful of memory leaks in the ASDF code --- SEFramework/SEFramework/ASDF/AsdfFile.h | 10 +++---- SEFramework/src/lib/ASDF/AsdfFile.cpp | 40 ++++++++++++++++--------- 2 files changed, 30 insertions(+), 20 deletions(-) diff --git a/SEFramework/SEFramework/ASDF/AsdfFile.h b/SEFramework/SEFramework/ASDF/AsdfFile.h index 928013b53..b80d61397 100644 --- a/SEFramework/SEFramework/ASDF/AsdfFile.h +++ b/SEFramework/SEFramework/ASDF/AsdfFile.h @@ -138,7 +138,9 @@ class AsdfFile { friend class AsdfFile; ~FitsWCS() { - asdf_gwcs_fits_destroy(m_gwcs_fits_ptr); + asdf_gwcs_destroy(m_gwcs_ptr); + m_gwcs_fits_ptr = nullptr; + m_gwcs_ptr = nullptr; } std::array crpix() const noexcept { @@ -172,12 +174,8 @@ class AsdfFile { * Private constructor for creating the `Ndarray` wrapper from a raw asdf_value_t * */ explicit FitsWCS(const AsdfFile& file, asdf_value_t *ptr); - /** - * Private constructor for creating the `Ndarray` wrapper from a raw asdf_ndarray_t * - */ - explicit FitsWCS(asdf_gwcs_fits_t *ptr) - : m_gwcs_fits_ptr(ptr) {} + asdf_gwcs_t *m_gwcs_ptr; asdf_gwcs_fits_t* m_gwcs_fits_ptr; }; diff --git a/SEFramework/src/lib/ASDF/AsdfFile.cpp b/SEFramework/src/lib/ASDF/AsdfFile.cpp index f2f89d17c..9decbd4e3 100644 --- a/SEFramework/src/lib/ASDF/AsdfFile.cpp +++ b/SEFramework/src/lib/ASDF/AsdfFile.cpp @@ -60,17 +60,17 @@ void AsdfFile::open() { if (ptr == nullptr) { throw Elements::Exception() << "Can't open ASDF file: "; - } else { - // Check if the file was created but has an error condition set - const char *error_message = asdf_error(ptr); - - if (error_message != nullptr) { - throw Elements::Exception() - << "Can't open ASDF file: " << m_path.native() << " reason: " << error_message; - } } m_asdf_ptr.reset(ptr); + + // Check if the file was created but has an error condition set + const char *error_message = asdf_error(ptr); + + if (error_message != nullptr) { + throw Elements::Exception() + << "Can't open ASDF file: " << m_path.native() << " reason: " << error_message; + } } // TODO: Support negative indexing as well @@ -90,8 +90,16 @@ std::unique_ptr AsdfFile::getNdarray(int index) { asdf_value_t* value = asdf_mapping_item_value(item); if (asdf_value_is_ndarray(value)) { if (count == index) { - return std::unique_ptr(new Ndarray(*this, value)); + // Have to clone the value before destroying the mapping item + // This is an unfortunate foot gun that needs to be fixed in the + // asdf_mapping_iter interface... + AsdfValuePtr value_clone = AsdfValuePtr(asdf_value_clone(value)); + asdf_mapping_item_destroy(item); + return std::unique_ptr(new Ndarray(*this, value_clone.get())); } else if (count > index) { + // Breaking the loop early means we have to manually clean up the + // mapping item + asdf_mapping_item_destroy(item); break; } count++; @@ -117,7 +125,9 @@ bool fitsWcsValuePredicate(asdf_value_t* value) { return false; } - return asdf_gwcs_is_fits((asdf_file_t*)asdf_value_file(value), gwcs); + bool is_fits = asdf_gwcs_is_fits((asdf_file_t*)asdf_value_file(value), gwcs); + asdf_gwcs_destroy(gwcs); + return is_fits; } @@ -138,7 +148,8 @@ std::unique_ptr AsdfFile::getFitsWCS() { return std::unique_ptr{}; } - asdf_value_t* value = asdf_find_item_value(item); + asdf_value_t* value = asdf_value_clone(asdf_find_item_value(item)); + asdf_find_item_destroy(item); return std::unique_ptr(new FitsWCS(*this, value)); } @@ -238,6 +249,7 @@ AsdfFile::Ndarray::Ndarray(const AsdfFile& file, asdf_value_t* value) asdf_ndarray_t* ndarray_ptr = nullptr; asdf_value_err_t err = asdf_value_as_ndarray(value, &ndarray_ptr); const char* path = asdf_value_path(value); + m_path = path; switch (err) { case ASDF_VALUE_OK: // Value exists and is an ndarray: OK @@ -252,8 +264,8 @@ AsdfFile::Ndarray::Ndarray(const AsdfFile& file, asdf_value_t* value) << file.m_path.native() << ": " << error_message; } } + m_ndarray_ptr = ndarray_ptr; - m_path = path; m_image_type = convertDatatypeToImageType(&ndarray_ptr->datatype); } @@ -312,8 +324,7 @@ bool AsdfFile::Ndarray::isSupportedImage() const { * * The full GWCS is needed in order to properly read the FITS WCS out of it. */ -AsdfFile::FitsWCS::FitsWCS(const AsdfFile& file, asdf_value_t* value) - : FitsWCS((asdf_gwcs_fits_t*)nullptr) { +AsdfFile::FitsWCS::FitsWCS(const AsdfFile& file, asdf_value_t* value) { asdf_gwcs_t* gwcs_ptr = nullptr; asdf_gwcs_fits_t* gwcs_fits_ptr = nullptr; asdf_value_err_t err = asdf_value_as_gwcs(value, &gwcs_ptr); @@ -345,6 +356,7 @@ AsdfFile::FitsWCS::FitsWCS(const AsdfFile& file, asdf_value_t* value) gwcs_fits_ptr = (asdf_gwcs_fits_t*)step0->transform; assert(gwcs_fits_ptr); + m_gwcs_ptr = gwcs_ptr; m_gwcs_fits_ptr = gwcs_fits_ptr; // We don't need the asdf_value_t anymore at this point and can release it. asdf_value_destroy(value); From 96bb8d943f25657474bbc2a823ea3841c9b33ac1 Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Mon, 3 Nov 2025 17:21:55 +0100 Subject: [PATCH 28/29] fix: improve handling of the wcsprm pointer and fix some minor memory leaks --- .../SEFramework/CoordinateSystem/WCS.h | 38 +++++++++- SEFramework/src/lib/CoordinateSystem/WCS.cpp | 76 ++++++++++--------- 2 files changed, 75 insertions(+), 39 deletions(-) diff --git a/SEFramework/SEFramework/CoordinateSystem/WCS.h b/SEFramework/SEFramework/CoordinateSystem/WCS.h index ec85e0304..389f018b4 100644 --- a/SEFramework/SEFramework/CoordinateSystem/WCS.h +++ b/SEFramework/SEFramework/CoordinateSystem/WCS.h @@ -27,6 +27,9 @@ #include #include +#include +#include + #include "SEFramework/CoordinateSystem/CoordinateSystem.h" #include "SEFramework/FITS/FitsImageSource.h" #include "SEFramework/Image/ImageSource.h" @@ -62,16 +65,43 @@ class WCS : public CoordinateSystem { void addOffset(PixelCoordinate pc); private: - explicit WCS(std::unique_ptr> wcs) - : m_wcs(std::move(wcs)) {} - void initFits(char* headers, int number_of_records); #ifdef WITH_ASDF void initAsdf(std::unique_ptr fits_wcs); #endif - std::unique_ptr> m_wcs; + struct WcsprmDestroy { + int nwcs; + bool owned; + + void operator()(wcsprm* wcs) { + if (!wcs) + return; + + if (nwcs > 0) { + wcsvfree(&nwcs, &wcs); + } else { + wcsfree(wcs); + if (owned) { + delete wcs; + } + } + } + }; + + using WcsprmPtr = std::unique_ptr; + + WcsprmPtr m_wcs; + + static WcsprmPtr make_wcsprm_ptr(); + static WcsprmPtr make_wcsprm_ptr(wcsprm* wcs); + static WcsprmPtr make_wcsprm_ptr(wcsprm* wcs, bool owned); + static WcsprmPtr make_wcsprm_ptr(wcsprm* wcs, bool owned, int nwcs); + + explicit WCS(WcsprmPtr wcs) + : m_wcs(std::move(wcs)) {} + }; } diff --git a/SEFramework/src/lib/CoordinateSystem/WCS.cpp b/SEFramework/src/lib/CoordinateSystem/WCS.cpp index 4f4027b93..07b67afec 100644 --- a/SEFramework/src/lib/CoordinateSystem/WCS.cpp +++ b/SEFramework/src/lib/CoordinateSystem/WCS.cpp @@ -169,14 +169,14 @@ static void wcsReportWarnings(const char *err_buffer) { } } -WCS::WCS(const FitsImageSource& fits_image_source) : m_wcs(nullptr, nullptr) { +WCS::WCS(const FitsImageSource& fits_image_source) { int number_of_records = 0; auto fits_headers = fits_image_source.getFitsHeaders(number_of_records); initFits(&(*fits_headers)[0], number_of_records); } -WCS::WCS(const WCS& original) : m_wcs(nullptr, nullptr) { +WCS::WCS(const WCS& original) { //FIXME Horrible hack: I couldn't figure out how to properly do a deep copy wcsprm so instead // of making a copy, I use the ascii headers output from the original to recreate a new one @@ -212,6 +212,9 @@ void WCS::initFits(char* headers, int number_of_records) { int ret = wcspih(headers, number_of_records, WCSHDR_strict, 2, &nreject_strict, &nwcs, &wcs); wcsRaiseOnParseError(ret); wcsReportWarnings(wcsprintf_buf()); + // It's still necessary to do wcsvfree before the second pass, alas, otherwise there + // are memory leaks; maybe there isanother way though. + wcsvfree(&nwcs, &wcs); // Do a second pass, in relaxed mode. We use the result. ret = wcspih(headers, number_of_records, WCSHDR_all, 0, &nreject, &nwcs, &wcs); @@ -221,11 +224,7 @@ void WCS::initFits(char* headers, int number_of_records) { // There are some things worth reporting about which WCS will not necessarily complain wcsCheckHeaders(wcs, headers, number_of_records); - m_wcs = decltype(m_wcs)(wcs, [nwcs](wcsprm* ptr) { - int nwcs_copy = nwcs; - wcsfree(ptr); - wcsvfree(&nwcs_copy, &ptr); - }); + m_wcs = make_wcsprm_ptr(wcs, false, nwcs); installSafeWcssub(); } @@ -241,8 +240,7 @@ void WCS::initAsdf(std::unique_ptr fits_wcs) { wcserr_enable(1); - wcsprm *wcs; - wcs = (wcsprm *)malloc(sizeof(*wcs)); + wcsprm *wcs = new wcsprm; if (!wcs) { throw Elements::Exception() << "failed to allocate memory for wcslib"; @@ -280,9 +278,7 @@ void WCS::initAsdf(std::unique_ptr fits_wcs) { wcsRaiseOnParseError(ret); wcsReportWarnings(wcsprintf_buf()); - m_wcs = decltype(m_wcs)(wcs, [](wcsprm* ptr) { - wcsfree(ptr); - }); + m_wcs = make_wcsprm_ptr(wcs, true); installSafeWcssub(); } @@ -297,14 +293,14 @@ void WCS::initAsdf(std::unique_ptr fits_wcs) { * Later we will figure out how to work in some config option(s) to explicitly * provide a path to the correct WCS to use if there is any ambiguity. */ -WCS::WCS(const AsdfImageSource& asdf_image_source) : m_wcs(nullptr, nullptr) { +WCS::WCS(const AsdfImageSource& asdf_image_source) { // First get whether we even have a WCS in the image auto fits_wcs = asdf_image_source.getFitsWCS(); initAsdf(std::move(fits_wcs)); } -WCS::WCS(const AsdfImageSource& asdf_image_source, std::optional wcs_path) : m_wcs(nullptr, nullptr) { +WCS::WCS(const AsdfImageSource& asdf_image_source, std::optional wcs_path) { // First get whether we even have a WCS in the image auto fits_wcs = asdf_image_source.getFitsWCS(wcs_path); initAsdf(std::move(fits_wcs)); @@ -327,17 +323,9 @@ WCS::~WCS() { WCS WCS::identity(int naxis) { - auto wcs = std::unique_ptr>( - new wcsprm, [](wcsprm* w) { - if (w) { - wcsfree(w); - delete w; - } - } - ); - - wcsini(1, naxis, wcs.get()); + WcsprmPtr wcs = make_wcsprm_ptr(); wcs->flag = -1; + wcsini(1, naxis, wcs.get()); for (int i = 0; i < naxis; i++) { wcs->crpix[i] = 1.0; wcs->crval[i] = 0.0; @@ -351,9 +339,9 @@ WCS WCS::identity(int naxis) { WorldCoordinate WCS::imageToWorld(ImageCoordinate image_coordinate) const { // wcsprm is in/out - wcsprm wcs_copy; - wcs_copy.flag = -1; - safe_wcssub(true, m_wcs.get(), nullptr, nullptr, &wcs_copy); + WcsprmPtr wcs_copy = make_wcsprm_ptr(); + wcs_copy->flag = -1; + safe_wcssub(true, m_wcs.get(), nullptr, nullptr, wcs_copy.get()); // +1 as fits standard coordinates start at 1 double pc_array[2] {image_coordinate.m_x + 1, image_coordinate.m_y + 1}; @@ -363,18 +351,17 @@ WorldCoordinate WCS::imageToWorld(ImageCoordinate image_coordinate) const { double phi, theta; int status = 0; - int ret_val = wcsp2s(&wcs_copy, 1, 1, pc_array, ic_array, &phi, &theta, wc_array, &status); - wcsRaiseOnTransformError(&wcs_copy, ret_val); - wcsfree(&wcs_copy); + int ret_val = wcsp2s(wcs_copy.get(), 1, 1, pc_array, ic_array, &phi, &theta, wc_array, &status); + wcsRaiseOnTransformError(wcs_copy.get(), ret_val); return WorldCoordinate(wc_array[0], wc_array[1]); } ImageCoordinate WCS::worldToImage(WorldCoordinate world_coordinate) const { // wcsprm is in/out - wcsprm wcs_copy; - wcs_copy.flag = -1; - safe_wcssub(true, m_wcs.get(), nullptr, nullptr, &wcs_copy); + WcsprmPtr wcs_copy = make_wcsprm_ptr(); + wcs_copy->flag = -1; + safe_wcssub(true, m_wcs.get(), nullptr, nullptr, wcs_copy.get()); double pc_array[2] {0, 0}; double ic_array[2] {0, 0}; @@ -382,13 +369,12 @@ ImageCoordinate WCS::worldToImage(WorldCoordinate world_coordinate) const { double phi, theta; int status = 0; - int ret_val = wcss2p(&wcs_copy, 1, 1, wc_array, &phi, &theta, ic_array, pc_array, &status); + int ret_val = wcss2p(wcs_copy.get(), 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(); } - wcsfree(&wcs_copy); return ImageCoordinate(pc_array[0] - 1, pc_array[1] - 1); // -1 as fits standard coordinates start at 1 } @@ -422,4 +408,24 @@ void WCS::addOffset(PixelCoordinate pc) { } +WCS::WcsprmPtr WCS::make_wcsprm_ptr() { + wcsprm *wcs = new wcsprm; + return WcsprmPtr(wcs, WcsprmDestroy{0, true}); +} + + +WCS::WcsprmPtr WCS::make_wcsprm_ptr(wcsprm *wcs) { + return WcsprmPtr(wcs, WcsprmDestroy{0, false}); +} + + +WCS::WcsprmPtr WCS::make_wcsprm_ptr(wcsprm *wcs, bool owned) { + return WcsprmPtr(wcs, WcsprmDestroy{0, owned}); +} + + +WCS::WcsprmPtr WCS::make_wcsprm_ptr(wcsprm *wcs, bool owned, int nwcs) { + return WcsprmPtr(wcs, WcsprmDestroy{nwcs, owned}); +} + } From ab4ab831081e5dea8e18e6363fa1c5ec10dddfe6 Mon Sep 17 00:00:00 2001 From: "E. Madison Bray" Date: Tue, 9 Dec 2025 11:34:28 +0100 Subject: [PATCH 29/29] delete stale TODO comment--this is now done --- SEFramework/src/lib/ASDF/AsdfImageSource.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/SEFramework/src/lib/ASDF/AsdfImageSource.cpp b/SEFramework/src/lib/ASDF/AsdfImageSource.cpp index d36afde13..ad10cc0e2 100644 --- a/SEFramework/src/lib/ASDF/AsdfImageSource.cpp +++ b/SEFramework/src/lib/ASDF/AsdfImageSource.cpp @@ -114,7 +114,6 @@ void AsdfImageSource::setLayer(int layer) { std::shared_ptr AsdfImageSource::getImageTile(int x, int y, int width, int height) const { - // TODO: (#5) support image data type conversion auto tile = ImageTile::create(m_image_type, x, y, width, height, std::const_pointer_cast(shared_from_this())); m_ndarray->fillImageTile(tile, m_current_layer);