fitsCompression.h

Go to the documentation of this file.

// -*- lsst-c++ -*-
#ifndef LSST_AFW_fitsCompression_h_INCLUDED
#define LSST_AFW_fitsCompression_h_INCLUDED
 
#include <string>
#include <limits>
 
#include "boost/cstdfloat.hpp"
 
#include "lsst/pex/exceptions.h"
#include "lsst/daf/base.h"
#include "ndarray.h"
#include "ndarray/eigen.h"
 
#include "lsst/afw/image/Image.h"
#include "lsst/afw/image/Mask.h"
 
namespace lsst {
namespace afw {
namespace fits {
 
// Forward declarations
class Fits;
 
namespace detail {

template <typename T>
struct Bitpix;
 
template <>
struct Bitpix<std::uint8_t> {
    static int const value = 8;
};
template <>
struct Bitpix<std::int16_t> {
    static int const value = 16;
};
template <>
struct Bitpix<std::int32_t> {
    static int const value = 32;
};
template <>
struct Bitpix<std::int64_t> {
    static int const value = 64;
};
template <>
struct Bitpix<std::uint16_t> {
    static int const value = 16;
};
template <>
struct Bitpix<std::uint32_t> {
    static int const value = 32;
};
template <>
struct Bitpix<std::uint64_t> {
    static int const value = 64;
};
template <>
struct Bitpix<float> {
    static int const value = -32;
};
template <>
struct Bitpix<double> {
    static int const value = -64;
};

class PixelArrayBase {
public:
    virtual ~PixelArrayBase() {}

    virtual void const* getData() const = 0;

    std::size_t getNumElements() const { return _num; }
 
protected:
    PixelArrayBase(std::size_t num) : _num(num) {}
 
private:
    std::size_t _num;  // Number of pixel values
};

template <typename T>
class PixelArray : public PixelArrayBase {
public:
    PixelArray() = delete;
    PixelArray(PixelArray const&) = delete;

    PixelArray(ndarray::Array<T, 1, 1> const& array)
            : PixelArrayBase(array.getNumElements()),
              _pixels(array.getData()),
              _manager(array.getManager()) {}

    template <typename U>
    PixelArray(ndarray::Array<U, 1, 1> const& array) : PixelArrayBase(array.getNumElements()) {
        auto mem = ndarray::SimpleManager<U>::allocate(getNumElements());
        _manager = mem.first;
        _pixels = mem.second;
        std::copy(array.begin(), array.end(),
                  const_cast<typename std::remove_const<T>::type*>(reinterpret_cast<T const*>(_pixels)));
    }
 
    ~PixelArray() override {}
 
    void const* getData() const override { return _pixels; }
 
private:
    void const* _pixels;             // The data
    ndarray::Manager::Ptr _manager;  // Memory manager; holds onto the data while we use it
};

template <typename T>
std::shared_ptr<PixelArrayBase> makePixelArray(int bitpix, ndarray::Array<T, 1, 1> const& array) {
    switch (bitpix) {
        case 0:
            return std::make_shared<PixelArray<T>>(array);
        case 8:
            return std::make_shared<PixelArray<std::uint8_t>>(array);
        case 16:
            return std::make_shared<PixelArray<std::int16_t>>(array);
        case 32:
            return std::make_shared<PixelArray<std::int32_t>>(array);
        case 64:
            return std::make_shared<PixelArray<std::int64_t>>(array);
        case -32:
            return std::make_shared<PixelArray<boost::float32_t>>(array);
        case -64:
            return std::make_shared<PixelArray<boost::float64_t>>(array);
        default:
            std::ostringstream os;
            os << "Unrecognized bitpix: " << bitpix;
            throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, os.str());
    }
}
 
}  // namespace detail

struct ImageCompressionOptions {
    enum CompressionAlgorithm {
        NONE,          
        GZIP,          
        GZIP_SHUFFLE,  
        RICE,          
        PLIO,          
    };
    using Tiles = ndarray::Array<long, 1, 1>;
 
    CompressionAlgorithm algorithm;  
    Tiles tiles;          
    float quantizeLevel;  

    explicit ImageCompressionOptions(CompressionAlgorithm algorithm_, Tiles tiles_,
                                     float quantizeLevel_ = 0.0)
            : algorithm(algorithm_), tiles(ndarray::copy(tiles_)), quantizeLevel(quantizeLevel_) {}
 
    explicit ImageCompressionOptions(CompressionAlgorithm algorithm_, std::vector<long> tiles_,
                                     float quantizeLevel_ = 0.0)
            : algorithm(algorithm_), tiles(ndarray::allocate(tiles_.size())), quantizeLevel(quantizeLevel_) {
        std::copy(tiles_.cbegin(), tiles_.cend(), tiles.begin());
    }

    explicit ImageCompressionOptions(CompressionAlgorithm algorithm_, int rows = 1,
                                     float quantizeLevel_ = 0.0);

    template <typename T>
    explicit ImageCompressionOptions(image::Image<T> const& image)
            : ImageCompressionOptions(image.getBBox().getArea() > 0 ? NONE : NONE) {}
    template <typename T>
    explicit ImageCompressionOptions(image::Mask<T> const& mask)
            : ImageCompressionOptions(mask.getBBox().getArea() > 0 ? NONE : NONE) {}
 
    // Disable compression for int64: cfitsio won't compress them
    explicit ImageCompressionOptions(image::Image<std::int64_t> const& image)
            : ImageCompressionOptions(NONE) {}
    explicit ImageCompressionOptions(image::Mask<std::int64_t> const& mask) : ImageCompressionOptions(NONE) {}
    explicit ImageCompressionOptions(image::Image<std::uint64_t> const& image)
            : ImageCompressionOptions(NONE) {}
    explicit ImageCompressionOptions(image::Mask<std::uint64_t> const& mask)
            : ImageCompressionOptions(NONE) {}
};

ImageCompressionOptions::CompressionAlgorithm compressionAlgorithmFromString(std::string const& name);

std::string compressionAlgorithmToString(ImageCompressionOptions::CompressionAlgorithm algorithm);

ImageCompressionOptions::CompressionAlgorithm compressionAlgorithmFromCfitsio(int cfitsio);

int compressionAlgorithmToCfitsio(ImageCompressionOptions::CompressionAlgorithm algorithm);

struct ImageScale {
    int bitpix;     
    double bscale;  
    double bzero;   
    long blank;     

    ImageScale(int bitpix_, double bscale_, double bzero_)
            : bitpix(bitpix_),
              bscale(bscale_),
              bzero(std::floor(bzero_ / bscale_ + 0.5) * bscale_),
              blank(bitpix > 0 ? (bitpix == 8 ? 255 : (1L << (bitpix - 1)) - 1) : 0) {}

    template <typename T>
    std::shared_ptr<detail::PixelArrayBase> toFits(
            ndarray::Array<T const, 2, 2> const& image, bool forceNonfiniteRemoval, bool fuzz = true,
            ndarray::Array<long, 1> const& tiles = ndarray::Array<long, 1, 1>(), int seed = 1) const;

    template <typename T>
    ndarray::Array<T, 2, 2> fromFits(ndarray::Array<T, 2, 2> const& image) const;
};

class ImageScalingOptions {
public:
    enum ScalingAlgorithm {
        NONE,            
        RANGE,           
        STDEV_POSITIVE,  
        STDEV_NEGATIVE,  
        STDEV_BOTH,      
        MANUAL,          
    };
    ScalingAlgorithm algorithm;           
    int bitpix;                           
    bool fuzz;                            
    int seed;                             
    std::vector<std::string> maskPlanes;  
    float quantizeLevel;                  
    float quantizePad;  
    double bscale;      
    double bzero;       

    explicit ImageScalingOptions()
            : ImageScalingOptions(NONE, 0, {}, 1, 4.0, 5.0, false, std::numeric_limits<double>::quiet_NaN(),
                                  std::numeric_limits<double>::quiet_NaN()) {}

    ImageScalingOptions(ScalingAlgorithm algorithm_, int bitpix_,
                        std::vector<std::string> const& maskPlanes_ = {}, int seed_ = 1,
                        float quantizeLevel_ = 4.0, float quantizePad_ = 5.0, bool fuzz_ = true,
                        double bscale_ = 1.0, double bzero_ = 0.0);

    ImageScalingOptions(int bitpix_, double bscale_ = 1.0, double bzero_ = 0.0)
            : ImageScalingOptions(MANUAL, bitpix_, {}, 1, 4.0, 5.0, false, bscale_, bzero_) {}


    template <typename T>
    ImageScale determine(image::ImageBase<T> const& image,
                         image::Mask<image::MaskPixel> const * mask = nullptr) const {
        auto const arrays = _toArray(image, mask);
        return determine(arrays.first, arrays.second);
    }
 
    template <typename T, int N>
    ImageScale determine(ndarray::Array<T const, N, N> const& image,
                         ndarray::Array<bool, N, N> const& mask) const;
 
private:
    template <typename T>
    std::pair<ndarray::Array<T const, 2, 2>, ndarray::Array<bool, 2, 2>> _toArray(
            image::ImageBase<T> const& image,
            image::Mask<image::MaskPixel> const * mask = nullptr) const {
        if (mask && image.getDimensions() != mask->getDimensions()) {
            std::ostringstream os;
            os << "Size mismatch between image and mask: ";
            os << image.getWidth() << "x" << image.getHeight();
            os << " vs ";
            os << mask->getWidth() << "x" << mask->getHeight();
            throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, os.str());
        }
        ndarray::Array<T const, 2, 2> imageArray = ndarray::dynamic_dimension_cast<2>(image.getArray());
        if (imageArray.empty()) imageArray = ndarray::copy(image.getArray());
        ndarray::Array<bool, 2, 2> maskArray = ndarray::allocate(imageArray.getShape());
        if (mask) {
            maskArray.deep() = (mask->getArray() & mask->getPlaneBitMask(maskPlanes));
        } else {
            maskArray.deep() = false;
        }
        return std::make_pair(imageArray, maskArray);
    }

    template <typename T, int N>
    ImageScale determineFromRange(ndarray::Array<T const, N, N> const& image,
                                  ndarray::Array<bool, N, N> const& mask, bool isUnsigned = false,
                                  bool cfitsioPadding = true) const;

    template <typename T, int N>
    ImageScale determineFromStdev(ndarray::Array<T const, N, N> const& image,
                                  ndarray::Array<bool, N, N> const& mask, bool isUnsigned = false,
                                  bool cfitsioPadding = true) const;
};

ImageScalingOptions::ScalingAlgorithm scalingAlgorithmFromString(std::string const& name);

std::string scalingAlgorithmToString(ImageScalingOptions::ScalingAlgorithm algorithm);
 
}  // namespace fits
}  // namespace afw
}  // namespace lsst
 
#endif  // ifndef LSST_AFW_fitsCompression_h_INCLUDED