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/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:
    using ImagePtr = std::shared_ptr<image::Image<ImagePixelT>>;
    using MaskPtr = std::shared_ptr<image::Mask<MaskPixelT>>;
    using VariancePtr = std::shared_ptr<image::Image<VariancePixelT>>;
    using MaskPlaneDict = typename Mask<MaskPixelT>::MaskPlaneDict;
 
    using Variance = lsst::afw::image::Image<VariancePixelT>;
    using Image = lsst::afw::image::Image<ImagePixelT>;
    using Mask = lsst::afw::image::Mask<MaskPixelT>;
 
    using image_category = detail::MaskedImage_tag;
 
    template <typename ImagePT = ImagePixelT, typename MaskPT = MaskPixelT, typename VarPT = VariancePixelT>
    struct ImageTypeFactory {
        using type = MaskedImage<ImagePT, MaskPT, VarPT>;
    };
 
    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;
 
    using Pixel = lsst::afw::image::pixel::Pixel<ImagePixelT, MaskPixelT, VariancePixelT>;
    using SinglePixel = lsst::afw::image::pixel::SinglePixel<ImagePixelT, MaskPixelT, VariancePixelT>;
 
    template <typename ImageIterator, typename MaskIterator, typename VarianceIterator,
              template <typename> class Ref = Reference>
    class MaskedImageIteratorBase {
        using IMV_iterator_tuple = boost::tuple<ImageIterator, MaskIterator, VarianceIterator>;
 
    public:
        using IMV_tuple = typename boost::zip_iterator<IMV_iterator_tuple>::reference;
        template <typename, typename, typename>
        friend class const_MaskedImageIterator;
        using type = Pixel;
 
        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> {
        using MaskedImageIteratorBase_t = MaskedImageIteratorBase<ImageIterator, MaskIterator, VarianceIterator>;
 
    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> {
        using const_ImageIterator = typename detail::const_iterator_type<ImageIterator>::type;
        using const_MaskIterator = typename detail::const_iterator_type<MaskIterator>::type;
        using const_VarianceIterator = typename detail::const_iterator_type<VarianceIterator>::type;
 
        using MaskedImageIteratorBase_t = MaskedImageIteratorBase<const_ImageIterator, const_MaskIterator, const_VarianceIterator, ConstReference>;
 
    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 {
        using IMVLocator = typename boost::tuple<ImageLocator, MaskLocator, VarianceLocator>;
        //
        // 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 {
            using IterT = typename LocT::x_iterator;
 
        public:
            apply_x(LocT& loc) : _loc(loc) {}
            IterT& operator()() { return _loc.x(); }
 
        private:
            LocT& _loc;
        };
 
        template <typename LocT>
        class apply_y {
            using IterT = typename LocT::y_iterator;
 
        public:
            apply_y(LocT& loc) : _loc(loc) {}
            IterT& operator()() { return _loc.y(); }
 
        private:
            LocT& _loc;
        };
 
    public:
        using IMVCachedLocation = typename boost::tuple<typename ImageLocator::cached_location_t, typename MaskLocator::cached_location_t, typename VarianceLocator::cached_location_t>;
        using x_iterator = _x_or_y_iterator<apply_x>;
        using y_iterator = _x_or_y_iterator<apply_y>;
        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
        //
        using PixelTVec = typename boost::mpl::vector<ImagePixelT, MaskPixelT, VariancePixelT>;
 
        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> {
        using MaskedImageLocatorBase_t = MaskedImageLocatorBase<ImageLocator, MaskLocator, VarianceLocator>;
 
    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> {
        using const_ImageLocator = typename detail::const_locator_type<ImageLocator>::type;
        using const_MaskLocator = typename detail::const_locator_type<MaskLocator>::type;
        using const_VarianceLocator = typename detail::const_locator_type<VarianceLocator>::type;
 
        using MaskedImageLocatorBase_t = MaskedImageLocatorBase<const_ImageLocator, const_MaskLocator, const_VarianceLocator, ConstReference>;
 
    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
    using iterator = MaskedImageIterator<typename Image::iterator, typename Mask::iterator, typename Variance::iterator>;
    // A const_iterator to a MaskedImage
    using const_iterator = const_MaskedImageIterator<typename Image::iterator, typename Mask::iterator, typename Variance::iterator>;
    // A reverse_iterator to a MaskedImage
    using reverse_iterator = MaskedImageIterator<typename Image::reverse_iterator, typename Mask::reverse_iterator, typename Variance::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
    using x_iterator = MaskedImageIterator<typename Image::x_iterator, typename Mask::x_iterator, typename Variance::x_iterator>;
    using const_x_iterator = const_MaskedImageIterator<typename Image::x_iterator, typename Mask::x_iterator, typename Variance::x_iterator>;
    using fast_iterator = x_iterator;
    using y_iterator = MaskedImageIterator<typename Image::y_iterator, typename Mask::y_iterator, typename Variance::y_iterator>;
    using const_y_iterator = const_MaskedImageIterator<typename Image::y_iterator, typename Mask::y_iterator, typename Variance::y_iterator>;
 
    using xy_locator = MaskedImageLocator<typename Image::xy_locator, typename Mask::xy_locator, typename Variance::xy_locator>;
    using const_xy_locator = const_MaskedImageLocator<typename Image::xy_locator, typename Mask::xy_locator, typename Variance::xy_locator>;
 
    using xy_x_iterator = typename MaskedImageLocator<typename Image::xy_locator, typename Mask::xy_locator, typename Variance::xy_locator>::x_iterator;
    using xy_y_iterator = typename MaskedImageLocator<typename Image::xy_locator, typename Mask::xy_locator, typename Variance::xy_locator>::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  
                )
            : _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); }
 
    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_ {
    using type = ImageT;
    static std::shared_ptr<type> getImage(std::shared_ptr<ImageT> image) { return image; }
};
 
template <typename ImageT>
struct GetImage_<ImageT, typename image::detail::MaskedImage_tag> {
    using type = typename ImageT::Image;
    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