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 <string>

#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) throw LSST_FITS_EXCEPT(lsst::afw::fits::FitsError, fitsObj, __VA_ARGS__)

template <typename T>
int getBitPix();

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

    explicit MemFileManager(std::size_t len) : _ptr(0), _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;
}

struct ImageWriteOptions {
    ImageCompressionOptions compression;  
    ImageScalingOptions scaling;          

    template <typename T>
    explicit ImageWriteOptions(image::Image<T> const& image) : compression(image) {}

    template <typename T>
    explicit ImageWriteOptions(image::Mask<T> const& mask) : compression(mask) {}

    explicit ImageWriteOptions(ImageCompressionOptions const& compression_ =
                                       ImageCompressionOptions(ImageCompressionOptions::NONE),
                               ImageScalingOptions const& scaling_ = ImageScalingOptions())
            : compression(compression_), scaling(scaling_) {}

    explicit ImageWriteOptions(ImageScalingOptions const& scaling_)
            : compression(ImageCompressionOptions::NONE), scaling(scaling_) {}

    ImageWriteOptions(daf::base::PropertySet const& config);

    static std::shared_ptr<daf::base::PropertySet> validate(daf::base::PropertySet const& config);
};

class Fits {
    void createImageImpl(int bitpix, int nAxis, long const* nAxes);
    template <typename T>
    void writeImageImpl(T const* data, int nElements);
    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);

    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(detail::Bitpix<PixelT>::value, 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(detail::Bitpix<PixelT>::value, 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, ImageWriteOptions const& options,
            std::shared_ptr<daf::base::PropertySet const> header = std::shared_ptr<daf::base::PropertyList>(),
            std::shared_ptr<image::Mask<image::MaskPixel> const> mask =
                    std::shared_ptr<image::Mask<image::MaskPixel>>());

    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(0), 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();

    void setImageCompression(ImageCompressionOptions const& options);

    ImageCompressionOptions getImageCompression();

    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(
        std::shared_ptr<const daf::base::PropertyList> first,
        std::shared_ptr<const daf::base::PropertyList> 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(fits::MemFileManager& manager, int hdu = DEFAULT_HDU,
                                                      bool strip = false);
std::shared_ptr<daf::base::PropertyList> readMetadata(fits::Fits& fitsfile, bool strip = false);

void setAllowImageCompression(bool allow);
bool getAllowImageCompression();



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