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/>.
 */

/*
 * Implementation of the Class MaskedImage
 */

#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/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)$";
}  // namespace detail
}  // namespace image
}  // namespace afw
}  // namespace lsst

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

namespace lsst {
namespace afw {
namespace image {

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

    typedef lsst::afw::image::Image<VariancePixelT> Variance;
    typedef lsst::afw::image::Image<ImagePixelT> Image;
    typedef lsst::afw::image::Mask<MaskPixelT> Mask;

    typedef detail::MaskedImage_tag image_category;

    template <typename ImagePT = ImagePixelT, typename MaskPT = MaskPixelT, typename VarPT = VariancePixelT>
    struct ImageTypeFactory {
        typedef MaskedImage<ImagePT, MaskPT, VarPT> type;
    };

    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;

    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(); }

        MaskedImageIteratorBase& operator+=(std::ptrdiff_t delta  
        ) {
            _iter += delta;
            return *this;
        }
        MaskedImageIteratorBase& operator-=(std::ptrdiff_t delta  
        ) {
            _iter -= delta;
            return *this;
        }
        MaskedImageIteratorBase& operator++() {  // prefix
            ++_iter;
            return *this;
        }
        MaskedImageIteratorBase operator++(int) {  // postfix
            MaskedImageIteratorBase tmp(*this);
            _iter++;
            return tmp;
        }
        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) {}

            _x_or_y_iterator& 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;
                return *this;
            }

            _x_or_y_iterator& 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>())();
                return *this;
            }

            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>())) {
            ;
        }
    };

    // 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(lsst::geom::Extent2I const& dimensions = lsst::geom::Extent2I(),
                         MaskPlaneDict const& planeDict = MaskPlaneDict());
    explicit MaskedImage(ImagePtr image, MaskPtr mask = MaskPtr(), VariancePtr variance = VariancePtr());
    explicit MaskedImage(lsst::geom::Box2I const& bbox, MaskPlaneDict const& planeDict = MaskPlaneDict());

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

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

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

    MaskedImage(MaskedImage const& rhs, bool const deep = false);
    MaskedImage(MaskedImage&& rhs);

    MaskedImage(MaskedImage const& rhs, lsst::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 = ImagePtr(new Image(*rhs.getImage(), deep));
        _mask = MaskPtr(new Mask(*rhs.getMask(), deep));
        _variance = VariancePtr(new Variance(*rhs.getVariance(), deep));
    }

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

    virtual ~MaskedImage() = default;

    void swap(MaskedImage& rhs);

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

    MaskedImage subset(lsst::geom::Box2I const & bbox, ImageOrigin origin=PARENT) const {
        return MaskedImage(*this, bbox, origin, false);
    }

    MaskedImage operator[](lsst::geom::Box2I const & bbox) const {
        return subset(bbox);
    }

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

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

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

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

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

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

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

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

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

    void writeFits(std::string const& fileName, fits::ImageWriteOptions const& imageOptions,
                   fits::ImageWriteOptions const& maskOptions, fits::ImageWriteOptions const& varianceOptions,
                   std::shared_ptr<daf::base::PropertySet const> metadata = nullptr,
                   std::shared_ptr<daf::base::PropertySet const> imageMetadata = nullptr,
                   std::shared_ptr<daf::base::PropertySet const> maskMetadata = nullptr,
                   std::shared_ptr<daf::base::PropertySet const> varianceMetadata = nullptr) const;

    void writeFits(fits::MemFileManager& manager, fits::ImageWriteOptions const& imageOptions,
                   fits::ImageWriteOptions const& maskOptions, fits::ImageWriteOptions const& varianceOptions,
                   std::shared_ptr<daf::base::PropertySet const> metadata = nullptr,
                   std::shared_ptr<daf::base::PropertySet const> imageMetadata = nullptr,
                   std::shared_ptr<daf::base::PropertySet const> maskMetadata = nullptr,
                   std::shared_ptr<daf::base::PropertySet const> varianceMetadata = nullptr) const;

    void writeFits(fits::Fits& fitsfile, fits::ImageWriteOptions const& imageOptions,
                   fits::ImageWriteOptions const& maskOptions, fits::ImageWriteOptions const& varianceOptions,
                   std::shared_ptr<daf::base::PropertySet const> metadata = nullptr,
                   std::shared_ptr<daf::base::PropertySet const> imageMetadata = nullptr,
                   std::shared_ptr<daf::base::PropertySet const> maskMetadata = nullptr,
                   std::shared_ptr<daf::base::PropertySet const> varianceMetadata = nullptr) 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() const { return _image; }

    void setImage(Image const& other) { _image->assign(other); }

    MaskPtr getMask() const { return _mask; }

    void setMask(Mask const& other) { _mask->assign(other); }

    void setVariance(Variance const& other) { _variance->assign(other); }

    VariancePtr getVariance() const { return _variance; }

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

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

    void setXY0(lsst::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 contiguous) const;
    fast_iterator end(bool contiguous) 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:
    void conformSizes();

    ImagePtr _image;
    MaskPtr _mask;
    VariancePtr _variance;
};

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

template <typename ImagePixelT1, typename ImagePixelT2>
bool imagesOverlap(MaskedImage<ImagePixelT1, MaskPixel, VariancePixel> const& image1,
                   MaskedImage<ImagePixelT2, MaskPixel, VariancePixel> const& image2);

/*
 * Some metafunctions to extract an Image pointer from a MaskedImage pointer (or return the original Image
 * pointer)
 *
 * 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(shared_ptr<ImageT> image) {
 *    typename shared_ptr<GetImage<ImageT>::type> im = GetImage<ImageT>::getImage(image);
 * }
 *
 * "im" is a shared_ptr<Image> irrespective of whether ImageT is Masked or not.
 */
namespace {
template <typename ImageT, typename TagT>
struct GetImage_ {
    typedef ImageT type;
    static std::shared_ptr<type> getImage(std::shared_ptr<ImageT> image) { return image; }
};

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

template <typename ImageT>
struct GetImage : public GetImage_<ImageT, typename ImageT::image_category> {};
}  // namespace image
}  // namespace afw
}  // namespace lsst

#endif  //  LSST_IMAGE_MASKEDIMAGE_H