MaskedImage.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_IMAGE_MASKEDIMAGE_H
#define LSST_IMAGE_MASKEDIMAGE_H

#include <list>
#include <map>
#include <memory>
#include <ostream>
#include <string>

#include "boost/mpl/at.hpp"
#include "boost/iterator/zip_iterator.hpp"
#include "boost/tuple/tuple.hpp" // cannot convert to std::tuple (yet) because of use with boost::gil

#include "lsst/daf/base/Citizen.h"
#include "lsst/daf/base/Persistable.h"
#include "lsst/afw/formatters/MaskedImageFormatter.h"
#include "lsst/afw/image/Image.h"
#include "lsst/afw/image/Mask.h"

namespace lsst {
namespace afw {
namespace image {
    namespace detail {
        struct MaskedImage_tag : public basic_tag { };
        struct MaskedImagePixel_tag { };

        std::string const fitsFile_RE = "\\.fits(\\.[fg]z)?$"; 
        std::string const compressedFileNoMEF_RE = "(\\.gz)$";   
    }
}}}

#include "lsst/afw/image/Pixel.h"
#include "lsst/afw/image/LsstImageTypes.h"

namespace lsst {
namespace afw {
    namespace formatters {
        template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT> class MaskedImageFormatter;
    }

namespace image {

template<typename ImagePixelT, typename MaskPixelT=lsst::afw::image::MaskPixel,
         typename VariancePixelT=lsst::afw::image::VariancePixel>
class MaskedImage : public lsst::daf::base::Persistable,
                    public lsst::daf::base::Citizen {
public:
    typedef typename Image<ImagePixelT>::Ptr ImagePtr;
    typedef typename Mask<MaskPixelT>::Ptr MaskPtr;
    typedef typename Image<VariancePixelT>::Ptr VariancePtr;
    typedef std::shared_ptr<MaskedImage> Ptr;
    typedef std::shared_ptr<const MaskedImage> ConstPtr;
    typedef typename Mask<MaskPixelT>::MaskPlaneDict MaskPlaneDict;

    typedef lsst::afw::image::Image<VariancePixelT> Variance; // These need to be here, and in this order, as
    typedef lsst::afw::image::Image<ImagePixelT> Image;       // "typedef Image::Ptr ImagePtr;" confuses swig (it can't
    typedef lsst::afw::image::Mask<MaskPixelT> Mask;          // find ImagePtr) and we can't use Image<> after these typedefs

    typedef detail::MaskedImage_tag image_category;

#if !defined(SWIG)
    template<typename ImagePT=ImagePixelT, typename MaskPT=MaskPixelT, typename VarPT=VariancePixelT>
    struct ImageTypeFactory {
        typedef MaskedImage<ImagePT, MaskPT, VarPT> type;
    };
#endif

    /************************************************************************************************************/
    template<typename, typename, typename> class MaskedImageIterator;
    template<typename, typename, typename> class const_MaskedImageIterator;
    template<typename, typename, typename> class MaskedImageLocator;
    template<typename, typename, typename> class const_MaskedImageLocator;

    /************************************************************************************************************/

#if !defined(SWIG)
    typedef lsst::afw::image::pixel::Pixel<ImagePixelT, MaskPixelT, VariancePixelT> Pixel;
    typedef lsst::afw::image::pixel::SinglePixel<ImagePixelT, MaskPixelT, VariancePixelT> SinglePixel;

    /************************************************************************************************************/
    template<typename ImageIterator, typename MaskIterator, typename VarianceIterator,
             template<typename> class Ref=Reference>
    class MaskedImageIteratorBase {
        typedef boost::tuple<ImageIterator, MaskIterator, VarianceIterator> IMV_iterator_tuple;

    public:
        typedef typename boost::zip_iterator<IMV_iterator_tuple>::reference IMV_tuple;
        template<typename, typename, typename> friend class const_MaskedImageIterator;
        typedef Pixel type;

        MaskedImageIteratorBase(ImageIterator const& img, MaskIterator const& msk, VarianceIterator const &var) :
            _iter(boost::make_zip_iterator(boost::make_tuple(img, msk, var))) {
        }
        typename Ref<ImagePixelT>::type image() {
            return _iter->template get<0>()[0];
        }

        typename Ref<MaskPixelT>::type mask() {
            return _iter->template get<1>()[0];
        }

        typename Ref<VariancePixelT>::type variance() {
            return _iter->template get<2>()[0];
        }

        const IMV_iterator_tuple get_iterator_tuple() const {
            return _iter.get_iterator_tuple();
        }

        void operator+=(std::ptrdiff_t delta 
                       ) {
            _iter += delta;
        }
        void operator-=(std::ptrdiff_t delta 
                       ) {
            _iter -= delta;
        }
        void operator++() {         // prefix
            ++_iter;
        }
        void operator++(int) {      // postfix
            _iter++;
        }
        std::ptrdiff_t operator-(MaskedImageIteratorBase const& rhs) {
            return &this->_iter->template get<0>() - &rhs._iter->template get<0>();
        }
        bool operator==(MaskedImageIteratorBase const& rhs) {
            return &this->_iter->template get<0>() == &rhs._iter->template get<0>();
        }
        bool operator!=(MaskedImageIteratorBase const& rhs) {
            return &this->_iter->template get<0>() != &rhs._iter->template get<0>();
        }
        bool operator<(MaskedImageIteratorBase const& rhs) {
            return &this->_iter->template get<0>() < &rhs._iter->template get<0>();
        }
        operator Pixel() const {
            return Pixel(_iter->template get<0>()[0],
                         _iter->template get<1>()[0],
                         _iter->template get<2>()[0]);
        }

        Pixel operator*() {
            return Pixel(image(), mask(), variance());
        }
        const Pixel operator*() const {
            return Pixel(image(), mask(), variance());
        }

    protected:
        typename boost::zip_iterator<IMV_iterator_tuple> _iter;
    };

    template<typename ImageIterator, typename MaskIterator, typename VarianceIterator>
    class MaskedImageIterator :
        public  MaskedImageIteratorBase<ImageIterator, MaskIterator, VarianceIterator> {
        typedef MaskedImageIteratorBase<ImageIterator, MaskIterator, VarianceIterator> MaskedImageIteratorBase_t;
    public:
        MaskedImageIterator(ImageIterator& img, MaskIterator& msk, VarianceIterator &var) :
            MaskedImageIteratorBase_t(img, msk, var) {
        }
        MaskedImageIterator operator+(std::ptrdiff_t delta) {
            MaskedImageIterator lhs = *this;
            lhs += delta;

            return lhs;
        }
    };

    template<typename ImageIterator, typename MaskIterator, typename VarianceIterator>
    class const_MaskedImageIterator :
        public  MaskedImageIteratorBase<typename detail::const_iterator_type<ImageIterator>::type,
                                        typename detail::const_iterator_type<MaskIterator>::type,
                                        typename detail::const_iterator_type<VarianceIterator>::type,
                                        ConstReference> {

        typedef typename detail::const_iterator_type<ImageIterator>::type const_ImageIterator;
        typedef typename detail::const_iterator_type<MaskIterator>::type const_MaskIterator;
        typedef typename detail::const_iterator_type<VarianceIterator>::type const_VarianceIterator;

        typedef MaskedImageIteratorBase<const_ImageIterator, const_MaskIterator,
                                        const_VarianceIterator, ConstReference> MaskedImageIteratorBase_t;
    public:
        const_MaskedImageIterator(MaskedImageIterator<ImageIterator, MaskIterator, VarianceIterator> const& iter) :
            MaskedImageIteratorBase_t(const_ImageIterator(iter.get_iterator_tuple().template get<0>()),
                                      const_MaskIterator(iter.get_iterator_tuple().template get<1>()),
                                      const_VarianceIterator(iter.get_iterator_tuple().template get<2>())
                                     ) {
            ;
        }
        const_MaskedImageIterator& operator+(std::ptrdiff_t delta) {
            const_MaskedImageIterator lhs = *this;
            lhs += delta;

            return lhs;
        }
    };

    /************************************************************************************************************/
    template<typename ImageLocator, typename MaskLocator, typename VarianceLocator,
             template<typename> class Ref=Reference>
    class MaskedImageLocatorBase {
        typedef typename boost::tuple<ImageLocator, MaskLocator, VarianceLocator> IMVLocator;
        //
        // A class to provide _[xy]_iterator for MaskedImageLocator.  We can't just use
        // a zip_iterator as moving this iterator must be the same as moving the locator
        // itself, for consistency with {Image,Mask}::xy_locator
        //
        template<template<typename> class X_OR_Y >
        class _x_or_y_iterator {
        public:
            _x_or_y_iterator(MaskedImageLocatorBase* mil) : _mil(mil) {}

            void operator+=(const int di) {
                // Equivalent to "_mil->_loc.template get<0>().x() += di;"
                X_OR_Y<ImageLocator>(_mil->_loc.template get<0>())() += di;
                X_OR_Y<MaskLocator>(_mil->_loc.template get<1>())() += di;
                X_OR_Y<VarianceLocator>(_mil->_loc.template get<2>())() += di;
            }

            void operator++() {     // prefix
                // Equivalent to "++_mil->_loc.template get<0>().x();"
                ++X_OR_Y<ImageLocator>(_mil->_loc.template get<0>())();
                ++X_OR_Y<MaskLocator>(_mil->_loc.template get<1>())();
                ++X_OR_Y<VarianceLocator>(_mil->_loc.template get<2>())();
            }

            bool operator==(_x_or_y_iterator const& rhs) {
                return
                    X_OR_Y<ImageLocator>(_mil->_loc.template get<0>())() ==
                    X_OR_Y<ImageLocator>(rhs._mil->_loc.template get<0>())();

            }
            bool operator!=(_x_or_y_iterator const& rhs) {
                return
                    X_OR_Y<ImageLocator>(_mil->_loc.template get<0>())() !=
                    X_OR_Y<ImageLocator>(rhs._mil->_loc.template get<0>())();

            }
            bool operator<(_x_or_y_iterator const& rhs) {
                return
                    X_OR_Y<ImageLocator>(_mil->_loc.template get<0>())() <
                    X_OR_Y<ImageLocator>(rhs._mil->_loc.template get<0>())();

            }

            Pixel operator*() {
                return Pixel((*(X_OR_Y<ImageLocator>(_mil->_loc.template get<0>())()))[0],
                             (*(X_OR_Y<MaskLocator>(_mil->_loc.template get<1>())()))[0],
                             (*(X_OR_Y<VarianceLocator>(_mil->_loc.template get<2>())()))[0]);
            }

            typename Ref<ImagePixelT>::type image() {
                // Equivalent to "return (*_mil->_loc.template get<0>().x())[0];"

                return (*(X_OR_Y<ImageLocator>(_mil->_loc.template get<0>())()))[0];
            }
            typename Ref<MaskPixelT>::type mask() {
                return (*(X_OR_Y<MaskLocator>(_mil->_loc.template get<1>())()))[0];
            }
            typename Ref<VariancePixelT>::type variance() {
                return (*(X_OR_Y<VarianceLocator>(_mil->_loc.template get<2>())()))[0];
            }
        protected:
            MaskedImageLocatorBase *_mil;
        };

        // Two classes to provide .x() and .y() in _x_or_y_iterator
        template<typename LocT>
        class apply_x {
            typedef typename LocT::x_iterator IterT;
        public:
            apply_x(LocT &loc) : _loc(loc) { }
            IterT& operator()() { return _loc.x(); }
        private:
            LocT& _loc;
        };

        template<typename LocT>
        class apply_y {
            typedef typename LocT::y_iterator IterT;
        public:
            apply_y(LocT &loc) : _loc(loc) { }
            IterT& operator()() { return _loc.y(); }
        private:
            LocT& _loc;
        };

    public:
        typedef typename boost::tuple<typename ImageLocator::cached_location_t,
                                      typename MaskLocator::cached_location_t,
                                      typename VarianceLocator::cached_location_t> IMVCachedLocation;
        typedef _x_or_y_iterator<apply_x> x_iterator;
        typedef _x_or_y_iterator<apply_y> y_iterator;
        class cached_location_t {
        public:
            //template<typename, typename, typename, template<typename> class> friend class MaskedImageLocatorBase;
            template<typename, typename, typename> friend class const_MaskedImageLocator;

            cached_location_t(IMVLocator const& loc, int x, int y) :
                _imv(loc.template get<0>().cache_location(x, y),
                     loc.template get<1>().cache_location(x, y),
                     loc.template get<2>().cache_location(x, y)) {
                ;
            }
            //protected:
            IMVCachedLocation _imv;
        };
        MaskedImageLocatorBase(ImageLocator const& img, MaskLocator const& msk, VarianceLocator const& var) :
            _loc(img, msk, var) {
            ;
        }

        Pixel operator*() {
            return Pixel(_loc.template get<0>().x()[0][0],
                         _loc.template get<1>().x()[0][0],
                         _loc.template get<2>().x()[0][0]);
        }

        Pixel operator()(int x, int y) {
            return Pixel(_loc.template get<0>()(x, y)[0],
                         _loc.template get<1>()(x, y)[0],
                         _loc.template get<2>()(x, y)[0]);
        }

        Pixel operator[](cached_location_t const& cached_loc) {
            return Pixel(_loc.template get<0>()[cached_loc._imv.template get<0>()][0],
                         _loc.template get<1>()[cached_loc._imv.template get<1>()][0],
                         _loc.template get<2>()[cached_loc._imv.template get<2>()][0]);
        }
        x_iterator x() {
            return x_iterator(this);
        }

        y_iterator y() {
            return y_iterator(this);
        }

        cached_location_t cache_location(int x, int y) const {
            return cached_location_t(_loc, x, y);
        }
        //
        // We don't want to duplicate code for image/mask/variance -- but the boost::mpl stuff isn't pretty
        // as we can't say int_<N> within a template<int N>.  So define a set of functions apply_IMV
        // to do the dirty work
        //
        typedef typename boost::mpl::vector<ImagePixelT, MaskPixelT, VariancePixelT> PixelTVec;

        template<typename N>
        typename Ref<typename boost::mpl::at<PixelTVec, N>::type>::type apply_IMV(cached_location_t const& cached_loc) {
            return _loc.template get<N::value>()[cached_loc._imv.template get<N::value>()][0];
        }

        template<typename N>
        typename Ref<typename boost::mpl::at<PixelTVec, N>::type>::type apply_IMV() {
            return _loc.template get<N::value>()[0][0];
        }

        template<typename N>
        typename Ref<typename boost::mpl::at<PixelTVec, N>::type>::type apply_IMV(int x, int y) {
            return _loc.template get<N::value>()(x, y)[0];
        }
        //
        // Use those templated classes to implement image/mask/variance
        //
        typename Ref<ImagePixelT>::type image(cached_location_t const& cached_loc) {
            return apply_IMV<boost::mpl::int_<0> >(cached_loc);
        }
        typename Ref<ImagePixelT>::type image() {
            return apply_IMV<boost::mpl::int_<0> >();
        }
        typename Ref<ImagePixelT>::type image(int x, int y) {
            return apply_IMV<boost::mpl::int_<0> >(x, y);
        }

        typename Ref<MaskPixelT>::type mask(cached_location_t const& cached_loc) {
            return apply_IMV<boost::mpl::int_<1> >(cached_loc);
        }
        typename Ref<MaskPixelT>::type mask() {
            return apply_IMV<boost::mpl::int_<1> >();
        }
        typename Ref<MaskPixelT>::type mask(int x, int y) {
            return apply_IMV<boost::mpl::int_<1> >(x, y);
        }

        typename Ref<VariancePixelT>::type variance(cached_location_t const& cached_loc) {
            return apply_IMV<boost::mpl::int_<2> >(cached_loc);
        }
        typename Ref<VariancePixelT>::type variance() {
            return apply_IMV<boost::mpl::int_<2> >();
        }
        typename Ref<VariancePixelT>::type variance(int x, int y) {
            return apply_IMV<boost::mpl::int_<2> >(x, y);
        }

        bool operator==(MaskedImageLocatorBase const& rhs) {
            return _loc.template get<0>() == rhs._loc.template get<0>();
        }
        bool operator!=(MaskedImageLocatorBase const& rhs) {
            return !(*this == rhs);
        }
        bool operator<(MaskedImageLocatorBase const& rhs) {
            return _loc.template get<0>() < rhs._loc.template get<0>();
        }

        MaskedImageLocatorBase& operator+=(pair2I const& p) {
            return operator+=(detail::difference_type(p.first, p.second));
        }

        MaskedImageLocatorBase& operator+=(detail::difference_type p) {
            _loc.template get<0>() += p;
            _loc.template get<1>() += p;
            _loc.template get<2>() += p;

            return *this;
        }

        // Workaround for DM-5590: clang-3.8 cannot access _loc from
        // friend class const_MaskedImageLocator.
        IMVLocator const & getLoc() const { return _loc; }

    protected:
        IMVLocator _loc;
    };

    template<typename ImageLocator, typename MaskLocator, typename VarianceLocator>
    class MaskedImageLocator :
        public  MaskedImageLocatorBase<ImageLocator, MaskLocator, VarianceLocator> {
        typedef MaskedImageLocatorBase<ImageLocator, MaskLocator, VarianceLocator> MaskedImageLocatorBase_t;
    public:
        MaskedImageLocator(ImageLocator& img, MaskLocator& msk, VarianceLocator &var) :
            MaskedImageLocatorBase_t(img, msk, var) {
        }
    };

    template<typename ImageLocator, typename MaskLocator, typename VarianceLocator>
    class const_MaskedImageLocator :
        public  MaskedImageLocatorBase<typename detail::const_locator_type<ImageLocator>::type,
                                       typename detail::const_locator_type<MaskLocator>::type,
                                       typename detail::const_locator_type<VarianceLocator>::type,
                                       ConstReference> {

        typedef typename detail::const_locator_type<ImageLocator>::type const_ImageLocator;
        typedef typename detail::const_locator_type<MaskLocator>::type const_MaskLocator;
        typedef typename detail::const_locator_type<VarianceLocator>::type const_VarianceLocator;

        typedef MaskedImageLocatorBase<const_ImageLocator, const_MaskLocator, const_VarianceLocator,
                                       ConstReference> MaskedImageLocatorBase_t;
    public:
        const_MaskedImageLocator(MaskedImageLocator<ImageLocator, MaskLocator, VarianceLocator> const& iter) :
            MaskedImageLocatorBase_t(const_ImageLocator(iter.getLoc().template get<0>()),
                                     const_MaskLocator(iter.getLoc().template get<1>()),
                                     const_VarianceLocator(iter.getLoc().template get<2>())
                                    ) {
            ;
        }
    };

#endif  // !defined(SWIG)

/************************************************************************************************************/
    // An iterator to a MaskedImage
    typedef MaskedImageIterator<typename Image::iterator,
                                typename Mask::iterator, typename Variance::iterator> iterator;
    // A const_iterator to a MaskedImage
    typedef const_MaskedImageIterator<typename Image::iterator,
                                typename Mask::iterator, typename Variance::iterator> const_iterator;
    // A reverse_iterator to a MaskedImage
    typedef MaskedImageIterator<typename Image::reverse_iterator,
                                typename Mask::reverse_iterator, typename Variance::reverse_iterator> reverse_iterator;
#if 0                                   // doesn't compile.  I should fix this, but it's low priority. RHL
    typedef const_MaskedImageIterator<typename Image::reverse_iterator,
                                typename Mask::reverse_iterator, typename Variance::reverse_iterator> const_reverse_iterator;
#endif
    typedef MaskedImageIterator<typename Image::x_iterator,
                                typename Mask::x_iterator, typename Variance::x_iterator> x_iterator;
    typedef const_MaskedImageIterator<typename Image::x_iterator,
                                typename Mask::x_iterator, typename Variance::x_iterator> const_x_iterator;
    typedef x_iterator fast_iterator;
    typedef MaskedImageIterator<typename Image::y_iterator,
                                typename Mask::y_iterator, typename Variance::y_iterator> y_iterator;
    typedef const_MaskedImageIterator<typename Image::y_iterator,
                                typename Mask::y_iterator, typename Variance::y_iterator> const_y_iterator;

    typedef MaskedImageLocator<typename Image::xy_locator,
                                typename Mask::xy_locator, typename Variance::xy_locator> xy_locator;
    typedef const_MaskedImageLocator<typename Image::xy_locator,
                                typename Mask::xy_locator, typename Variance::xy_locator> const_xy_locator;

    typedef typename MaskedImageLocator<typename Image::xy_locator,
                               typename Mask::xy_locator, typename Variance::xy_locator>::x_iterator xy_x_iterator;
    typedef typename MaskedImageLocator<typename Image::xy_locator,
                               typename Mask::xy_locator, typename Variance::xy_locator>::y_iterator xy_y_iterator;

    /************************************************************************************************************/

    // Constructors
    explicit MaskedImage(
        unsigned int width, unsigned int height,
        MaskPlaneDict const& planeDict=MaskPlaneDict()
    );
    explicit MaskedImage(
        geom::Extent2I const & dimensions=geom::Extent2I(),
        MaskPlaneDict const& planeDict=MaskPlaneDict()
    );
    explicit MaskedImage(
        ImagePtr image,
        MaskPtr mask = MaskPtr(),
        VariancePtr variance = VariancePtr()
    );
    explicit MaskedImage(
        geom::Box2I const & bbox,
        MaskPlaneDict const& planeDict=MaskPlaneDict()
    );

    explicit MaskedImage(
        std::string const & fileName,
        PTR(daf::base::PropertySet) metadata=PTR(daf::base::PropertySet)(),
        geom::Box2I const & bbox=geom::Box2I(), ImageOrigin origin=PARENT,
        bool conformMasks=false, bool needAllHdus=false,
        PTR(daf::base::PropertySet) imageMetadata=PTR(daf::base::PropertySet)(),
        PTR(daf::base::PropertySet) maskMetadata=PTR(daf::base::PropertySet)(),
        PTR(daf::base::PropertySet) varianceMetadata=PTR(daf::base::PropertySet)()
    );

    explicit MaskedImage(
        fits::MemFileManager & manager,
        PTR(daf::base::PropertySet) metadata=PTR(daf::base::PropertySet)(),
        geom::Box2I const & bbox=geom::Box2I(), ImageOrigin origin=PARENT,
        bool conformMasks=false, bool needAllHdus=false,
        PTR(daf::base::PropertySet) imageMetadata=PTR(daf::base::PropertySet)(),
        PTR(daf::base::PropertySet) maskMetadata=PTR(daf::base::PropertySet)(),
        PTR(daf::base::PropertySet) varianceMetadata=PTR(daf::base::PropertySet)()
    );

    explicit MaskedImage(
        fits::Fits & fitsfile,
        PTR(daf::base::PropertySet) metadata=PTR(daf::base::PropertySet)(),
        geom::Box2I const & bbox=geom::Box2I(), ImageOrigin origin=PARENT,
        bool conformMasks=false, bool needAllHdus=false,
        PTR(daf::base::PropertySet) imageMetadata=PTR(daf::base::PropertySet)(),
        PTR(daf::base::PropertySet) maskMetadata=PTR(daf::base::PropertySet)(),
        PTR(daf::base::PropertySet) varianceMetadata=PTR(daf::base::PropertySet)()
    );

    MaskedImage(
        MaskedImage const& rhs,
        bool const deep=false
    );
    MaskedImage(
        MaskedImage const & rhs,
        geom::Box2I const & bbox,
        ImageOrigin const origin=PARENT,
        bool const deep=false
    );
    template<typename OtherPixelT>
    MaskedImage(
        MaskedImage<OtherPixelT, MaskPixelT, VariancePixelT> const& rhs, 

        const bool deep     
    ) :
        lsst::daf::base::Citizen(typeid(this)), _image(), _mask(), _variance() {
        if (!deep) {
            throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
                "Only deep copies are permitted for MaskedImages with different pixel types");
        }

        _image =    typename Image::Ptr(new Image(*rhs.getImage(), deep));
        _mask =     typename Mask::Ptr(new Mask(*rhs.getMask(), deep));
        _variance = typename Variance::Ptr(new Variance(*rhs.getVariance(), deep));
    }

#if defined(DOXYGEN)
    MaskedImage& operator=(MaskedImage const& rhs);
#endif

    virtual ~MaskedImage() {}

    void swap(MaskedImage &rhs);

    // Operators
    MaskedImage& operator=(Pixel const& rhs);
    MaskedImage& operator=(SinglePixel const& rhs);

    void operator<<=(MaskedImage const& rhs);

    void assign(MaskedImage const &rsh, geom::Box2I const &bbox = geom::Box2I(), ImageOrigin origin=PARENT);

    void operator+=(ImagePixelT const rhs);
    void operator+=(MaskedImage const& rhs);
    void operator+=(lsst::afw::image::Image<ImagePixelT> const& rhs) {
        *_image += rhs;
    }
    void operator+=(lsst::afw::math::Function2<double> const& function) {
        *_image += function;
    }
    void scaledPlus(double const c, MaskedImage const& rhs);

    void operator-=(ImagePixelT const rhs);
    void operator-=(MaskedImage const& rhs);
    void operator-=(lsst::afw::image::Image<ImagePixelT> const& rhs) {
        *_image -= rhs;
    }
    void operator-=(lsst::afw::math::Function2<double> const& function) {
        *_image -= function;
    }
    void scaledMinus(double const c, MaskedImage const& rhs);

    void operator*=(ImagePixelT const rhs);
    void operator*=(MaskedImage const& rhs);
    void operator*=(lsst::afw::image::Image<ImagePixelT> const& rhs) {
        *_image *= rhs;
        *_variance *= rhs;           // yes, multiply twice
        *_variance *= rhs;
    }
    void scaledMultiplies(double const c, MaskedImage const& rhs);

    void operator/=(ImagePixelT const rhs);
    void operator/=(MaskedImage const& rhs);
    void operator/=(lsst::afw::image::Image<ImagePixelT> const& rhs) {
        *_image /= rhs;
        *_variance /= rhs; // yes, divide twice
        *_variance /= rhs;
    }
    void scaledDivides(double const c, MaskedImage const& rhs);

    void writeFits(
        std::string const & fileName,
        CONST_PTR(daf::base::PropertySet) metadata = CONST_PTR(daf::base::PropertySet)(),
        CONST_PTR(daf::base::PropertySet) imageMetadata = CONST_PTR(daf::base::PropertySet)(),
        CONST_PTR(daf::base::PropertySet) maskMetadata = CONST_PTR(daf::base::PropertySet)(),
        CONST_PTR(daf::base::PropertySet) varianceMetadata = CONST_PTR(daf::base::PropertySet)()
    ) const;

    void writeFits(
        fits::MemFileManager & manager,
        CONST_PTR(daf::base::PropertySet) metadata = CONST_PTR(daf::base::PropertySet)(),
        CONST_PTR(daf::base::PropertySet) imageMetadata = CONST_PTR(daf::base::PropertySet)(),
        CONST_PTR(daf::base::PropertySet) maskMetadata = CONST_PTR(daf::base::PropertySet)(),
        CONST_PTR(daf::base::PropertySet) varianceMetadata = CONST_PTR(daf::base::PropertySet)()
    ) const;

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

    static MaskedImage readFits(std::string const & filename) {
        return MaskedImage<ImagePixelT, MaskPixelT, VariancePixelT>(filename);
    }

    static MaskedImage readFits(fits::MemFileManager & manager) {
        return MaskedImage<ImagePixelT, MaskPixelT, VariancePixelT>(manager);
    }

    // Getters
    ImagePtr getImage(bool const noThrow=false) const {
        if (!_image && !noThrow) {
            throw LSST_EXCEPT(
                lsst::pex::exceptions::RuntimeError,
                "MaskedImage's Image is NULL"
            );
        }
        return _image;
    }
    MaskPtr getMask(bool const noThrow=false) const {
        if (!_mask && !noThrow) {
            throw LSST_EXCEPT(
                lsst::pex::exceptions::RuntimeError,
                "MaskedImage's Mask is NULL"
            );
        }

        return _mask;
    }
    VariancePtr getVariance(bool const noThrow=false) const {
        if (!_variance && !noThrow) {
            throw LSST_EXCEPT(
                lsst::pex::exceptions::RuntimeError,
                "MaskedImage's Variance is NULL"
            );
        }

        return _variance;
    }
    int getWidth() const { return _image->getWidth(); }
    int getHeight() const { return _image->getHeight(); }
    geom::Extent2I getDimensions() const {return _image->getDimensions();}
    geom::Box2I getBBox(ImageOrigin const origin=PARENT) const {return _image->getBBox(origin);}
    int getX0() const { return _image->getX0(); }
    int getY0() const { return _image->getY0(); }
    geom::Point2I getXY0() const { return _image->getXY0(); }

    void setXY0(int const x0, int const y0) {
        setXY0(geom::Point2I(x0,y0));
    }

    void setXY0(geom::Point2I const origin) {
        if (_image) {
            _image->setXY0(origin);
        }

        if (_mask) {
            _mask->setXY0(origin);
        }

        if (_variance) {
            _variance->setXY0(origin);
        }
    }


    inline double indexToPosition(
            double ind, 
            lsst::afw::image::xOrY const xy 
    ) const {
        return getImage()->indexToPosition(ind, xy);
    }

    std::pair<int, double> positionToIndex(
            double const pos, 
            lsst::afw::image::xOrY const xy 
    ) const {
        return getImage()->positionToIndex(pos, xy);
    }

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

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

    x_iterator row_begin(int y) const;
    x_iterator row_end(int y) const;

    x_iterator x_at(int x, int y) const {
#if 0
        typename Image::x_iterator imageEnd = getImage()->x_at(x, y);
        typename Mask::x_iterator maskEnd = getMask()->x_at(x, y);
        typename Variance::x_iterator varianceEnd = getVariance()->x_at(x, y);
#else  // bypass checks for non-NULL pointers
        typename Image::x_iterator imageEnd = _image->x_at(x, y);
        typename Mask::x_iterator maskEnd = _mask->x_at(x, y);
        typename Variance::x_iterator varianceEnd = _variance->x_at(x, y);
#endif

        return x_iterator(imageEnd, maskEnd, varianceEnd);
    }

    y_iterator col_begin(int x) const;
    y_iterator col_end(int x) const;

    y_iterator y_at(int x, int y) const {
#if 0
        typename Image::y_iterator imageEnd = getImage()->y_at(x, y);
        typename Mask::y_iterator maskEnd = getMask()->y_at(x, y);
        typename Variance::y_iterator varianceEnd = getVariance()->y_at(x, y);
#else  // bypass checks for non-NULL pointers
        typename Image::y_iterator imageEnd = _image->y_at(x, y);
        typename Mask::y_iterator maskEnd = _mask->y_at(x, y);
        typename Variance::y_iterator varianceEnd = _variance->y_at(x, y);
#endif
        return y_iterator(imageEnd, maskEnd, varianceEnd);
    }


    xy_locator xy_at(int x, int y) const {
#if 0
        typename Image::xy_locator imageEnd = getImage()->xy_at(x, y);
        typename Mask::xy_locator maskEnd = getMask()->xy_at(x, y);
        typename Variance::xy_locator varianceEnd = getVariance()->xy_at(x, y);
#else  // bypass checks for non-NULL pointers
        typename Image::xy_locator imageEnd = _image->xy_at(x, y);
        typename Mask::xy_locator maskEnd = _mask->xy_at(x, y);
        typename Variance::xy_locator varianceEnd = _variance->xy_at(x, y);
#endif

        return xy_locator(imageEnd, maskEnd, varianceEnd);
    }

private:

    LSST_PERSIST_FORMATTER(lsst::afw::formatters::MaskedImageFormatter<ImagePixelT, MaskPixelT, VariancePixelT>)
    void conformSizes();

    ImagePtr _image;
    MaskPtr _mask;
    VariancePtr _variance;
};

template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
MaskedImage<ImagePixelT, MaskPixelT, VariancePixelT>* makeMaskedImage(
    typename Image<ImagePixelT>::Ptr image, 
    typename Mask<MaskPixelT>::Ptr mask = typename Mask<MaskPixelT>::Ptr(),    
    typename Image<VariancePixelT>::Ptr variance = typename Image<VariancePixelT>::Ptr() 
                                                                     ) {
    return new MaskedImage<ImagePixelT, MaskPixelT, VariancePixelT>(image, mask, variance);
}

/*
 * Some metafunctions to extract an Image::Ptr from a MaskedImage::Ptr (or return the original Image::Ptr)
 *
 * GetImage is the public interface (it forwards the tag --- just for the sake of the UI); the real work
 * is in GetImage_ which defines a typedef for the Image and a static function, getImage
 *
 * E.g.
 * In the function
 *
 * template<typename ImageT>
 * void func(typename ImageT::Ptr image) {
 *    typename GetImage<ImageT>::type::Ptr im = GetImage<ImageT>::getImage(image);
 * }
 *
 * "im" is an Image::Ptr irrespective of whether ImageT is Masked or not.
 */
namespace {
template<typename ImageT, typename TagT>
struct GetImage_ {
    typedef ImageT type;
    static typename type::Ptr getImage(typename ImageT::Ptr image) {
        return image;
    }
};

template<typename ImageT>
struct GetImage_<ImageT, typename image::detail::MaskedImage_tag> {
    typedef typename ImageT::Image type;
    static typename type::Ptr getImage(typename ImageT::Ptr image) {
        return image->getImage();
    }
};
} // anonymous namespace

template<typename ImageT>
struct GetImage : public GetImage_<ImageT, typename ImageT::image_category> {
};

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

#endif //  LSST_IMAGE_MASKEDIMAGE_H