fitsIoFunctions.hpp

Go to the documentation of this file.

 
#ifndef RTCTK_COMPONENTFRAMEWORK_FITSIOFUNCTIONS_HPP
#define RTCTK_COMPONENTFRAMEWORK_FITSIOFUNCTIONS_HPP
 
#include <rtctk/componentFramework/exceptions.hpp>
#include <rtctk/componentFramework/logger.hpp>
#include <rtctk/componentFramework/matrixBuffer.hpp>
#include <rtctk/componentFramework/typeTraits.hpp>
 
#include <fitsio.h>
 
#include <fmt/format.h>
 
#include <any>
#include <array>
#include <cassert>
#include <filesystem>
#include <limits>
#include <string>
#include <type_traits>
#include <typeinfo>
#include <unistd.h>
 
namespace {

template <typename T>
void IdentifyCfitsioTypes(int& bitpix, int& datatype) {
    if constexpr (std::is_integral_v<T>) {
        if constexpr (std::is_same_v<T, bool>) {
            bitpix = SBYTE_IMG;
            datatype = TSBYTE;
        } else if constexpr (sizeof(T) == 1) {
            if constexpr (std::is_signed_v<T>) {
                bitpix = SBYTE_IMG;
                datatype = TSBYTE;
            } else {
                bitpix = BYTE_IMG;
                datatype = TBYTE;
            }
        } else if constexpr (sizeof(T) == 2) {
            if constexpr (std::is_signed_v<T>) {
                bitpix = SHORT_IMG;
                datatype = TSHORT;
            } else {
                bitpix = USHORT_IMG;
                datatype = TUSHORT;
            }
        } else if constexpr (sizeof(T) == 4) {
            if constexpr (std::is_signed_v<T>) {
                bitpix = LONG_IMG;
                if constexpr (sizeof(int) == 4) {
                    datatype = TINT;
                } else if constexpr (sizeof(long) == 4) {
                    datatype = TLONG;
                } else {
                    static_assert(sizeof(int) == 4 or sizeof(long) == 4,
                                  "Require either int or long type to be 32 bits wide.");
                }
            } else {
                bitpix = ULONG_IMG;
                if constexpr (sizeof(unsigned int) == 4) {
                    datatype = TUINT;
                } else if constexpr (sizeof(unsigned long) == 4) {
                    datatype = TULONG;
                } else {
                    static_assert(
                        sizeof(unsigned int) == 4 or sizeof(unsigned long) == 4,
                        "Require either unsigned int or unsigned long type to be 32 bits wide.");
                }
            }
        } else if constexpr (sizeof(T) == 8) {
            if constexpr (std::is_signed_v<T>) {
                bitpix = LONGLONG_IMG;
                datatype = TLONGLONG;
            } else {
#ifdef ULONGLONG_IMG
                bitpix = ULONGLONG_IMG;
                datatype = TULONGLONG;
#else
                // If we are using an older Cfitsio library that does not have ULONGLONG_IMG defined
                // then simply alias to using the signed version for the image code and data type.
                bitpix = LONGLONG_IMG;
                datatype = TLONGLONG;
#endif
            }
        } else {
            static_assert(
                sizeof(T) == 1 or sizeof(T) == 2 or sizeof(T) == 4 or sizeof(T) == 8,
                "The integer type used as the template parameter must be one of 8, 16, 32, or 64"
                " bits wide.");
        }
    }
    if constexpr (std::is_floating_point_v<T>) {
        if constexpr (sizeof(T) == 4) {
            bitpix = FLOAT_IMG;
            datatype = TFLOAT;
        } else if constexpr (sizeof(T) == 8) {
            bitpix = DOUBLE_IMG;
            datatype = TDOUBLE;
        } else {
            static_assert(
                sizeof(T) == 4 or sizeof(T) == 8,
                "The floating-point type used as the template parameter must be float or double.");
        }
    } else {
        static_assert(std::is_arithmetic_v<T>,
                      "The template parameter used must be an integer or floating-point type.");
    }
}
 
}  // namespace
 
namespace rtctk::componentFramework {

std::string GetCfitsioErrorMsg(int status);

std::string CfitsioImageTypeToString(int bitpix);

std::string CfitsioDataTypeToString(int datatype);
 
// Some versions of GCC do not recognise that returning const std::type_info& does not result in
// a dangling reference. Therefore this is a false-positive warning.
#if defined(__GNUC__)
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wdangling-reference"
#endif

const std::type_info& GetFitsImageType(const std::string& filename);
 
#if defined(__GNUC__)
#pragma GCC diagnostic pop
#endif

std::vector<uint64_t> GetFitsImageShape(const std::string& filename);

std::size_t GetFitsImageSize(const std::string& filename);

template <typename T>
void WriteDataToFits(const std::string& filename, const T& buffer) {
    static_assert(
        IS_MATRIX_BUFFER_TYPE<T> || IS_MATRIX_SPAN_TYPE<T> || IS_VECTOR_TYPE<T> || IS_SPAN_TYPE<T>,
        "The buffer argument must be a matrix or a vector (or a corresponding span).");
 
    // Determine the image dimensions and total number of elements.
    std::vector<long> dims;
    long n_elements = buffer.size();
 
    if constexpr (IS_MATRIX_BUFFER_TYPE<T> || IS_MATRIX_SPAN_TYPE<T>) {
        dims.push_back(static_cast<long>(buffer.GetNcols()));
        dims.push_back(static_cast<long>(buffer.GetNrows()));
    } else {  // vector/span
        dims.push_back(static_cast<long>(buffer.size()));
    }
 
    using U = typename T::value_type;
    // Workout the appropriate Cfitsio image type code and data type code to use, based on the type
    // of the buffer elements, i.e. the type of U.
    int bitpix = 0, datatype = 0;
    IdentifyCfitsioTypes<U>(bitpix, datatype);
 
    std::string tmpfilename = filename + ".new-" + std::to_string(getpid());
 
    // Create the new FITS file and open it for modifications.
    fitsfile* fptr = nullptr;
    int status = 0;  // It is important to initialise this to zero before each Cfitsio API call.
    if (fits_create_file(&fptr, tmpfilename.c_str(), &status) != 0) {
        std::string msg = fmt::format(
            "Failed to create FITS file '{}'. {}", tmpfilename, GetCfitsioErrorMsg(status));
        CII_THROW(RtctkException, msg);
    }
 
    try {
        // Create the image header in the HDU.
        status = 0;
        if (fits_create_img(fptr, bitpix, dims.size(), dims.data(), &status) != 0) {
            std::string msg = fmt::format("Failed to create image in FITS file '{}'. {}",
                                          tmpfilename,
                                          GetCfitsioErrorMsg(status));
            CII_THROW(RtctkException, msg);
        }
 
        // Set the pixel start location to 1 for all dimensions. Cfitsio counts from 1, not 0.
        std::vector<long> fpixel(dims.size(), 1);
 
        // check bool by provided type in template
        if constexpr (std::is_same_v<U, bool>) {
            // For boolean element types we need to set the BOOLITEM keyword to indicate the special
            // data type handling.
            std::string value;
            value.reserve(FLEN_CARD);  // Must use maximum keyword card size according to Cfitsio.
            value = "true\0";
            assert(value.capacity() == FLEN_CARD);
            if (fits_write_key(fptr, TSTRING, "BOOLITEM", value.data(), nullptr, &status) != 0) {
                std::string msg = fmt::format("Failed to write keyword BOOLITEM to '{}'. {}",
                                              filename,
                                              GetCfitsioErrorMsg(status));
                CII_THROW(RtctkException, msg);
            }
 
            // For booleans we need to convert these into a type that can be handled by Cfitsio.
            std::vector<int8_t> tmp_buffer;
            tmp_buffer.resize(buffer.size());
            std::transform(
                buffer.begin(), buffer.end(), tmp_buffer.begin(), [](bool element) -> int8_t {
                    return element == true ? 1 : 0;
                });
 
            // We are forced to use const_cast, since the fits_write_pix function does not accept
            // const void*.
            void* array_ptr = const_cast<void*>(reinterpret_cast<const void*>(tmp_buffer.data()));
            status = 0;
            if (fits_write_pix(fptr, datatype, fpixel.data(), n_elements, array_ptr, &status) !=
                0) {
                std::string msg = fmt::format("Failed to write buffer to FITS file '{}'. {}",
                                              tmpfilename,
                                              GetCfitsioErrorMsg(status));
                CII_THROW(RtctkException, msg);
            }
        } else {
            // For all non-boolean element types we can write the buffer directly as pixel data.
            void* array_ptr = const_cast<void*>(reinterpret_cast<const void*>(buffer.data()));
            status = 0;
            if (fits_write_pix(fptr, datatype, fpixel.data(), n_elements, array_ptr, &status) !=
                0) {
                std::string msg = fmt::format("Failed to write buffer to FITS file '{}'. {}",
                                              tmpfilename,
                                              GetCfitsioErrorMsg(status));
                CII_THROW(RtctkException, msg);
            }
        }
    } catch (...) {
        // Attempt to close the FITS file if an exception was thrown. But do not throw an additional
        // exception if closing the file failed. Just log this error and rethrow the original
        // exception.
        status = 0;
        if (fits_close_file(fptr, &status) != 0) {
            LOG4CPLUS_ERROR(GetLogger("rtctk"),
                            fmt::format("Failed to close file '{}' while handling an exception. {}",
                                        tmpfilename,
                                        GetCfitsioErrorMsg(status)));
        }
        throw;
    }
 
    status = 0;
    if (fits_close_file(fptr, &status) != 0) {
        std::string msg = fmt::format(
            "Failed to close FITS file '{}'. {}", tmpfilename, GetCfitsioErrorMsg(status));
        CII_THROW(RtctkException, msg);
    }
 
    // Rename the temporary file to its final value. Updating the file in this manner will make sure
    // that different processes or threads will see changes to the file in an atomic manner. Either
    // the old file will be visible, or the new file, but not a partially modified file.
    try {
        std::filesystem::rename(tmpfilename, filename);
    } catch (const std::exception& error) {
        CII_THROW_WITH_NESTED(RtctkException,
                              error,
                              fmt::format("Failed to rename '{}' to '{}'.", tmpfilename, filename));
    }
}

template <typename T>
void ReadDataFromFits(const std::string& filename, T& buffer) {
    static_assert(
        IS_MATRIX_BUFFER_TYPE<T> || IS_MATRIX_SPAN_TYPE<T> || IS_VECTOR_TYPE<T> || IS_SPAN_TYPE<T>,
        "The buffer argument must be a matrix or a vector (or a corresponding span).");
 
    using U = typename T::value_type;
    // Identify the Cfitsio data type code to be used and the expected Cfitsio image type code,
    // based on the buffer element type, i.e. type of U.
    int expected_bitpix = 0;
    int datatype = 0;
    IdentifyCfitsioTypes<U>(expected_bitpix, datatype);
 
    fitsfile* fptr = nullptr;
    int status = 0;  // It is important to initialise this to zero before every Cfitsio API call.
    if (fits_open_image(&fptr, filename.c_str(), READONLY, &status) != 0) {
        std::string msg =
            fmt::format("Failed to open FITS file '{}'. {}", filename, GetCfitsioErrorMsg(status));
        CII_THROW(RtctkException, msg);
    }
 
    try {
        // Read and check that the pixel format used in the FITS file corresponds to the type used
        // for the buffer's elements passed to this function.
        int bitpix = -1;
        status = 0;
        if (fits_get_img_equivtype(fptr, &bitpix, &status) != 0) {
            std::string msg = fmt::format("Failed to read the image type from FITS file '{}'. {}",
                                          filename,
                                          GetCfitsioErrorMsg(status));
            CII_THROW(RtctkException, msg);
        }
 
        if (bitpix == LONGLONG_IMG) {
            // NOTE: Correcting the bitpix for ULONGLONG_IMG. There is a bug in Cfitsio that returns
            // LONGLONG_IMG instead of ULONGLONG_IMG.
            std::array<char, FLEN_CARD> value;  // Must use maximum possible FITS card size.
            status = 0;
            int result = fits_read_keyword(fptr, "BZERO", value.data(), nullptr, &status);
            if (result != 0 && status != KEY_NO_EXIST) {
                std::string msg = fmt::format("Failed to read keyword BZERO from '{}'. {}",
                                              filename,
                                              GetCfitsioErrorMsg(status));
                CII_THROW(RtctkException, msg);
            }
            if (std::string(value.data()) ==
                std::to_string(std::numeric_limits<uint64_t>::max() / 2 + 1)) {
#ifdef ULONGLONG_IMG
                bitpix = ULONGLONG_IMG;
#else
                // If we are dealing with an older Cfitsio library that does not have ULONGLONG_IMG
                // defined, then simply alias to using the signed version instead.
                bitpix = LONGLONG_IMG;
#endif
            }
        }
        if (bitpix != expected_bitpix) {
            std::string msg = fmt::format(
                "The FITS file '{}' has the wrong image format of {}. "
                "Expected a FITS image of type {}.",
                filename,
                CfitsioImageTypeToString(bitpix),
                CfitsioImageTypeToString(expected_bitpix));
            CII_THROW(RtctkException, msg);
        }
 
        // Read and check that the FITS file contains the expected number of image axes.
        int naxis = -1;
        status = 0;
        if (fits_get_img_dim(fptr, &naxis, &status) != 0) {
            std::string msg =
                fmt::format("Failed to read the number of image axes from FITS file '{}'. {}",
                            filename,
                            GetCfitsioErrorMsg(status));
            CII_THROW(RtctkException, msg);
        }
 
        constexpr bool is_matrix = IS_MATRIX_BUFFER_TYPE<T> or IS_MATRIX_SPAN_TYPE<T>;
        const int expected_naxis = is_matrix ? 2 : 1;
 
        if (naxis != expected_naxis) {
            std::string msg = fmt::format(
                "The FITS file '{}' has image dimensions that we cannot handle. "
                "Expected a {}D image.",
                filename,
                expected_naxis);
            CII_THROW(RtctkException, msg);
        }
 
        using size_type = typename T::size_type;
 
        // Read the shape information of the image data, i.e. the dimensions.
        std::vector<long> dims;
        dims.resize(static_cast<size_type>(naxis));
        status = 0;
        if (fits_get_img_size(fptr, naxis, dims.data(), &status) != 0) {
            std::string msg = fmt::format("Failed to read the image size from FITS file '{}'. {}",
                                          filename,
                                          GetCfitsioErrorMsg(status));
            CII_THROW(RtctkException, msg);
        }
 
        long n_elements = 1;
        for (auto dim : dims) {
            n_elements *= dim;
        }
 
        if constexpr (IS_VECTOR_TYPE<T>) {
            buffer.resize(static_cast<size_type>(dims[0]));
        } else if constexpr (IS_MATRIX_BUFFER_TYPE<T>) {
            auto nrows = static_cast<size_type>(dims[1]);
            auto ncols = static_cast<size_type>(dims[0]);
            buffer.resize(nrows, ncols);
        } else {
            if (buffer.size() < static_cast<size_type>(n_elements)) {
                CII_THROW(BufferTooSmall, buffer.size(), n_elements);
            }
        }
 
        // Set the pixel start location to 1 for all dimensions. Cfitsio counts from 1, not 0.
        std::vector<long> fpixel(dims.size(), 1);
        int anynull = 0;
 
        if constexpr (std::is_same_v<U, bool>) {
            // For booleans we need to first read the FITS data as integers, then convert these
            // into actual boolean values and fill the buffer.
            std::vector<int8_t> tmp_buffer;
            tmp_buffer.resize(buffer.size());
 
            status = 0;
            if (fits_read_pix(fptr,
                              datatype,
                              fpixel.data(),
                              n_elements,
                              nullptr,
                              tmp_buffer.data(),
                              &anynull,
                              &status) != 0) {
                std::string msg = fmt::format("Failed to read the image from FITS file '{}'. {}",
                                              filename,
                                              GetCfitsioErrorMsg(status));
                CII_THROW(RtctkException, msg);
            }
 
            std::transform(
                tmp_buffer.begin(), tmp_buffer.end(), buffer.begin(), [](int8_t element) {
                    return element == 0 ? false : true;
                });
        } else {
            // For non-booleans we read the FITS data directly into the buffer.
            status = 0;
            if (fits_read_pix(fptr,
                              datatype,
                              fpixel.data(),
                              n_elements,
                              nullptr,
                              buffer.data(),
                              &anynull,
                              &status) != 0) {
                std::string msg = fmt::format("Failed to read the image from FITS file '{}'. {}",
                                              filename,
                                              GetCfitsioErrorMsg(status));
                CII_THROW(RtctkException, msg);
            }
        }
    } catch (...) {
        // Attempt to close the FITS file if an exception was thrown. But do not throw an additional
        // exception if closing the file failed. Just log this error and rethrow the original
        // exception.
        status = 0;
        if (fits_close_file(fptr, &status) != 0) {
            LOG4CPLUS_ERROR(GetLogger("rtctk"),
                            fmt::format("Failed to close file '{}' while handling an exception. {}",
                                        filename,
                                        GetCfitsioErrorMsg(status)));
        }
        throw;
    }
 
    status = 0;
    if (fits_close_file(fptr, &status) != 0) {
        std::string msg = fmt::format(
            "Failed to close the FITS file '{}'. {}", filename, GetCfitsioErrorMsg(status));
        CII_THROW(RtctkException, msg);
    }
}
 
template <>
void WriteDataToFits(const std::string& filename, const std::any& buffer);
 
template <>
void ReadDataFromFits(const std::string& filename, std::any& buffer);

template <typename T>
[[deprecated("Replaced by WriteDataToFits")]]
void WriteMatrixToFits(const std::string& filename, const T& matrix, bool boolean_items = false) {
    static_assert(IS_MATRIX_BUFFER_TYPE<T> or IS_MATRIX_SPAN_TYPE<T>,
                  "The matrix argument must be of type MatrixBuffer or MatrixSpan.");
    assert(uint64_t(matrix.GetNrows()) <= uint64_t(std::numeric_limits<long>::max()));
    assert(uint64_t(matrix.GetNcols()) <= uint64_t(std::numeric_limits<long>::max()));
 
    using U = typename T::value_type;
    // Workout the appropriate Cfitsio image type code and data type code to use, based on the type
    // used with MatrixBuffer or MatrixSpan, i.e. the type of U.
    int bitpix = 0;
    int datatype = 0;
    IdentifyCfitsioTypes<U>(bitpix, datatype);
 
    std::string tmpfilename = filename + ".new-" + std::to_string(getpid());
 
    // Create the new FITS file and open it for modifications.
    fitsfile* fptr = nullptr;
    int status = 0;  // It is important to initialise this to zero before each Cfitsio API call.
    if (fits_create_file(&fptr, tmpfilename.c_str(), &status) != 0) {
        std::string msg = fmt::format(
            "Failed to create FITS file '{}'. {}", tmpfilename, GetCfitsioErrorMsg(status));
        CII_THROW(RtctkException, msg);
    }
 
    try {
        // Create the image header in the HDU.
        long nrows = static_cast<long>(matrix.GetNrows());
        long ncols = static_cast<long>(matrix.GetNcols());
        std::array<long, 2> shape = {ncols, nrows};
        status = 0;
        if (fits_create_img(fptr, bitpix, shape.size(), shape.data(), &status) != 0) {
            std::string msg = fmt::format("Failed to create the image type from FITS file '{}'. {}",
                                          tmpfilename,
                                          GetCfitsioErrorMsg(status));
            CII_THROW(RtctkException, msg);
        }
 
        if (boolean_items) {
            std::string value;
            value.reserve(FLEN_CARD);  // Must use maximum keyword card size according to Cfitsio.
            value = "true\0";
            assert(value.capacity() == FLEN_CARD);
            if (fits_write_key(fptr, TSTRING, "BOOLITEM", value.data(), nullptr, &status) != 0) {
                std::string msg = fmt::format("Failed to write keyword BOOLITEM to '{}'. {}",
                                              filename,
                                              GetCfitsioErrorMsg(status));
                CII_THROW(RtctkException, msg);
            }
        }
 
        // Write the matrix data as pixels.
        std::array<long, 2> fpixel = {1, 1};  // Pixel start location. Cfitsio counts from 1, not 0.
        long nelements = matrix.size();
        // We are forced to use const_cast, since the fits_write_pix function does not accept
        // const void*.
        void* array = const_cast<void*>(reinterpret_cast<const void*>(matrix.data()));
        status = 0;
        if (fits_write_pix(fptr, datatype, fpixel.data(), nelements, array, &status) != 0) {
            std::string msg = fmt::format("Failed to write the matrix to FITS file '{}'. {}",
                                          tmpfilename,
                                          GetCfitsioErrorMsg(status));
            CII_THROW(RtctkException, msg);
        }
    } catch (...) {
        // Attempt to close the FITS file if an exception was thrown. But do not throw an additional
        // exception if closing the file failed. Just log this error and rethrow the original
        // exception.
        status = 0;
        if (fits_close_file(fptr, &status) != 0) {
            LOG4CPLUS_ERROR(GetLogger("rtctk"),
                            fmt::format("Failed to close file '{}' while handling an exception. {}",
                                        tmpfilename,
                                        GetCfitsioErrorMsg(status)));
        }
        throw;
    }
 
    status = 0;
    if (fits_close_file(fptr, &status) != 0) {
        std::string msg = fmt::format(
            "Failed to close the FITS file '{}'. ", tmpfilename, GetCfitsioErrorMsg(status));
        CII_THROW(RtctkException, msg);
    }
 
    // Rename the temporary file to its final value. Updating the file in this manner will make sure
    // that different processes or threads will see changes to the file in an atomic manner. Either
    // the old file will be visible, or the new file, but not a partially modified file.
    try {
        std::filesystem::rename(tmpfilename, filename);
    } catch (const std::exception& error) {
        CII_THROW_WITH_NESTED(
            rtctk::componentFramework::RtctkException,
            error,
            fmt::format("Failed to rename FITS file '{}' to '{}'.", tmpfilename, filename))
    }
}

template <typename T>
[[deprecated("Replaced by ReadDataFromFits")]]
void ReadMatrixFromFits(const std::string& filename, T& matrix) {
    static_assert(IS_MATRIX_BUFFER_TYPE<T> or IS_MATRIX_SPAN_TYPE<T>,
                  "The matrix argument must be of type MatrixBuffer or MatrixSpan.");
 
    using U = typename T::value_type;
    // Identify the Cfitsio data type code to be used and the expected Cfitsio image type code,
    // based on the type used with MatrixBuffer or MatrixSpan, i.e. type of U.
    int expected_bitpix = 0;
    int datatype = 0;
    IdentifyCfitsioTypes<U>(expected_bitpix, datatype);
 
    fitsfile* fptr = nullptr;
    int status = 0;  // It is important to initialise this to zero before every Cfitsio API call.
    if (fits_open_image(&fptr, filename.c_str(), READONLY, &status) != 0) {
        std::string msg =
            fmt::format("Failed to open FITS file '{}'. {}", filename, GetCfitsioErrorMsg(status));
        CII_THROW(RtctkException, msg);
    }
 
    try {
        // Read and check that the pixel format used in the FITS file corresponds to the type used
        // for the matrix object passed to this function.
        int bitpix = -1;
        status = 0;
        if (fits_get_img_equivtype(fptr, &bitpix, &status) != 0) {
            std::string msg = fmt::format("Failed to read the image type from FITS file '{}'. {}",
                                          filename,
                                          GetCfitsioErrorMsg(status));
            CII_THROW(RtctkException, msg);
        }
        if (bitpix == LONGLONG_IMG) {
            // NOTE: Correcting the bitpix for ULONGLONG_IMG. There is a bug in Cfitsio that returns
            // LONGLONG_IMG instead of ULONGLONG_IMG.
            char value[80];  // Use maximum possible FITS card size of 80 bytes.
            status = 0;
            int result = fits_read_keyword(fptr, "BZERO", value, nullptr, &status);
            if (result != 0 and status != KEY_NO_EXIST) {
                std::string msg = fmt::format("Failed to read keyword BZERO from '{}'. {}",
                                              filename,
                                              GetCfitsioErrorMsg(status));
                CII_THROW(RtctkException, msg);
            }
            if (std::string(value) ==
                std::to_string(std::numeric_limits<uint64_t>::max() / 2 + 1)) {
#ifdef ULONGLONG_IMG
                bitpix = ULONGLONG_IMG;
#else
                // If we are dealing with an older Cfitsio library that does not have ULONGLONG_IMG
                // defined, then simply alias to using the signed version instead.
                bitpix = LONGLONG_IMG;
#endif
            }
        }
        if (bitpix != expected_bitpix) {
            std::string msg = fmt::format(
                "The FITS file '{}' has the wrong image format of {}. Expected a FITS image of "
                "type {}.",
                filename,
                CfitsioImageTypeToString(bitpix),
                CfitsioImageTypeToString(expected_bitpix));
            CII_THROW(RtctkException, msg);
        }
 
        // Read and check that the number of image axes is 2.
        int naxis = -1;
        status = 0;
        if (fits_get_img_dim(fptr, &naxis, &status) != 0) {
            std::string msg =
                fmt::format("Failed to read the number of image axes from FITS file '{}'. {}",
                            filename,
                            GetCfitsioErrorMsg(status));
            CII_THROW(RtctkException, msg);
        }
        if (naxis != 2) {
            std::string msg = fmt::format(
                "The FITS file '{}' has image dimensions that we cannot handle. Expect a 2D image "
                "when loading as a matrix.",
                filename);
            CII_THROW(RtctkException, msg);
        }
 
        long size[2] = {-1, -1};
        status = 0;
        if (fits_get_img_size(fptr, 2, size, &status) != 0) {
            std::string msg = fmt::format("Failed to read the image size from FITS file '{}'. {}",
                                          filename,
                                          GetCfitsioErrorMsg(status));
            CII_THROW(RtctkException, msg);
        }
 
        using size_type = typename T::size_type;
        auto nrows = static_cast<size_type>(size[1]);
        auto ncols = static_cast<size_type>(size[0]);
        long nelements = size[0] * size[1];  // Total number of pixels, i.e. rows * columns.
 
        long fpixel[2] = {1, 1};  // Pixel start location. Cfitsio counts from 1, not 0.
        int anynull = 0;  // Boolean flag indicating if there were any null/nil pixel values.
        if constexpr (IS_MATRIX_BUFFER_TYPE<T>) {
            matrix.resize(nrows, ncols);
        } else {
            if (matrix.size() < static_cast<size_type>(nrows * ncols)) {
                CII_THROW(BufferTooSmall, matrix.size(), nrows * ncols);
            }
        }
        void* array = matrix.data();
        status = 0;
        int result =
            fits_read_pix(fptr, datatype, fpixel, nelements, nullptr, array, &anynull, &status);
        if (result != 0) {
            std::string msg = fmt::format("Failed to read the image from FITS file '{}'. {}",
                                          filename,
                                          GetCfitsioErrorMsg(status));
            CII_THROW(RtctkException, msg);
        }
    } catch (...) {
        // Attempt to close the FITS file if an exception was thrown. But do not throw an additional
        // exception if closing the file failed. Just log this error and rethrow the original
        // exception.
        status = 0;
        if (fits_close_file(fptr, &status) != 0) {
            LOG4CPLUS_ERROR(GetLogger("rtctk"),
                            fmt::format("Failed to close file '{}' while handling an exception. {}",
                                        filename,
                                        GetCfitsioErrorMsg(status)));
        }
        throw;
    }
 
    status = 0;
    if (fits_close_file(fptr, &status) != 0) {
        std::string msg = fmt::format(
            "Failed to close the FITS file '{}'. {}", filename, GetCfitsioErrorMsg(status));
        CII_THROW(RtctkException, msg);
    }
}

template <typename T>
[[deprecated("Replaced by WriteDataToFits")]]
void WriteVectorToFits(const std::string& filename, const T& vector, bool boolean_items = false) {
    static_assert(IS_VECTOR_TYPE<T> or IS_SPAN_TYPE<T>,
                  "The vector argument must be of type std::vector or gsl::span.");
    assert(uint64_t(vector.size()) <= uint64_t(std::numeric_limits<long>::max()));
 
    using U = typename T::value_type;
    // Workout the appropriate Cfitsio image type code and data type code to use, based on the type
    // used with std::vector or gsl::span, i.e. the type of U.
    int bitpix = 0;
    int datatype = 0;
    IdentifyCfitsioTypes<U>(bitpix, datatype);
 
    std::string tmpfilename = filename + ".new-" + std::to_string(getpid());
 
    // Create the new FITS file and open it for modifications.
    fitsfile* fptr = nullptr;
    int status = 0;  // It is important to initialise this to zero before each Cfitsio API call.
    if (fits_create_file(&fptr, tmpfilename.c_str(), &status) != 0) {
        std::string msg = fmt::format(
            "Failed to create FITS file '{}'. {}", tmpfilename, GetCfitsioErrorMsg(status));
        CII_THROW(RtctkException, msg);
    }
 
    try {
        // Create the image header in the HDU.
        long size = static_cast<long>(vector.size());
        status = 0;
        if (fits_create_img(fptr, bitpix, 1, &size, &status) != 0) {
            std::string msg = fmt::format("Failed to create the image type from FITS file '{}'. {}",
                                          tmpfilename,
                                          GetCfitsioErrorMsg(status));
            CII_THROW(RtctkException, msg);
        }
 
        if (boolean_items) {
            std::string value;
            value.reserve(FLEN_CARD);  // Must use maximum keyword card size according to Cfitsio.
            value = "true\0";
            assert(value.capacity() == FLEN_CARD);
            if (fits_write_key(fptr, TSTRING, "BOOLITEM", value.data(), nullptr, &status) != 0) {
                std::string msg = fmt::format("Failed to write keyword BOOLITEM to '{}'. {}",
                                              filename,
                                              GetCfitsioErrorMsg(status));
                CII_THROW(RtctkException, msg);
            }
        }
 
        // Write the vector data as pixels.
        long fpixel = 1;  // Pixel start location. Cfitsio counts from 1, not 0.
        long nelements = vector.size();
        // We are forced to use const_cast, since the fits_write_pix function does not accept
        // const void*.
        void* array = const_cast<void*>(reinterpret_cast<const void*>(vector.data()));
        status = 0;
        if (fits_write_pix(fptr, datatype, &fpixel, nelements, array, &status) != 0) {
            std::string msg = fmt::format("Failed to write the vector to FITS file '{}'. {}",
                                          tmpfilename,
                                          GetCfitsioErrorMsg(status));
            CII_THROW(RtctkException, msg);
        }
    } catch (...) {
        // Attempt to close the FITS file if an exception was thrown. But do not throw an additional
        // exception if closing the file failed. Just log this error and rethrow the original
        // exception.
        status = 0;
        if (fits_close_file(fptr, &status) != 0) {
            LOG4CPLUS_ERROR(GetLogger("rtctk"),
                            fmt::format("Failed to close file '{}' while handling an exception. {}",
                                        tmpfilename,
                                        GetCfitsioErrorMsg(status)));
        }
        throw;
    }
 
    status = 0;
    if (fits_close_file(fptr, &status) != 0) {
        std::string msg = fmt::format(
            "Failed to close the FITS file '{}'. {}", tmpfilename, GetCfitsioErrorMsg(status));
        CII_THROW(RtctkException, msg);
    }
 
    // Rename the temporary file to its final value. Updating the file in this manner will make sure
    // that different processes or threads will see changes to the file in an atomic manner. Either
    // the old file will be visible, or the new file, but not a partially modified file.
    try {
        std::filesystem::rename(tmpfilename, filename);
    } catch (const std::exception& error) {
        CII_THROW_WITH_NESTED(
            rtctk::componentFramework::RtctkException,
            error,
            fmt::format("Failed to rename FITS file '{}' to '{}'.", tmpfilename, filename));
    }
}

template <typename T>
[[deprecated("Replaced by ReadDataFromFits")]]
void ReadVectorFromFits(const std::string& filename, T& vector) {
    static_assert(IS_VECTOR_TYPE<T> or IS_SPAN_TYPE<T>,
                  "The vector argument must be of type std::vector or gsl::span.");
 
    using U = typename T::value_type;
    // Identify the Cfitsio data type code to be used and the expected Cfitsio image type code,
    // based on the type used with std::vector or gsl::span, i.e. type of U.
    int expected_bitpix = 0;
    int datatype = 0;
    IdentifyCfitsioTypes<U>(expected_bitpix, datatype);
 
    fitsfile* fptr = nullptr;
    int status = 0;  // It is important to initialise this to zero before every Cfitsio API call.
    if (fits_open_image(&fptr, filename.c_str(), READONLY, &status) != 0) {
        std::string msg =
            fmt::format("Failed to open FITS file '{}'. {}", filename, GetCfitsioErrorMsg(status));
        CII_THROW(RtctkException, msg);
    }
 
    try {
        // Read and check that the pixel format used in the FITS file corresponds to the type used
        // for the vector object passed to this function.
        int bitpix = -1;
        status = 0;
        if (fits_get_img_equivtype(fptr, &bitpix, &status) != 0) {
            std::string msg = fmt::format("Failed to read the image type from FITS file '{}'. {}",
                                          filename,
                                          GetCfitsioErrorMsg(status));
            CII_THROW(RtctkException, msg);
        }
        if (bitpix == LONGLONG_IMG) {
            // NOTE: Correcting the bitpix for ULONGLONG_IMG. There is a bug in Cfitsio that returns
            // LONGLONG_IMG instead of ULONGLONG_IMG.
            char value[80];  // Use maximum possible FITS card size of 80 bytes.
            status = 0;
            int result = fits_read_keyword(fptr, "BZERO", value, nullptr, &status);
            if (result != 0 and status != KEY_NO_EXIST) {
                std::string msg = fmt::format("Failed to read keyword BZERO from '{}'. {}",
                                              filename,
                                              GetCfitsioErrorMsg(status));
                CII_THROW(RtctkException, msg);
            }
            if (std::string(value) ==
                std::to_string(std::numeric_limits<uint64_t>::max() / 2 + 1)) {
#ifdef ULONGLONG_IMG
                bitpix = ULONGLONG_IMG;
#else
                // If we are dealing with an older Cfitsio library that does not have ULONGLONG_IMG
                // defined, then simply alias to using the signed version instead.
                bitpix = LONGLONG_IMG;
#endif
            }
        }
        if (bitpix != expected_bitpix) {
            std::string msg = fmt::format(
                "The FITS file '{}' has the wrong image format of {}. Expected a FITS image of "
                "type {}.",
                filename,
                CfitsioImageTypeToString(bitpix),
                CfitsioImageTypeToString(expected_bitpix));
            CII_THROW(RtctkException, msg);
        }
 
        // Read and check that the number of image axes is 2.
        int naxis = -1;
        status = 0;
        if (fits_get_img_dim(fptr, &naxis, &status) != 0) {
            std::string msg =
                fmt::format("Failed to read the number of image axes from FITS file '{}'. {}",
                            filename,
                            GetCfitsioErrorMsg(status));
            CII_THROW(RtctkException, msg);
        }
        if (naxis != 1) {
            std::string msg = fmt::format(
                "The FITS file '{}' has image dimensions that we cannot handle. Expect a 1D image "
                "when loading as a vector.",
                filename);
            CII_THROW(RtctkException, msg);
        }
 
        long nelements = -1;
        status = 0;
        if (fits_get_img_size(fptr, 1, &nelements, &status) != 0) {
            std::string msg =
                fmt::format("Failed to read the 1D image size from FITS file '{}'. {}",
                            filename,
                            GetCfitsioErrorMsg(status));
            CII_THROW(RtctkException, msg);
        }
 
        using size_type = typename T::size_type;
        long fpixel = 1;  // Pixel start location. Cfitsio counts from 1, not 0.
        int anynull = 0;  // Boolean flag indicating if there were any null/nil pixel values.
        if constexpr (IS_VECTOR_TYPE<T>) {
            vector.resize(static_cast<size_type>(nelements));
        } else {
            if (vector.size() < static_cast<size_type>(nelements)) {
                CII_THROW(BufferTooSmall, vector.size(), nelements);
            }
        }
        void* array = vector.data();
        status = 0;
        int result =
            fits_read_pix(fptr, datatype, &fpixel, nelements, nullptr, array, &anynull, &status);
        if (result != 0) {
            std::string msg = fmt::format("Failed to read the image from FITS file '{}'. {}",
                                          filename,
                                          GetCfitsioErrorMsg(status));
            CII_THROW(RtctkException, msg);
        }
    } catch (...) {
        // Attempt to close the FITS file if an exception was thrown. But do not throw an additional
        // exception if closing the file failed. Just log this error and rethrow the original
        // exception.
        status = 0;
        if (fits_close_file(fptr, &status) != 0) {
            LOG4CPLUS_ERROR(GetLogger("rtctk"),
                            fmt::format("Failed to close file '{}' while handling an exception. {}",
                                        filename,
                                        GetCfitsioErrorMsg(status)));
        }
        throw;
    }
 
    status = 0;
    if (fits_close_file(fptr, &status) != 0) {
        std::string msg = fmt::format(
            "Failed to close the FITS file '{}'. {}", filename, GetCfitsioErrorMsg(status));
        CII_THROW(RtctkException, msg);
    }
}
 
// The following are template specialisations for vectors and matrices of boolean values.

template <typename A>
[[deprecated("Replaced by WriteDataToFits")]]
void WriteMatrixToFits(const std::string& filename, const MatrixBuffer<bool, A>& matrix) {
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wdeprecated-declarations"
    // Convert the boolean matrix to a matrix of 8 bit integers and write that to file instead.
    MatrixBuffer<int8_t> int_matrix;
    int_matrix.resize(matrix.GetNrows(), matrix.GetNcols());
    for (MatrixBuffer<int8_t>::size_type n = 0; n < int_matrix.GetNrows(); ++n) {
        for (MatrixBuffer<int8_t>::size_type m = 0; m < int_matrix.GetNcols(); ++m) {
            int_matrix(n, m) = matrix(n, m) ? 1 : 0;
        }
    }
    WriteMatrixToFits(filename, int_matrix, true);
#pragma GCC diagnostic pop
}

template <typename A>
[[deprecated("Replaced by ReadDataFromFits")]]
void ReadMatrixFromFits(const std::string& filename, MatrixBuffer<bool, A>& matrix) {
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wdeprecated-declarations"
    // Load the matrix from file as 8 bit integers and convert it to a boolean matrix.
    MatrixBuffer<int8_t> int_matrix;
    ReadMatrixFromFits(filename, int_matrix);
    matrix.resize(int_matrix.GetNrows(), int_matrix.GetNcols());
    for (MatrixBuffer<int8_t>::size_type n = 0; n < int_matrix.GetNrows(); ++n) {
        for (MatrixBuffer<int8_t>::size_type m = 0; m < int_matrix.GetNcols(); ++m) {
            matrix(n, m) = int_matrix(n, m) == 0 ? false : true;
        }
    }
#pragma GCC diagnostic pop
}

template <typename A>
[[deprecated("Replaced by WriteDataToFits")]]
void WriteVectorToFits(const std::string& filename, const std::vector<bool, A>& vector) {
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wdeprecated-declarations"
    // Concert the boolean vector to a vector of 8 bit integers and write that to file instead.
    std::vector<int8_t> int_vector;
    int_vector.resize(vector.size());
    for (std::vector<int8_t>::size_type n = 0; n < int_vector.size(); ++n) {
        int_vector[n] = vector[n] ? 1 : 0;
    }
    WriteVectorToFits(filename, int_vector, true);
#pragma GCC diagnostic pop
}

template <typename A>
[[deprecated("Replaced by ReadDataFromFits")]]
void ReadVectorFromFits(const std::string& filename, std::vector<bool, A>& vector) {
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wdeprecated-declarations"
    // Load the vector from file as 8 bit integers and convert it to a boolean vector.
    std::vector<int8_t> int_vector;
    ReadVectorFromFits(filename, int_vector);
    vector.resize(int_vector.size());
    for (std::vector<int8_t>::size_type n = 0; n < int_vector.size(); ++n) {
        vector[n] = int_vector[n] == 0 ? false : true;
    }
#pragma GCC diagnostic pop
}
 
}  // namespace rtctk::componentFramework
 
#endif  // RTCTK_COMPONENTFRAMEWORK_FITSIOFUNCTIONS_HPP