diff --git a/SEFramework/CMakeLists.txt b/SEFramework/CMakeLists.txt index 775552542..0d3172eee 100644 --- a/SEFramework/CMakeLists.txt +++ b/SEFramework/CMakeLists.txt @@ -31,10 +31,25 @@ 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) +#=============================================================================== +# 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) + add_definitions(-DWITH_ASDF) +endif() + + + #=============================================================================== # Declare the library dependencies here # Example: @@ -57,10 +72,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 +204,18 @@ 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) + 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 # Examples : diff --git a/SEFramework/SEFramework/ASDF/AsdfFile.h b/SEFramework/SEFramework/ASDF/AsdfFile.h new file mode 100644 index 000000000..b80d61397 --- /dev/null +++ b/SEFramework/SEFramework/ASDF/AsdfFile.h @@ -0,0 +1,229 @@ +/** 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 +#include + +#include + +#include "SEFramework/Image/ImageSourceWithMetadata.h" +#include "SEFramework/Image/ImageTile.h" + +namespace SourceXtractor { + +/** + * @class AsdfFile + * @brief represents access to a whole ASDF file + * + */ +class AsdfFile { +public: + /** + * 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 {}; + + /** + * 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) {} + + AsdfFile(AsdfFile&&) = default; + + virtual ~AsdfFile(); + + asdf_file_t* getAsdfFilePtr() const { return m_asdf_ptr.get(); } + + /** + * 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 { + 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. + */ + 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 * + */ + explicit Ndarray(const AsdfFile& file, asdf_value_t *ptr); + /** + * Private constructor for creating the `Ndarray` wrapper from a raw asdf_ndarray_t * + */ + explicit Ndarray(asdf_ndarray_t *ptr) + : 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; + }; + + /** + * 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_destroy(m_gwcs_ptr); + m_gwcs_fits_ptr = nullptr; + m_gwcs_ptr = nullptr; + } + + 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); + + asdf_gwcs_t *m_gwcs_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. + */ + std::unique_ptr 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. + */ + 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. */ + /** + * 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 { + 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; + + void open(); +}; + +} // namespace SourceXtractor + +#endif /* _SEFRAMEWORK_ASDF_ASDFFILE_H_ */ diff --git a/SEFramework/SEFramework/ASDF/AsdfImageSource.h b/SEFramework/SEFramework/ASDF/AsdfImageSource.h new file mode 100644 index 000000000..88725f9f9 --- /dev/null +++ b/SEFramework/SEFramework/ASDF/AsdfImageSource.h @@ -0,0 +1,161 @@ +/** + * 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 + */ +/* + * 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" +#include "SEFramework/Image/ImageTile.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, + ImageTile::ImageType image_type = ImageTile::AutoType, + std::shared_ptr manager = FileManager::getDefault()) : + AsdfImageSource(filename, image_index, std::nullopt, image_type, 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, + ImageTile::ImageType image_type = ImageTile::AutoType, + std::shared_ptr manager = FileManager::getDefault()) : + AsdfImageSource(filename, 0, image_path, image_type, 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; + } + + int getNdim() const { + return m_ndarray->ndim(); + } + + /** + * 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; + } + + /** + * 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_image_type; + } + + void setLayer(int layer) override; + + void saveTile(ImageTile& tile) override; + + /** 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, + 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; + 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..4671566dd --- /dev/null +++ b/SEFramework/SEFramework/ASDF/AsdfReader.h @@ -0,0 +1,59 @@ +/** Copyright © 2019-2025 Université de Genève, LMU Munich - Faculty of Physics, IAP-CNRS/Sorbonne Université + * + * This library is free software; you can redistribute it and/or modify it under + * the terms of the GNU Lesser General Public License as published by the Free + * Software Foundation; either version 3.0 of the License, or (at your option) + * any later version. + * + * This library is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + */ +/** + * @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/Image/ImageTile.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, + 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 */ + + +#endif /* _SEFRAMEWORK_ASDF_ASDFREADER_H */ + diff --git a/SEFramework/SEFramework/CoordinateSystem/WCS.h b/SEFramework/SEFramework/CoordinateSystem/WCS.h index 16f1ac8f7..389f018b4 100644 --- a/SEFramework/SEFramework/CoordinateSystem/WCS.h +++ b/SEFramework/SEFramework/CoordinateSystem/WCS.h @@ -27,8 +27,16 @@ #include #include +#include +#include + #include "SEFramework/CoordinateSystem/CoordinateSystem.h" #include "SEFramework/FITS/FitsImageSource.h" +#include "SEFramework/Image/ImageSource.h" +#ifdef WITH_ASDF +#include "SEFramework/ASDF/AsdfFile.h" +#include "SEFramework/ASDF/AsdfImageSource.h" +#endif struct wcsprm; @@ -37,10 +45,18 @@ 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 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,9 +65,43 @@ class WCS : public CoordinateSystem { void addOffset(PixelCoordinate pc); private: - void init(char* headers, int number_of_records); + void initFits(char* headers, int number_of_records); + +#ifdef WITH_ASDF + void initAsdf(std::unique_ptr fits_wcs); +#endif + + 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)) {} - std::unique_ptr> m_wcs; }; } diff --git a/SEFramework/SEFramework/FITS/FitsImageSource.h b/SEFramework/SEFramework/FITS/FitsImageSource.h index a66107d63..da186d359 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 @@ -58,7 +71,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, @@ -88,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; @@ -120,7 +139,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..b9a099048 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,15 +15,16 @@ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ /** - * @file SEFramework/Image/FitsReader.h - * @date 06/14/16 - * @author nikoapos + * @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" +#include "SEFramework/Image/ImageTile.h" #include "SEFramework/FITS/FitsImageSource.h" namespace SourceXtractor { @@ -33,24 +34,38 @@ namespace SourceXtractor { * @brief * */ -template -class FitsReader { +class FitsReader : public ImageFileReader { public: + using ImageFileReader::ImageFileReader; + static bool test(std::istream& stream); + + using ImageFileReader::get; /** - * @brief Destructor + * 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 */ - virtual ~FitsReader() = default; + 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, + ImageTile::ImageType image_type = ImageTile::AutoType); - 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; }; /* End of FitsReader class */ } /* namespace SourceXtractor */ -#endif +#endif /* _SEFRAMEWORK_FITS_FITSREADER_H */ + diff --git a/SEFramework/SEFramework/Image/ImageFileReader.h b/SEFramework/SEFramework/Image/ImageFileReader.h new file mode 100644 index 000000000..411b4864b --- /dev/null +++ b/SEFramework/SEFramework/Image/ImageFileReader.h @@ -0,0 +1,267 @@ +/** 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/BufferedImage.h" +#include "SEFramework/Image/ImageSource.h" +#include "SEFramework/Image/ImageTile.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; + + /** + * 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 { + 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() noexcept : m_reader(nullptr), m_index(0), m_image_type(ImageTile::AutoType) {} + + Iterator(ImageFileReader* reader, int index, + ImageTile::ImageType image_type = ImageTile::AutoType) + : m_reader(reader), m_index(index), m_image_type(image_type) { + 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); + } + + /* 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" + if (nullptr == m_reader) { + m_current = nullptr; + return; + } + + try { + m_current = m_reader->get(m_index, m_image_type); + } catch (...) { + m_current = nullptr; + } + } + + 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, ImageTile::AutoType); + } + + Iterator end() { + return Iterator(); + } + + /** + * 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(); + + /** + * 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) + * @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/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/SEFramework/auxdir/multiple_hdu.asdf b/SEFramework/auxdir/multiple_hdu.asdf new file mode 100644 index 000000000..8cca55c74 Binary files /dev/null and b/SEFramework/auxdir/multiple_hdu.asdf differ diff --git a/SEFramework/auxdir/wcs_header.asdf b/SEFramework/auxdir/wcs_header.asdf new file mode 100644 index 000000000..9f5bad3b8 Binary files /dev/null and b/SEFramework/auxdir/wcs_header.asdf differ diff --git a/SEFramework/auxdir/with_primary.asdf b/SEFramework/auxdir/with_primary.asdf new file mode 100644 index 000000000..14622e356 Binary files /dev/null and b/SEFramework/auxdir/with_primary.asdf differ diff --git a/SEFramework/src/lib/ASDF/AsdfFile.cpp b/SEFramework/src/lib/ASDF/AsdfFile.cpp new file mode 100644 index 000000000..9decbd4e3 --- /dev/null +++ b/SEFramework/src/lib/ASDF/AsdfFile.cpp @@ -0,0 +1,366 @@ +/** 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, 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() {} + + +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: "; + } + + 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 +std::unique_ptr 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) { + // 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++; + } + } + + throw AsdfValueNotFoundException() << "No ndarray at index " << index << " in file: " + << 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; + } + + bool is_fits = asdf_gwcs_is_fits((asdf_file_t*)asdf_value_file(value), gwcs); + asdf_gwcs_destroy(gwcs); + return is_fits; +} + + +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_value_clone(asdf_find_item_value(item)); + asdf_find_item_destroy(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()); + if (!v) { + throw AsdfValueNotFoundException() << "No value at path " << path << " in file: " + << m_path.native(); + } + + return AsdfValuePtr(v); +} + + +static ImageTile::ImageType convertDatatypeToImageType(const asdf_datatype_t* datatype) { + ImageTile::ImageType image_type; + + switch (datatype->type) { + 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->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; + 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 + break; + case ASDF_VALUE_ERR_TYPE_MISMATCH: { + 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; + m_image_type = convertDatatypeToImageType(&ndarray_ptr->datatype); +} + + +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, + convertImageTypeToDatatype(image_tile->getType()), + &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"; + } +} + + +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; +} + + +/** + * 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) { + 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_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); +} + + +} // namespace SourceXtractor diff --git a/SEFramework/src/lib/ASDF/AsdfImageSource.cpp b/SEFramework/src/lib/ASDF/AsdfImageSource.cpp new file mode 100644 index 000000000..ad10cc0e2 --- /dev/null +++ b/SEFramework/src/lib/ASDF/AsdfImageSource.cpp @@ -0,0 +1,201 @@ +/** + * 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 + +#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, + 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)) { + + 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_image_type, 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"; +}; + + +/** + * 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; + 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; + + 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/src/lib/ASDF/AsdfReader.cpp b/SEFramework/src/lib/ASDF/AsdfReader.cpp new file mode 100644 index 000000000..3f55369c0 --- /dev/null +++ b/SEFramework/src/lib/ASDF/AsdfReader.cpp @@ -0,0 +1,66 @@ +/** 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, ImageTile::ImageType image_type) { + return std::make_shared(m_filename, image_index, image_type); +} + + +std::shared_ptr AsdfReader::get(const std::string& extname, + ImageTile::ImageType image_type) { + return std::make_shared(m_filename, extname, image_type); +} + + +static ImageFileReader::StaticFileType asdf_reader{ + "asdf", + {"asdf"} +}; + +} + diff --git a/SEFramework/src/lib/CoordinateSystem/WCS.cpp b/SEFramework/src/lib/CoordinateSystem/WCS.cpp index ac14c0612..07b67afec 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" @@ -42,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 */ @@ -140,30 +169,22 @@ 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) { +WCS::WCS(const FitsImageSource& fits_image_source) { 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) { +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 + // (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 +192,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; @@ -191,6 +212,9 @@ void WCS::init(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); @@ -200,30 +224,124 @@ void WCS::init(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); - 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 WITH_ASDF +void WCS::initAsdf(std::unique_ptr fits_wcs) { + if (!fits_wcs) { + auto tmp = WCS::identity(2); + m_wcs = std::move(tmp.m_wcs); + return; } + + wcserr_enable(1); + + wcsprm *wcs = new wcsprm; + + 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 = make_wcsprm_ptr(wcs, true); + + 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) { + // 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) { + // 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 */ + + +/** + * 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() { } + +WCS WCS::identity(int naxis) { + 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; + 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; - 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}; @@ -233,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}; @@ -252,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 } @@ -292,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}); +} + } diff --git a/SEFramework/src/lib/FITS/FitsImageSource.cpp b/SEFramework/src/lib/FITS/FitsImageSource.cpp index 776ac420b..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; @@ -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; + throw UnknownHduException() << "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 UnknownHduException() << "Can't get the active HDU from the FITS file: " << filename; + } + } else { switchHdu(fptr, m_hdu_number); } @@ -104,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; } @@ -161,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; } @@ -270,14 +275,29 @@ 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; + } +} + + +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 UnknownHduException() << "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..f608f0f05 --- /dev/null +++ b/SEFramework/src/lib/FITS/FitsReader.cpp @@ -0,0 +1,100 @@ +/** 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/FITS/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(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) + // so should just use that instead. + auto it = m_image_hdu_map.find(image_index); + + if (it != m_image_hdu_map.end()) { + return std::make_shared(m_filename, it->second, image_type); + } + + // 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 { + auto image_source = getHdu(++hdu_num, image_type); + 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; + } + } + + // Shouldn't get here + throw FitsImageSource::UnknownHduException(); +} + + +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, ImageTile::ImageType image_type) { + return std::make_shared(m_filename, hdu_num, image_type); +} + + +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..c1c3f905e --- /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, ImageTile::AutoType); + } else if (m_image_path) { + return get(*m_image_path, ImageTile::AutoType); + } + + 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/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/AsdfFile_test.cpp b/SEFramework/tests/src/ASDF/AsdfFile_test.cpp new file mode 100644 index 000000000..b973f930a --- /dev/null +++ b/SEFramework/tests/src/ASDF/AsdfFile_test.cpp @@ -0,0 +1,137 @@ +/** 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/AsdfFile_test.cpp + * @date 08/14/25 + * @author embray + */ + +#include + +#include +#include + +#include "SEFramework/ASDF/AsdfFile.h" + +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(); + } +}; + +//----------------------------------------------------------------------------- + +BOOST_AUTO_TEST_SUITE (AsdfFile_test) + + +//----------------------------------------------------------------------------- + +BOOST_FIXTURE_TEST_CASE( open_file, AsdfFileFixture ) { + BOOST_CHECK_NO_THROW({ + AsdfFile asdf_file(primary_path); + }); +} + +//----------------------------------------------------------------------------- + +BOOST_AUTO_TEST_CASE( missing_file ) { + BOOST_CHECK_THROW({ + AsdfFile asdf_file("/does/not/exist"); + }, Elements::Exception); +} + +//----------------------------------------------------------------------------- + +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(), 2); + auto shape = ndarray->shape(); + std::vector expected_shape{1, 1}; + BOOST_CHECK_EQUAL_COLLECTIONS( + shape.begin(), shape.end(), + expected_shape.begin(), expected_shape.end() + ); +} + +//----------------------------------------------------------------------------- + +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{1, 1}; + BOOST_CHECK_EQUAL_COLLECTIONS( + shape.begin(), shape.end(), + expected_shape.begin(), expected_shape.end() + ); + BOOST_CHECK_EQUAL(ndarray->getPath(), "/PRIMARY"); +} + + +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() + ); +} + +//----------------------------------------------------------------------------- + + +BOOST_AUTO_TEST_SUITE_END () diff --git a/SEFramework/tests/src/ASDF/AsdfImageSource_test.cpp b/SEFramework/tests/src/ASDF/AsdfImageSource_test.cpp new file mode 100644 index 000000000..e63fe0abb --- /dev/null +++ b/SEFramework/tests/src/ASDF/AsdfImageSource_test.cpp @@ -0,0 +1,159 @@ +/** Copyright © 2020 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/AsdfImageSource_test.cpp + * @date 10/03/2025 + * @author embray + */ + +#include + +#include + +#include "SEFramework/ASDF/AsdfImageSource.h" +#include "SEFramework/Image/ImageAccessor.h" + +using namespace SourceXtractor; + +struct AsdfImageSourceFixture { + 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(); + } +}; + +//----------------------------------------------------------------------------- + +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_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() + +//----------------------------------------------------------------------------- + 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 () 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() //----------------------------------------------------------------------------- diff --git a/SEFramework/tests/src/FITS/FitsReader_test.cpp b/SEFramework/tests/src/FITS/FitsReader_test.cpp index f667773ec..c90e5dd32 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::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); @@ -66,8 +68,42 @@ BOOST_FIXTURE_TEST_CASE ( image_source, FitsReaderFixture ) { //----------------------------------------------------------------------------- +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 ( 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( + 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::readImage("/not/existing/path"), Elements::Exception); } //----------------------------------------------------------------------------- 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 54241f27d..92896a120 100644 --- a/SEImplementation/SEImplementation/Configuration/DetectionImageConfig.h +++ b/SEImplementation/SEImplementation/Configuration/DetectionImageConfig.h @@ -24,9 +24,13 @@ #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" +#ifdef WITH_ASDF +#include "SEFramework/ASDF/AsdfImageSource.h" +#endif namespace SourceXtractor { @@ -47,7 +51,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 +59,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 +94,25 @@ 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); +#ifdef WITH_ASDF + DetectionImageExtension(std::shared_ptr fits_image_source, double gain, + double saturation, double flux_scale, int interpolation_gap, + std::optional wcs_path); +#endif + + 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..2d47ba100 100644 --- a/SEImplementation/src/lib/Configuration/DetectionImageConfig.cpp +++ b/SEImplementation/src/lib/Configuration/DetectionImageConfig.cpp @@ -22,15 +22,20 @@ #include "Configuration/ConfigManager.h" #include +#include using boost::regex; using boost::regex_match; using boost::smatch; -#include -#include +#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" +#ifdef WITH_ASDF +#include "SEFramework/ASDF/AsdfImageSource.h" +#endif #include "SEImplementation/Configuration/DetectionImageConfig.h" @@ -47,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(), @@ -66,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 }}}; } @@ -78,9 +107,21 @@ 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); +#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; @@ -95,109 +136,146 @@ 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->iter(ImageTile::FloatImage)) { + 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(*image_source); + 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)); +#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) { + 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); + } +#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); +} + + +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; } diff --git a/SEImplementation/src/lib/Configuration/MeasurementImageConfig.cpp b/SEImplementation/src/lib/Configuration/MeasurementImageConfig.cpp index ba8140805..0f7c79519 100644 --- a/SEImplementation/src/lib/Configuration/MeasurementImageConfig.cpp +++ b/SEImplementation/src/lib/Configuration/MeasurementImageConfig.cpp @@ -32,8 +32,9 @@ #include #include +#include +#include #include -#include #include #include @@ -81,8 +82,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 +118,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 +153,7 @@ WeightImage::PixelType extractWeightThreshold(const PyMeasurementImage& py_image } else { threshold = std::numeric_limits::max(); } - break; + break; } return threshold; } @@ -185,15 +187,15 @@ 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); + 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; diff --git a/SEImplementation/src/lib/Configuration/SegmentationConfig.cpp b/SEImplementation/src/lib/Configuration/SegmentationConfig.cpp index 04a6438f9..8dcc582a5 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::readImage(filename); // give some feedback on the filter segConfigLogger.info() << "Loaded segmentation filter: " << filename << " height: " << convolution_kernel->getHeight() << " width: " << convolution_kernel->getWidth(); 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); 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