fits.h

Go to the documentation of this file.

// -*- lsst-c++ -*-
#ifndef LSST_AFW_fits_h_INCLUDED
#define LSST_AFW_fits_h_INCLUDED
 
/*
 *  Utilities for working with FITS files.  These are mostly thin wrappers around
 *  cfitsio calls, and their main purpose is to transform functions signatures from
 *  void pointers and cfitsio's preprocessor type enums to a more type-safe and
 *  convenient interface using overloads and templates.
 *
 *  This was written as part of implementing the afw/table library.  Someday
 *  the afw/image FITS I/O should be modified to use some of these with the goal
 *  of eliminating a lot of code between the two.
 */
 
#include <climits>
#include <optional>
#include <string>
#include <set>
 
#include <boost/format.hpp>
 
#include "lsst/base.h"
#include "lsst/pex/exceptions.h"
#include "lsst/daf/base.h"
#include "ndarray.h"
#include "lsst/afw/fitsCompression.h"
#include "lsst/afw/fitsDefaults.h"
 
namespace lsst {
namespace afw {
namespace fits {

LSST_EXCEPTION_TYPE(FitsError, lsst::pex::exceptions::IoError, lsst::afw::fits::FitsError)
 

LSST_EXCEPTION_TYPE(FitsTypeError, lsst::afw::fits::FitsError, lsst::afw::fits::FitsTypeError)

class HeaderIterationFunctor {
public:
    virtual void operator()(std::string const& key, std::string const& value, std::string const& comment) = 0;
 
    virtual ~HeaderIterationFunctor() {}
};

std::string makeErrorMessage(std::string const& fileName = "", int status = 0, std::string const& msg = "");
inline std::string makeErrorMessage(std::string const& fileName, int status, boost::format const& msg) {
    return makeErrorMessage(fileName, status, msg.str());
}

std::string makeErrorMessage(void* fptr, int status = 0, std::string const& msg = "");
inline std::string makeErrorMessage(void* fptr, int status, boost::format const& msg) {
    return makeErrorMessage(fptr, status, msg.str());
}

std::string makeLimitedFitsHeader(lsst::daf::base::PropertySet const& metadata,
                                  std::set<std::string> const& excludeNames = {});

#define LSST_FITS_EXCEPT(type, fitsObj, ...) \
    type(LSST_EXCEPT_HERE, lsst::afw::fits::makeErrorMessage((fitsObj).fptr, (fitsObj).status, __VA_ARGS__))

#define LSST_FITS_CHECK_STATUS(fitsObj, ...)                                              \
    if ((fitsObj).status != 0) {                                                          \
        auto except = LSST_FITS_EXCEPT(lsst::afw::fits::FitsError, fitsObj, __VA_ARGS__); \
        (fitsObj).status = 0;                                                             \
        throw except;                                                                     \
    }

template <typename T>
int getBitPix();

class MemFileManager {
public:
    MemFileManager() : _ptr(nullptr), _len(0), _managed(true) {}

    explicit MemFileManager(std::size_t len) : _ptr(nullptr), _len(0), _managed(true) { reset(len); }

    MemFileManager(void* ptr, std::size_t len) : _ptr(ptr), _len(len), _managed(false) {}

    void reset();

    void reset(std::size_t len);

    void reset(void* ptr, std::size_t len) {
        reset();
        _ptr = ptr;
        _len = len;
        _managed = false;
    }
 
    ~MemFileManager() { reset(); }
 
    // No copying
    MemFileManager(const MemFileManager&) = delete;
    MemFileManager& operator=(const MemFileManager&) = delete;
 
    // No moving
    MemFileManager(MemFileManager&&) = delete;
    MemFileManager& operator=(MemFileManager&&) = delete;

    void* getData() const { return _ptr; }

    std::size_t getLength() const { return _len; }
 
private:
    friend class Fits;
 
    void* _ptr;
    std::size_t _len;
    bool _managed;
};

template <typename T, int N, int C>
ndarray::Array<T const, N, N> const makeContiguousArray(ndarray::Array<T, N, C> const& array) {
    ndarray::Array<T const, N, N> contiguous = ndarray::dynamic_dimension_cast<N>(array);
    if (contiguous.empty()) contiguous = ndarray::copy(array);
    return contiguous;
}

enum class HduType : int { IMAGE = 0, ASCII_TABLE = 1, BIN_TABLE = 2, ANY = -1 };

class Fits {
    void createImageImpl(int bitpix, int nAxis, long const* nAxes);
    template <typename T>
    void writeImageImpl(T const* data, int nElements, std::optional<T> explicit_null = std::nullopt);
    template <typename T>
    void readImageImpl(int nAxis, T* data, long* begin, long* end, long* increment);
    void getImageShapeImpl(int maxDim, long* nAxes);
 
public:
    enum BehaviorFlags {
        AUTO_CLOSE = 0x01,  // Close files when the Fits object goes out of scope if fptr != NULL
        AUTO_CHECK = 0x02   // Call LSST_FITS_CHECK_STATUS after every cfitsio call
    };

    std::string getFileName() const;

    int getHdu();

    void setHdu(int hdu, bool relative = false);

    void setHdu(std::string const& name, HduType hdutype = HduType::ANY, int hduver = 0);

    int countHdus();


    template <typename T>
    void updateKey(std::string const& key, T const& value, std::string const& comment);
    void updateKey(std::string const& key, char const* value, std::string const& comment) {
        updateKey(key, std::string(value), comment);
    }
    template <typename T>
    void updateKey(std::string const& key, T const& value);
    void updateKey(std::string const& key, char const* value) { updateKey(key, std::string(value)); }


    template <typename T>
    void writeKey(std::string const& key, T const& value, std::string const& comment);
    void writeKey(std::string const& key, char const* value, std::string const& comment) {
        writeKey(key, std::string(value), comment);
    }
    template <typename T>
    void writeKey(std::string const& key, T const& value);
    void writeKey(std::string const& key, char const* value) { writeKey(key, std::string(value)); }


    template <typename T>
    void updateColumnKey(std::string const& prefix, int n, T const& value, std::string const& comment);
    void updateColumnKey(std::string const& prefix, int n, char const* value, std::string const& comment) {
        updateColumnKey(prefix, n, std::string(value), comment);
    }
    template <typename T>
    void updateColumnKey(std::string const& prefix, int n, T const& value);
    void updateColumnKey(std::string const& prefix, int n, char const* value) {
        updateColumnKey(prefix, n, std::string(value));
    }



    template <typename T>
    void writeColumnKey(std::string const& prefix, int n, T const& value, std::string const& comment);
    void writeColumnKey(std::string const& prefix, int n, char const* value, std::string const& comment) {
        writeColumnKey(prefix, n, std::string(value), comment);
    }
    template <typename T>
    void writeColumnKey(std::string const& prefix, int n, T const& value);
    void writeColumnKey(std::string const& prefix, int n, char const* value) {
        writeColumnKey(prefix, n, std::string(value));
    }


    void writeMetadata(daf::base::PropertySet const& metadata);

    void readMetadata(daf::base::PropertySet& metadata, bool strip = false);

    template <typename T>
    void readKey(std::string const& key, T& value);

    void forEachKey(HeaderIterationFunctor& functor);

    void createEmpty();


    template <typename PixelT, int N>
    void createImage(ndarray::Vector<ndarray::Size, N> const& shape) {
        ndarray::Vector<long, N> nAxes(shape.reverse());
        createImageImpl(getBitPix<PixelT>(), N, nAxes.elems);
    }
 
    template <int N>
    void createImage(int bitpix, ndarray::Vector<ndarray::Size, N> const& shape) {
        ndarray::Vector<long, N> nAxes(shape.reverse());
        createImageImpl(bitpix, N, nAxes.elems);
    }


    template <typename PixelT>
    void createImage(long x, long y) {
        long naxes[2] = {x, y};
        createImageImpl(getBitPix<PixelT>(), 2, naxes);
    }

    template <typename T, int N, int C>
    void writeImage(ndarray::Array<T const, N, C> const& array) {
        writeImageImpl(makeContiguousArray(array).getData(), array.getNumElements());
    }


    template <typename T>
    void writeImage(
            image::ImageBase<T> const& image, CompressionOptions const * compression = nullptr,
            daf::base::PropertySet const * header = nullptr,
            image::Mask<image::MaskPixel> const * mask = nullptr);

    int getImageDim();

    template <int N>
    ndarray::Vector<ndarray::Size, N> getImageShape() {
        ndarray::Vector<long, N> nAxes(1);
        getImageShapeImpl(N, nAxes.elems);
        ndarray::Vector<ndarray::Size, N> shape;
        for (int i = 0; i < N; ++i) shape[i] = nAxes[N - i - 1];
        return shape;
    }

    template <typename T>
    bool checkImageType();

    std::string getImageDType();

    template <typename T, int N>
    void readImage(ndarray::Array<T, N, N> const& array, ndarray::Vector<int, N> const& offset) {
        ndarray::Vector<long, N> begin(offset.reverse());
        ndarray::Vector<long, N> end(begin);
        end += array.getShape().reverse();
        ndarray::Vector<long, N> increment(1);
        begin += increment;  // first FITS pixel is 1, not 0
        readImageImpl(N, array.getData(), begin.elems, end.elems, increment.elems);
    }

    void createTable();

    template <typename T>
    int addColumn(std::string const& ttype, int size, std::string const& comment);

    template <typename T>
    int addColumn(std::string const& ttype, int size);

    std::size_t addRows(std::size_t nRows);

    std::size_t countRows();

    template <typename T>
    void writeTableArray(std::size_t row, int col, int nElements, T const* value);

    template <typename T>
    void writeTableScalar(std::size_t row, int col, T value) {
        writeTableArray(row, col, 1, &value);
    }

    void writeTableScalar(std::size_t row, int col, std::string const& value);

    template <typename T>
    void readTableArray(std::size_t row, int col, int nElements, T* value);

    template <typename T>
    void readTableScalar(std::size_t row, int col, T& value) {
        readTableArray(row, col, 1, &value);
    }

    void readTableScalar(std::size_t row, int col, std::string& value, bool isVariableLength);

    long getTableArraySize(int col);

    long getTableArraySize(std::size_t row, int col);

    Fits() : fptr(nullptr), status(0), behavior(0) {}

    Fits(std::string const& filename, std::string const& mode, int behavior);

    Fits(MemFileManager& manager, std::string const& mode, int behavior);

    void closeFile();

    bool checkCompressedImagePhu();
 
    ~Fits() {
        if ((fptr) && (behavior & AUTO_CLOSE)) closeFile();
    }
 
    // No copying
    Fits(const Fits&) = delete;
    Fits& operator=(const Fits&) = delete;
 
    // No moving
    Fits(Fits&&) = delete;
    Fits& operator=(Fits&&) = delete;
 
    void* fptr;    // the actual cfitsio fitsfile pointer; void to avoid including fitsio.h here.
    int status;    // the cfitsio status indicator that gets passed to every cfitsio call.
    int behavior;  // bitwise OR of BehaviorFlags
};


std::shared_ptr<daf::base::PropertyList> combineMetadata(
        daf::base::PropertyList const & first,
        daf::base::PropertyList const & second);

std::shared_ptr<daf::base::PropertyList> readMetadata(std::string const& fileName, int hdu = DEFAULT_HDU,
                                                      bool strip = false);

std::shared_ptr<daf::base::PropertyList> readMetadata(std::string const& fileName, std::string const& hduname,
                                                      HduType type = HduType::ANY, int hduver = 0,
                                                      bool strip = false);

std::shared_ptr<daf::base::PropertyList> readMetadata(fits::MemFileManager& manager, int hdu = DEFAULT_HDU,
                                                      bool strip = false);

std::shared_ptr<daf::base::PropertyList> readMetadata(fits::MemFileManager& manager,
                                                      std::string const& hduname, HduType type = HduType::ANY,
                                                      int hduver = 0, bool strip = false);

std::shared_ptr<daf::base::PropertyList> readMetadata(fits::Fits& fitsfile, bool strip = false);
 

class HduMoveGuard {
public:
 
    HduMoveGuard() = delete;
 
    HduMoveGuard(HduMoveGuard const &) = delete;
    HduMoveGuard(HduMoveGuard &&) = delete;
 
    HduMoveGuard & operator=(HduMoveGuard const &) = delete;
    HduMoveGuard & operator=(HduMoveGuard &&) = delete;

    HduMoveGuard(Fits & fits, int hdu, bool relative=false);
 
    ~HduMoveGuard();

    void disable() { _enabled = false; }
 
private:
    Fits & _fits;
    int _oldHdu;
    bool _enabled;
};
 
 
}  // namespace fits
}  // namespace afw
}  // namespace lsst
 
#endif  // !LSST_AFW_fits_h_INCLUDED