fitsIoFunctions.hpp

Go to the documentation of this file.

 
#ifndef RTCTK_COMPONENTFRAMEWORK_FITSIOFUNCTIONS_HPP
#define RTCTK_COMPONENTFRAMEWORK_FITSIOFUNCTIONS_HPP
 
#include <cassert>
#include <type_traits>
#include <limits>
#include <string>
#include <boost/filesystem.hpp>
#include <unistd.h>
#include <fitsio.h>
#include <rtctk/componentFramework/logger.hpp>
#include <rtctk/componentFramework/exceptions.hpp>
#include <rtctk/componentFramework/matrixBuffer.hpp>
 
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 = BYTE_IMG;
            datatype = TBYTE;
        } 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 anonymous
 
namespace rtctk::componentFramework {
 
std::string GetCfitsioErrorMsg(int status);
std::string CfitsioImageTypeToString(int bitpix);
std::string CfitsioDataTypeToString(int datatype);
 
template <typename T, typename A>
void WriteMatrixToFits(const std::string& filename, const MatrixBuffer<T, A>& matrix) {
    assert(int64_t(matrix.get_nrows()) <= int64_t(std::numeric_limits<long>::max()));
    assert(int64_t(matrix.get_ncols()) <= int64_t(std::numeric_limits<long>::max()));
 
    // Workout the appropriate Cfitsio image type code and data type code to use, based on the type
    // used with MatrixBuffer, i.e. the type of T.
    int bitpix = 0;
    int datatype = 0;
    IdentifyCfitsioTypes<T>(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 = "Failed to create FITS file '" + tmpfilename + "'. "
            + GetCfitsioErrorMsg(status);
        CII_THROW(std::runtime_error, msg);
    }
 
    try {
        // Create the image header in the HDU.
        long nrows = long(matrix.get_nrows());
        long ncols = long(matrix.get_ncols());
        long size[2] = {ncols, nrows};
        status = 0;
        if (fits_create_img(fptr, bitpix, 2, size, &status) != 0) {
            std::string msg = "Failed to create the image type from FITS file '" + tmpfilename
                + "'. " + GetCfitsioErrorMsg(status);
            CII_THROW(std::runtime_error, msg);
        }
 
        // Write the matrix data as pixels.
        long fpixel[2] = {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, nelements, array, &status) != 0) {
            std::string msg = "Failed to write the matrix to FITS file '" + tmpfilename + "'. "
                + GetCfitsioErrorMsg(status);
            CII_THROW(std::runtime_error, 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(), "Failed to close file '" + tmpfilename
                + "' while handling an exception. " + GetCfitsioErrorMsg(status));
        }
        throw;
    }
 
    status = 0;
    if (fits_close_file(fptr, &status) != 0) {
        std::string msg = "Failed to close the FITS file '" + tmpfilename + "'. "
            + GetCfitsioErrorMsg(status);
        CII_THROW(std::runtime_error, 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 {
        boost::filesystem::rename(tmpfilename, filename);
    } catch (const std::exception& error) {
        CII_THROW_WITH_NESTED(std::runtime_error, error,
            "Failed to rename FITS file '" + tmpfilename + "' to '" + filename + "'.")
    }
}
 
template <typename T, typename A>
void ReadMatrixFromFits(const std::string& filename, MatrixBuffer<T, A>& matrix) {
    // Identify the Cfitsio data type code to be used and the expected Cfitsio image type code,
    // based on the type used with MatrixBuffer, i.e. type of T.
    int expected_bitpix = 0;
    int datatype = 0;
    IdentifyCfitsioTypes<T>(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 = "Failed to open FITS file '" + filename + "'. "
            + GetCfitsioErrorMsg(status);
        CII_THROW(std::runtime_error, 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 = "Failed to read the image type from FITS file '" + filename + "'. "
                + GetCfitsioErrorMsg(status);
            CII_THROW(std::runtime_error, 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 = "Failed to read keyword BZERO from '" + filename + "'. "
                    + GetCfitsioErrorMsg(status);
                CII_THROW(std::runtime_error, 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 = "The FITS file '" + filename + "' has the wrong image format of " +
                CfitsioImageTypeToString(bitpix) + ". Expected a"
                " FITS image of type " + CfitsioImageTypeToString(expected_bitpix) + ".";
            CII_THROW(std::runtime_error, 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 = "Failed to read the number of image axes from FITS file '"
                + filename + "'. " + GetCfitsioErrorMsg(status);
            CII_THROW(std::runtime_error, msg);
        }
        if (naxis != 2) {
            std::string msg = "The FITS file '" + filename + "' has image dimensions that we cannot"
                " handle. Expect a 2D image when loading as a matrix.";
            CII_THROW(std::runtime_error, msg);
        }
 
        long size[2] = {-1, -1};
        status = 0;
        if (fits_get_img_size(fptr, 2, size, &status) != 0) {
            std::string msg = "Failed to read the image size from FITS file '" + filename + "'. "
                + GetCfitsioErrorMsg(status);
            CII_THROW(std::runtime_error, msg);
        }
 
        using size_type = typename MatrixBuffer<T, A>::size_type;
        auto nrows = size_type(size[1]);
        auto ncols = 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.
        matrix.resize(nrows, ncols);
        void* array = matrix.data();
        status = 0;
        if (fits_read_pix(fptr, datatype, fpixel, nelements, nullptr, array, &anynull,
                          &status) != 0) {
            std::string msg = "Failed to read the image from FITS file '" + filename + "'. "
                + GetCfitsioErrorMsg(status);
            CII_THROW(std::runtime_error, 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(), "Failed to close file '" + filename
                + "' while handling an exception. "+ GetCfitsioErrorMsg(status));
        }
        throw;
    }
 
    status = 0;
    if (fits_close_file(fptr, &status) != 0) {
        std::string msg = "Failed to close the FITS file '" + filename + "'. "
            + GetCfitsioErrorMsg(status);
        CII_THROW(std::runtime_error, msg);
    }
}
 
template <typename T, typename A>
void WriteVectorToFits(const std::string& filename, const std::vector<T, A>& vector) {
    assert(int64_t(vector.size()) <= int64_t(std::numeric_limits<long>::max()));
 
    // Workout the appropriate Cfitsio image type code and data type code to use, based on the type
    // used with std::vector, i.e. the type of T.
    int bitpix = 0;
    int datatype = 0;
    IdentifyCfitsioTypes<T>(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 = "Failed to create FITS file '" + tmpfilename + "'. "
            + GetCfitsioErrorMsg(status);
        CII_THROW(std::runtime_error, msg);
    }
 
    try {
        // Create the image header in the HDU.
        long size = long(vector.size());
        status = 0;
        if (fits_create_img(fptr, bitpix, 1, &size, &status) != 0) {
            std::string msg = "Failed to create the image type from FITS file '" + tmpfilename
                + "'. " + GetCfitsioErrorMsg(status);
            CII_THROW(std::runtime_error, 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 = "Failed to write the vector to FITS file '" + tmpfilename + "'. "
                + GetCfitsioErrorMsg(status);
            CII_THROW(std::runtime_error, 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(), "Failed to close file '" + tmpfilename
                + "' while handling an exception. " + GetCfitsioErrorMsg(status));
        }
        throw;
    }
 
    status = 0;
    if (fits_close_file(fptr, &status) != 0) {
        std::string msg = "Failed to close the FITS file '" + tmpfilename + "'. "
            + GetCfitsioErrorMsg(status);
        CII_THROW(std::runtime_error, 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 {
        boost::filesystem::rename(tmpfilename, filename);
    } catch (const std::exception& error) {
        CII_THROW_WITH_NESTED(std::runtime_error, error,
            "Failed to rename FITS file '" + tmpfilename + "' to '" + filename + "'.")
    }
}
 
template <typename T, typename A>
void ReadVectorFromFits(const std::string& filename, std::vector<T, A>& vector) {
    // Identify the Cfitsio data type code to be used and the expected Cfitsio image type code,
    // based on the type used with std::vector, i.e. type of T.
    int expected_bitpix = 0;
    int datatype = 0;
    IdentifyCfitsioTypes<T>(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 = "Failed to open FITS file '" + filename + "'. "
            + GetCfitsioErrorMsg(status);
        CII_THROW(std::runtime_error, 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 = "Failed to read the image type from FITS file '" + filename + "'. "
                + GetCfitsioErrorMsg(status);
            CII_THROW(std::runtime_error, 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 = "Failed to read keyword BZERO from '" + filename + "'. "
                    + GetCfitsioErrorMsg(status);
                CII_THROW(std::runtime_error, 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 = "The FITS file '" + filename + "' has the wrong image format of " +
                CfitsioImageTypeToString(bitpix) + ". Expected a"
                " FITS image of type " + CfitsioImageTypeToString(expected_bitpix) + ".";
            CII_THROW(std::runtime_error, 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 = "Failed to read the number of image axes from FITS file '"
                + filename + "'. " + GetCfitsioErrorMsg(status);
            CII_THROW(std::runtime_error, msg);
        }
        if (naxis != 1) {
            std::string msg = "The FITS file '" + filename + "' has image dimensions that we cannot"
                " handle. Expect a 1D image when loading as a vector.";
            CII_THROW(std::runtime_error, msg);
        }
 
        long nelements = -1;
        status = 0;
        if (fits_get_img_size(fptr, 1, &nelements, &status) != 0) {
            std::string msg = "Failed to read the 1D image size from FITS file '" + filename + "'. "
                + GetCfitsioErrorMsg(status);
            CII_THROW(std::runtime_error, msg);
        }
 
        using size_type = typename std::vector<T, A>::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.
        vector.resize(size_type(nelements));
        void* array = vector.data();
        status = 0;
        if (fits_read_pix(fptr, datatype, &fpixel, nelements, nullptr, array, &anynull,
                          &status) != 0) {
            std::string msg = "Failed to read the image from FITS file '" + filename + "'. "
                + GetCfitsioErrorMsg(status);
            CII_THROW(std::runtime_error, 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(), "Failed to close file '" + filename
                + "' while handling an exception. "+ GetCfitsioErrorMsg(status));
        }
        throw;
    }
 
    status = 0;
    if (fits_close_file(fptr, &status) != 0) {
        std::string msg = "Failed to close the FITS file '" + filename + "'. "
            + GetCfitsioErrorMsg(status);
        CII_THROW(std::runtime_error, msg);
    }
}
 
// The following are template specialisations for vectors and matrices of boolean values.
 
template <typename A>
void WriteMatrixToFits(const std::string& filename, const MatrixBuffer<bool, A>& matrix) {
    // Concert 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.get_nrows(), matrix.get_ncols());
    for (MatrixBuffer<int8_t>::size_type n = 0; n < int_matrix.get_nrows(); ++n) {
        for (MatrixBuffer<int8_t>::size_type m = 0; m < int_matrix.get_ncols(); ++m) {
            int_matrix(n, m) = matrix(n, m) ? 1 : 0;
        }
    }
    WriteMatrixToFits(filename, int_matrix);
}
 
template <typename A>
void ReadMatrixFromFits(const std::string& filename, MatrixBuffer<bool, A>& matrix) {
    // 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.get_nrows(), int_matrix.get_ncols());
    for (MatrixBuffer<int8_t>::size_type n = 0; n < int_matrix.get_nrows(); ++n) {
        for (MatrixBuffer<int8_t>::size_type m = 0; m < int_matrix.get_ncols(); ++m) {
            matrix(n, m) = int_matrix(n, m) == 1 ? true : false;
        }
    }
}
 
template <typename A>
void WriteVectorToFits(const std::string& filename, const std::vector<bool, A>& vector) {
    // 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);
}
 
template <typename A>
void ReadVectorFromFits(const std::string& filename, std::vector<bool, A>& vector) {
    // 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] == 1 ? true : false;
    }
}
 
} // namespace rtctk::componentFramework
 
#endif // RTCTK_COMPONENTFRAMEWORK_FITSIOFUNCTIONS_HPP