Image.h

Go to the documentation of this file.

// -*- lsst-c++ -*-

/* 
 * LSST Data Management System
 * Copyright 2008, 2009, 2010 LSST Corporation.
 * 
 * This product includes software developed by the
 * LSST Project (http://www.lsst.org/).
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the LSST License Statement and 
 * the GNU General Public License along with this program.  If not, 
 * see <http://www.lsstcorp.org/LegalNotices/>.
 */
 
#ifndef LSST_AFW_IMAGE_IMAGE_H
#define LSST_AFW_IMAGE_IMAGE_H

#include <list>
#include <map>
#include <string>
#include <utility>
#include <functional>

#include "boost/mpl/bool.hpp"
#include "boost/shared_ptr.hpp"

#include "lsst/afw/geom/Extent.h"
#include "lsst/afw/geom/Box.h"
#include "lsst/afw/geom/Point.h"
#include "lsst/afw/image/lsstGil.h"
#include "lsst/afw/image/Utils.h"
#include "lsst/afw/image/ImageUtils.h"
#include "lsst/afw/math/Function.h"
#include "lsst/daf/base.h"
#include "lsst/daf/base/Citizen.h"
#include "lsst/pex/exceptions.h"
#include "ndarray.h"

namespace lsst { namespace afw {

namespace formatters {
    template <typename PixelT> class ImageFormatter;
    template <typename PixelT> class DecoratedImageFormatter;
}

namespace fits {
class Fits;
class MemFileManager;
}

namespace image {
    namespace detail {
        //
        // Traits for image types
        //
        struct basic_tag { };
        struct Image_tag : public basic_tag { };
        template<typename ImageT>
        struct image_traits {
            typedef typename ImageT::image_category image_category;
        };
        //
        std::string const wcsNameForXY0 = "A"; // the name of the WCS to use to save (X0, Y0) to FITS files; e.g. "A"
    }

    /*********************************************************************************************************/
    class CheckIndices {
    public:
        explicit CheckIndices(bool check=true) : _check(check) {}
        operator bool() const { return _check; }
    private:
        bool _check;
    };

    /*********************************************************************************************************/
    template<typename PixelT>
    struct Reference {
        typedef typename boost::gil::channel_traits<PixelT>::reference type; 
    };
    template<typename PixelT>
    struct ConstReference {
        typedef typename boost::gil::channel_traits<PixelT>::const_reference type; 
    };
    
    enum ImageOrigin {PARENT, LOCAL};

    //
    // You are not expected to use this class directly in your own code; use one of the
    // specialised subclasses
    //
    template<typename PixelT>
    class ImageBase : public lsst::daf::base::Persistable,
                      public lsst::daf::base::Citizen {
    private:
        typedef typename lsst::afw::image::detail::types_traits<PixelT>::view_t _view_t;
        typedef typename lsst::afw::image::detail::types_traits<PixelT>::const_view_t _const_view_t;


        typedef ndarray::Manager Manager;
    public:        

        typedef boost::shared_ptr<ImageBase<PixelT> > Ptr; 
        typedef boost::shared_ptr<const ImageBase<PixelT> > ConstPtr; 

        typedef detail::basic_tag image_category; 

        typedef PixelT SinglePixel;
        typedef PixelT Pixel;
        typedef typename Reference<PixelT>::type PixelReference;
        typedef typename ConstReference<PixelT>::type PixelConstReference;
        typedef typename _view_t::xy_locator xy_locator;
        typedef typename _view_t::xy_locator::const_t const_xy_locator;
        typedef typename _view_t::iterator iterator;
        typedef typename _const_view_t::iterator const_iterator;
        typedef typename _view_t::reverse_iterator reverse_iterator;
        typedef typename _const_view_t::reverse_iterator const_reverse_iterator;
        typedef typename _view_t::x_iterator x_iterator;
        typedef x_iterator fast_iterator;
        typedef typename _view_t::x_iterator xy_x_iterator;
        typedef typename _const_view_t::x_iterator const_x_iterator;
        typedef typename _view_t::y_iterator y_iterator;
        typedef typename _view_t::y_iterator xy_y_iterator;
        typedef typename _const_view_t::y_iterator const_y_iterator;
        typedef typename ndarray::Array<PixelT, 2, 1> Array;
        typedef typename ndarray::Array<PixelT const, 2, 1> ConstArray;

        template<typename OtherPixelT> friend class ImageBase; // needed by generalised copy constructors
        //
        //
        template<typename SinglePixelT>
        static SinglePixel PixelCast(SinglePixelT rhs) {
            return SinglePixel(rhs);
        }
        //
        // DecoratedImage needs enough access to ImageBase to read data from disk; we might be able to design around this
        //
        template<typename> friend class DecoratedImage;
        template<typename, typename, typename> friend class MaskedImage;
        explicit ImageBase(const geom::Extent2I  & dimensions=geom::Extent2I());
        explicit ImageBase(const geom::Box2I &bbox);
        ImageBase(const ImageBase& src, const bool deep=false);
        explicit ImageBase(const ImageBase& src, const geom::Box2I& bbox,
                           const ImageOrigin origin=PARENT, const bool deep=false);
        template<typename OtherPixelT>
        ImageBase(const ImageBase<OtherPixelT>& rhs, const bool deep) :
            lsst::daf::base::Citizen(typeid(this)) {
            if (!deep) {
                throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
                    "Only deep copies are permitted for ImageBases with different pixel types");
            }

            ImageBase<PixelT> tmp(rhs.getBBox());
            copy_and_convert_pixels(rhs._gilView, tmp._gilView); // from boost::gil

            using std::swap;                           // See Meyers, Effective C++, Item 25
            ImageBase<PixelT>::swap(tmp);                  // See Meyers, Effective C++, Items 11 and 43
        }

        explicit ImageBase(
            Array const & array, bool deep = false, geom::Point2I const & xy0 = geom::Point2I()
        );

        virtual ~ImageBase() { }
        ImageBase& operator=(const ImageBase& rhs);
        ImageBase& operator=(const PixelT rhs);
        void operator<<=(const ImageBase& rhs);
        //
        // Operators etc.
        //
        PixelReference operator()(int x, int y);
        PixelReference operator()(int x, int y, CheckIndices const&);
        PixelConstReference operator()(int x, int y) const;
        PixelConstReference operator()(int x, int y, CheckIndices const&) const;

        PixelConstReference get0(int x, int y) const {
            return operator()(x-getX0(), y-getY0());
        }
        PixelConstReference get0(int x, int y, CheckIndices const& check) const {
            return operator()(x-getX0(), y-getY0(), check);
        }
        void set0(int x, int y, const PixelT v) {
            operator()(x-getX0(), y-getY0()) = v;
        }
        void set0(int x, int y, const PixelT v, CheckIndices const& check) {
            operator()(x-getX0(), y-getY0(), check) = v;
        }

        int getWidth() const { return _gilView.width(); }
        int getHeight() const { return _gilView.height(); }
        int getX0() const { return _origin.getX(); }
        int getY0() const { return _origin.getY(); }

        geom::Point2I getXY0() const { return _origin; }
        
        std::pair<int, double> positionToIndex(
                double const pos, 
                lsst::afw::image::xOrY const xy 
        ) const {
            double const fullIndex = pos - PixelZeroPos - (xy == X ? getX0() : getY0());
            int const roundedIndex = static_cast<int>(fullIndex + 0.5);
            double const residual = fullIndex - roundedIndex;
            return std::pair<int, double>(roundedIndex, residual);
        }

        inline double indexToPosition(
                double ind, 
                lsst::afw::image::xOrY const xy 
        ) const {
            return ind + PixelZeroPos + (xy == X ? getX0() : getY0());
        }
        
        geom::Extent2I getDimensions() const { return geom::Extent2I(getWidth(), getHeight()); }
        
        void swap(ImageBase &rhs);

        Array getArray();
        ConstArray getArray() const;
        //
        // Iterators and Locators
        //
        iterator begin() const;
        iterator end() const;
        reverse_iterator rbegin() const;
        reverse_iterator rend() const;
        iterator at(int x, int y) const;

        fast_iterator begin(bool) const;
        fast_iterator end(bool) const;

        x_iterator row_begin(int y) const {
            return _gilView.row_begin(y);
        }

        x_iterator row_end(int y) const {
            return _gilView.row_end(y);
        }

        x_iterator x_at(int x, int y) const { return _gilView.x_at(x, y); }

        y_iterator col_begin(int x) const {
            return _gilView.col_begin(x);
        }
        
        y_iterator col_end(int x) const {
            return _gilView.col_end(x);
        }

        y_iterator y_at(int x, int y) const {
            return _gilView.y_at(x, y);
        }

        xy_locator xy_at(int x, int y) const {
            return xy_locator(_gilView.xy_at(x, y));
        }
        void setXY0(geom::Point2I const origin) {
            _origin=origin;
        }
        void setXY0(int const x0, int const y0) {
            setXY0(geom::Point2I(x0,y0));
        }

        geom::Box2I getBBox(ImageOrigin origin=PARENT) const {
            if (origin == PARENT) {
                return geom::Box2I(_origin, getDimensions());
            }
            else return geom::Box2I(geom::Point2I(0,0), getDimensions());
        }
    private:
        geom::Point2I _origin;
        Manager::Ptr _manager;
        _view_t _gilView;

        //oring of ImageBase in some larger image as returned to and manipulated
        //by the user


    protected:
#if !defined(SWIG)
        static _view_t _allocateView(geom::Extent2I const & dimensions, Manager::Ptr & manager);
        static _view_t _makeSubView(
            geom::Extent2I const & dimensions, 
            geom::Extent2I const & offset, 
            const _view_t & view
        );

        _view_t _getRawView() const { return _gilView; }

#endif
        inline bool isContiguous() const {
            return begin()+getWidth()*getHeight() == end();
        }    
    };

    template<typename PixelT>
    void swap(ImageBase<PixelT>& a, ImageBase<PixelT>& b);

    /************************************************************************************************************/
    template<typename PixelT>
    class Image : public ImageBase<PixelT> {
    public:
        template<typename, typename, typename> friend class MaskedImage;
        typedef boost::shared_ptr<Image<PixelT> > Ptr;
        typedef boost::shared_ptr<const Image<PixelT> > ConstPtr;

        typedef detail::Image_tag image_category;

#if !defined(SWIG)
        template<typename ImagePT=PixelT>
        struct ImageTypeFactory {
            typedef Image<ImagePT> type;
        };
#endif
        template<typename OtherPixelT> friend class Image; // needed by generalised copy constructors
        
        explicit Image(unsigned int width, unsigned int height, PixelT initialValue=0);
        explicit Image(geom::Extent2I const & dimensions=geom::Extent2I(), PixelT initialValue=0);
        explicit Image(geom::Box2I const & bbox, PixelT initialValue=0);

        explicit Image(Image const & rhs, geom::Box2I const & bbox, ImageOrigin const origin=PARENT, 
                       const bool deep=false);
        Image(const Image& rhs, const bool deep=false);

        explicit Image(
            std::string const & fileName, int hdu=0,
            PTR(lsst::daf::base::PropertySet) metadata=PTR(lsst::daf::base::PropertySet)(),
            geom::Box2I const & bbox=geom::Box2I(), 
            ImageOrigin origin=PARENT
        );

        explicit Image(
            fits::MemFileManager & manager, int hdu=0,
            PTR(lsst::daf::base::PropertySet) metadata=PTR(lsst::daf::base::PropertySet)(),
            geom::Box2I const & bbox=geom::Box2I(), 
            ImageOrigin origin=PARENT
        );

        explicit Image(
            fits::Fits & fitsfile,
            PTR(lsst::daf::base::PropertySet) metadata=PTR(lsst::daf::base::PropertySet)(),
            geom::Box2I const & bbox=geom::Box2I(), 
            ImageOrigin origin=PARENT
        );

        // generalised copy constructor
        template<typename OtherPixelT>
        Image(Image<OtherPixelT> const& rhs, const bool deep) :
            image::ImageBase<PixelT>(rhs, deep) {}

        explicit Image(ndarray::Array<PixelT,2,1> const & array, bool deep = false,
                       geom::Point2I const & xy0 = geom::Point2I()) :
            image::ImageBase<PixelT>(array, deep, xy0) {}

        virtual ~Image() { }
        //
        // Assignment operators are not inherited
        //
        Image& operator=(const PixelT rhs);
        Image& operator=(const Image& rhs);

        void writeFits(
            std::string const& fileName,
            CONST_PTR(lsst::daf::base::PropertySet) metadata = CONST_PTR(lsst::daf::base::PropertySet)(),
            std::string const& mode="w"
        ) const;

        void writeFits(
            fits::MemFileManager & manager,
            CONST_PTR(lsst::daf::base::PropertySet) metadata = CONST_PTR(lsst::daf::base::PropertySet)(),
            std::string const& mode="w"
        ) const;

        void writeFits(
            fits::Fits & fitsfile,
            CONST_PTR(lsst::daf::base::PropertySet) metadata = CONST_PTR(lsst::daf::base::PropertySet)()
        ) const;

        static Image readFits(std::string const & filename, int hdu=0) {
            return Image<PixelT>(filename, hdu);
        }

        static Image readFits(fits::MemFileManager & manager, int hdu=0) {
            return Image<PixelT>(manager, hdu);
        }

        void swap(Image &rhs);
        //
        // Operators etc.
        //
        void operator+=(PixelT const rhs);
        virtual void operator+=(Image<PixelT>const & rhs);
        void operator+=(lsst::afw::math::Function2<double> const& function);
        void scaledPlus(double const c, Image<PixelT>const & rhs);
        void operator-=(PixelT const rhs);
        void operator-=(Image<PixelT> const& rhs);
        void operator-=(lsst::afw::math::Function2<double> const& function);
        void scaledMinus(double const c, Image<PixelT>const & rhs);
        void operator*=(PixelT const rhs);
        void operator*=(Image<PixelT> const& rhs);
        void scaledMultiplies(double const c, Image<PixelT>const & rhs);
        void operator/=(PixelT const rhs);
        void operator/=(Image<PixelT> const& rhs);
        void scaledDivides(double const c, Image<PixelT>const & rhs);

        // In-place per-pixel sqrt().  Useful when handling variance planes.
        void sqrt();
    protected:
        using ImageBase<PixelT>::_getRawView;
    private:
        LSST_PERSIST_FORMATTER(lsst::afw::formatters::ImageFormatter<PixelT>)
    };
    
    template<typename LhsPixelT, typename RhsPixelT>
    void operator+=(Image<LhsPixelT> &lhs, Image<RhsPixelT> const& rhs);
    template<typename LhsPixelT, typename RhsPixelT>
    void operator-=(Image<LhsPixelT> &lhs, Image<RhsPixelT> const& rhs);
    template<typename LhsPixelT, typename RhsPixelT>
    void operator*=(Image<LhsPixelT> &lhs, Image<RhsPixelT> const& rhs);
    template<typename LhsPixelT, typename RhsPixelT>
    void operator/=(Image<LhsPixelT> &lhs, Image<RhsPixelT> const& rhs);

    template<typename PixelT>
    void swap(Image<PixelT>& a, Image<PixelT>& b);

/************************************************************************************************************/
    template<typename PixelT>
    class DecoratedImage : public lsst::daf::base::Persistable,
                           public lsst::daf::base::Citizen {
    public:
        typedef boost::shared_ptr<DecoratedImage> Ptr;
        typedef boost::shared_ptr<const DecoratedImage> ConstPtr;
        typedef PTR(Image<PixelT>) ImagePtr;
        typedef CONST_PTR(Image<PixelT>) ImageConstPtr;

        explicit DecoratedImage(const geom::Extent2I & dimensions=geom::Extent2I());
        explicit DecoratedImage(const geom::Box2I & bbox);
        explicit DecoratedImage(PTR(Image<PixelT>) rhs);
        DecoratedImage(DecoratedImage const& rhs, const bool deep=false);
        explicit DecoratedImage(
            std::string const& fileName, 
            const int hdu=0, 
            geom::Box2I const& bbox=geom::Box2I(), 
            ImageOrigin const origin = PARENT
        );

        DecoratedImage& operator=(const DecoratedImage& image);

        PTR(lsst::daf::base::PropertySet) getMetadata() const { return _metadata; }
        void setMetadata(PTR(lsst::daf::base::PropertySet) metadata) { _metadata = metadata; }

        int getWidth() const { return _image->getWidth(); }
        int getHeight() const { return _image->getHeight(); }
        
        int getX0() const { return _image->getX0(); }
        int getY0() const { return _image->getY0(); }

        const geom::Extent2I getDimensions() const { return _image->getDimensions(); }

        void swap(DecoratedImage &rhs);
        
        void writeFits(
            std::string const& fileName,
            CONST_PTR(lsst::daf::base::PropertySet) metadata = CONST_PTR(lsst::daf::base::PropertySet)(),
            std::string const& mode="w"
        ) const;

        ImagePtr      getImage()       { return _image; }
        ImageConstPtr getImage() const { return _image; }

        double getGain() const { return _gain; }
        void setGain(double gain) { _gain = gain; }
    private:
        LSST_PERSIST_FORMATTER(lsst::afw::formatters::DecoratedImageFormatter<PixelT>)
        PTR(Image<PixelT>) _image;
        PTR(lsst::daf::base::PropertySet) _metadata;
        
        double _gain;

        void init();
    };

    template<typename PixelT>
    void swap(DecoratedImage<PixelT>& a, DecoratedImage<PixelT>& b);

}}}  // lsst::afw::image

#endif