Pixel.h

Go to the documentation of this file.

/*
 * LSST Data Management System
 * Copyright 2008-2016  AURA/LSST.
 *
 * 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/>.
 */

#if !defined(LSST_AFW_IMAGE_PIXEL_H)
#define LSST_AFW_IMAGE_PIXEL_H

#include <cmath>
#include <functional>
#include <iostream>
#include <type_traits>

namespace lsst { namespace afw { namespace image { namespace pixel {

template <typename, typename, typename, typename, typename> class BinaryExpr;

template <typename> struct exprTraits;

template <typename> struct bitwise_or;
template <typename> struct variance_divides;
template <typename> struct variance_multiplies;
template <typename> struct variance_plus;

/************************************************************************************************************/
template<typename _ImagePixelT, typename _MaskPixelT, typename _VariancePixelT=double>
class SinglePixel : public detail::MaskedImagePixel_tag {
public:
    template<typename, typename, typename> friend class Pixel;
    template<typename T> friend class PixelTypeTraits;

    typedef _ImagePixelT ImagePixelT;
    typedef _MaskPixelT MaskPixelT;
    typedef _VariancePixelT VariancePixelT;

    SinglePixel(double const image, int mask=0, double const variance=0) :
        _image(image), _mask(mask), _variance(variance) {}
    SinglePixel(int const image, int mask=0, double const variance=0) :
        _image(image), _mask(mask), _variance(variance) {}

    template<typename rhsExpr>
    SinglePixel(rhsExpr const& rhs) : _image(rhs.image()), _mask(rhs.mask()), _variance(rhs.variance()) {}

    ImagePixelT image() const { return _image; }
    MaskPixelT mask() const { return _mask; }
    VariancePixelT variance() const { return _variance; }

private:
    SinglePixel() :
        _image(std::numeric_limits<_ImagePixelT>::has_quiet_NaN ?
               std::numeric_limits<_ImagePixelT>::quiet_NaN() : 0),
        _mask(0),
        _variance(std::numeric_limits<_VariancePixelT>::has_quiet_NaN ?
                  std::numeric_limits<_VariancePixelT>::quiet_NaN() : 0)
        {}

    ImagePixelT _image;
    MaskPixelT _mask;
    VariancePixelT _variance;
};

template<typename PixelT>
struct PixelTypeTraits
{
    static inline const PixelT padValue() {
        return
            std::numeric_limits<PixelT>::has_quiet_NaN ?
            std::numeric_limits<PixelT>::quiet_NaN()
            : 0;
    }
};

template<typename _ImagePixelT, typename _MaskPixelT, typename _VariancePixelT>
struct PixelTypeTraits<SinglePixel<_ImagePixelT, _MaskPixelT, _VariancePixelT> >
{
    typedef SinglePixel<_ImagePixelT, _MaskPixelT, _VariancePixelT> PixelT;

    static inline const PixelT padValue() {
        return PixelT();
    }
};

template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
SinglePixel<ImagePixelT, MaskPixelT, VariancePixelT> makeSinglePixel(ImagePixelT x, MaskPixelT m, VariancePixelT v) {
    return SinglePixel<ImagePixelT, MaskPixelT, VariancePixelT>(x, m, v);
}

template<typename _ImagePixelT, typename _MaskPixelT, typename _VariancePixelT=double>
class Pixel : public detail::MaskedImagePixel_tag {
public:
    typedef _ImagePixelT ImagePixelT;
    typedef _MaskPixelT MaskPixelT;
    typedef _VariancePixelT VariancePixelT;

#if 0
    Pixel(ImagePixelT& image, MaskPixelT& mask, VariancePixelT& variance) :
        _image(image), _mask(mask), _variance(variance) {}
#else
    //
    // This constructor casts away const.  This should be fixed by making const Pixels.
    //
    Pixel(ImagePixelT const& image, MaskPixelT const& mask=0x0, VariancePixelT const& variance=0) :
        _image(const_cast<ImagePixelT&>(image)),
        _mask(const_cast<MaskPixelT&>(mask)),
        _variance(const_cast<VariancePixelT&>(variance)) {
    }
#endif

    Pixel(SinglePixel<ImagePixelT, MaskPixelT, VariancePixelT>& rhs) :
        _image(rhs._image), _mask(rhs._mask), _variance(rhs._variance) {}

    Pixel operator=(Pixel const& rhs) { // the following template won't stop the compiler trying to generate operator=
        _variance = rhs.variance();     // evaluate before we update image()
        _image = rhs.image();
        _mask = rhs.mask();

        return *this;
    }
    template<typename rhsExpr>
    Pixel operator=(rhsExpr const& rhs) {
        _variance = rhs.variance();     // evaluate before we update image()
        _image = rhs.image();
        _mask = rhs.mask();

        return *this;
    }

    Pixel operator=(double const& rhs_image) {
        _image = rhs_image;
        _mask = 0;
        _variance = 0;

        return *this;
    }

    Pixel operator=(int const& rhs_image) {
        _image = rhs_image;
        _mask = 0;
        _variance = 0;

        return *this;
    }
    ImagePixelT image() const { return _image; }
    MaskPixelT mask() const { return _mask; }
    VariancePixelT variance() const { return _variance; }
    //
    // Logical operators.  We don't need to construct BinaryExpr for them
    // as efficiency isn't a concern.
    //
    template<typename T1>
    friend bool operator==(Pixel const& lhs, T1 const& rhs) {
        return lhs.image() == rhs.image() && lhs.mask() == rhs.mask() && lhs.variance() == rhs.variance();
    }

    template<typename T1>
    friend bool operator!=(Pixel const& lhs, T1 const& rhs) {
        return !(lhs == rhs);
    }

    //
    // Provide friend versions of the op= operators to permit argument promotion on their first arguments
    //
    template<typename ExprT>
    friend Pixel operator+=(Pixel const& e1, ExprT const& e2) {
        Pixel tmp(e1);                  // n.b. shares storage with e1 but gets around "const" (which is required)
        tmp = BinaryExpr<Pixel, ExprT,
            std::plus<ImagePixelT>, bitwise_or<MaskPixelT>, variance_plus<VariancePixelT> >(tmp, e2);
        return tmp;
    }

    template<typename ExprT>
    friend Pixel operator-=(Pixel const& e1, ExprT const& e2) {
        Pixel tmp(e1);                  // n.b. shares storage with e1 but gets around "const" (which is required)
        tmp = BinaryExpr<Pixel, ExprT,
            std::minus<ImagePixelT>, bitwise_or<MaskPixelT>, variance_plus<VariancePixelT> >(tmp, e2);
        return tmp;
    }

    template<typename ExprT>
    friend Pixel operator*=(Pixel const& e1, ExprT const& e2) {
        Pixel tmp(e1);                  // n.b. shares storage with e1 but gets around "const" (which is required)
        tmp = BinaryExpr<Pixel, ExprT,
            std::multiplies<ImagePixelT>, bitwise_or<MaskPixelT>, variance_multiplies<VariancePixelT> >(tmp, e2);
        return tmp;
    }

    template<typename ExprT>
    friend Pixel operator/=(Pixel const& e1, ExprT const& e2) {
        Pixel tmp(e1);                  // n.b. shares storage with e1 but gets around "const" (which is required)
        tmp = BinaryExpr<Pixel, ExprT,
            std::divides<ImagePixelT>, bitwise_or<MaskPixelT>, variance_divides<VariancePixelT> >(tmp, e2);
        return tmp;
    }

private:
    ImagePixelT& _image;
    MaskPixelT& _mask;
    VariancePixelT& _variance;
};

/************************************************************************************************************/
template <typename ExprT>
struct exprTraits {
    typedef ExprT expr_type;
    typedef typename ExprT::ImagePixelT ImagePixelT;
    typedef typename ExprT::MaskPixelT MaskPixelT;
    typedef typename ExprT::VariancePixelT VariancePixelT;
};

template <>
struct exprTraits<double> {
    typedef double ImagePixelT;
    typedef int MaskPixelT;
    typedef double VariancePixelT;
    typedef SinglePixel<ImagePixelT, MaskPixelT, VariancePixelT> expr_type;
};

template <>
struct exprTraits<float> {
    typedef float ImagePixelT;
    typedef exprTraits<double>::MaskPixelT MaskPixelT;
    typedef exprTraits<double>::VariancePixelT VariancePixelT;
    typedef SinglePixel<ImagePixelT, MaskPixelT, VariancePixelT> expr_type;
};

template <>
struct exprTraits<int> {
    typedef int ImagePixelT;
    typedef exprTraits<double>::MaskPixelT MaskPixelT;
    typedef exprTraits<double>::VariancePixelT VariancePixelT;
    typedef SinglePixel<ImagePixelT, MaskPixelT, VariancePixelT> expr_type;
};

template <>
struct exprTraits<unsigned short> {
    typedef int ImagePixelT;
    typedef exprTraits<double>::MaskPixelT MaskPixelT;
    typedef exprTraits<double>::VariancePixelT VariancePixelT;
    typedef SinglePixel<ImagePixelT, MaskPixelT, VariancePixelT> expr_type;
};

/************************************************************************************************************/
template <typename T1>
struct noop : public std::unary_function<T1, T1> {
    T1 operator()(const T1& x) const {
        return x;
    }
};

template <typename T1>
struct bitwise_or : public std::binary_function<T1, T1, T1> {
    T1 operator()(const T1& x, const T1& y) const {
        return (x | y);
    }
    T1 operator()(const T1& x) const {
        return x;
    }
};

template <typename T1>
struct variance_divides {
    T1 operator()(T1 const& x, T1 const& y, T1 const& vx, T1 const& vy) const {
        T1 const x2 = x*x;
        T1 const y2 = y*y;
        T1 const iy2 = 1.0/y2;
        return x2*vy*iy2*iy2 + vx*iy2;
    }

    T1 operator()(T1 const&, T1 const& y, T1 const& vx) const {
        return vx/(y*y);
    }
};
template <typename T1>
struct variance_multiplies {
    T1 operator()(T1 const& x, T1 const& y, T1 const& vx, T1 const& vy) const {
        T1 const x2 = x*x;
        T1 const y2 = y*y;
        return x2*vy + y2*vx;
    }

    T1 operator()(T1 const&, T1 const& y, T1 const& vx) const {
        return vx*y*y;
    }
};
template <typename T1>
struct variance_plus {
    T1 operator()(T1 const&, T1 const&, T1 const& vx, T1 const& vy) const {
        return vx + vy;
    }

    T1 operator()(T1 const&, T1 const&, T1 const& vx) const {
        return vx;
    }
};
template <typename T1>
struct variance_plus_covar {
    variance_plus_covar(double alpha=0) : _alpha(alpha) {}

    T1 operator()(T1 const&, T1 const&, T1 const& vx, T1 const& vy) const {
        return vx + vy + 2*_alpha*sqrt(vx*vy);
    }
    T1 operator()(T1 const&, T1 const&, T1 const& vx) const {
        return vx;
    }
private:
    double _alpha;
};

/************************************************************************************************************/
template <typename ExprT1, typename ImageBinOp, typename MaskBinOp, typename VarianceBinOp>
class UnaryExpr {
public:
    typedef typename exprTraits<ExprT1>::ImagePixelT ImagePixelT;
    typedef typename exprTraits<ExprT1>::MaskPixelT MaskPixelT;
    typedef typename exprTraits<ExprT1>::VariancePixelT VariancePixelT;
    UnaryExpr(ExprT1 e1,
              ImageBinOp imageOp=ImageBinOp(), MaskBinOp maskOp=MaskBinOp(), VarianceBinOp varOp=VarianceBinOp()) :
        _expr1(e1), _imageOp(imageOp), _maskOp(maskOp), _varOp(varOp) {}

    ImagePixelT image() const {
        return _imageOp(_expr1.image());
    }

    MaskPixelT mask() const {
        return _maskOp(_expr1.mask());
    }

    VariancePixelT variance() const {
        return _varOp(_expr1.variance());
    }
private:
    typename exprTraits<ExprT1>::expr_type _expr1;
    ImageBinOp _imageOp;
    MaskBinOp _maskOp;
    VarianceBinOp _varOp;
};

template <typename ExprT1, typename ExprT2, typename ImageBinOp, typename MaskBinOp, typename VarianceBinOp>
class BinaryExpr {
public:
    typedef typename exprTraits<ExprT1>::ImagePixelT ImagePixelT;
    typedef typename exprTraits<ExprT1>::MaskPixelT MaskPixelT;
    typedef typename exprTraits<ExprT1>::VariancePixelT VariancePixelT;
    BinaryExpr(ExprT1 e1, ExprT2 e2,
               ImageBinOp imageOp=ImageBinOp(), MaskBinOp maskOp=MaskBinOp(), VarianceBinOp varOp=VarianceBinOp()) :
        _expr1(e1), _expr2(e2), _imageOp(imageOp), _maskOp(maskOp), _varOp(varOp) {}

    BinaryExpr(ExprT1 e1, ExprT2 e2, double const alpha,
               ImageBinOp imageOp=ImageBinOp(), MaskBinOp maskOp=MaskBinOp(), VarianceBinOp =VarianceBinOp()) :
        _expr1(e1), _expr2(e2), _imageOp(imageOp), _maskOp(maskOp), _varOp(VarianceBinOp(alpha)) {}
    ImagePixelT image() const {
        return _imageOp(_expr1.image(), _expr2.image());
    }

    MaskPixelT mask() const {
        return _maskOp(_expr1.mask(), _expr2.mask());
    }

    VariancePixelT variance() const {
        return _varOp(_expr1.image(), _expr2.image(), _expr1.variance(), _expr2.variance());
    }
private:
    typename exprTraits<ExprT1>::expr_type _expr1;
    typename exprTraits<ExprT2>::expr_type _expr2;
    ImageBinOp _imageOp;
    MaskBinOp _maskOp;
    VarianceBinOp _varOp;
};

template <typename ExprT1, typename ImageBinOp, typename MaskBinOp, typename VarianceBinOp>
class BinaryExpr<ExprT1, double, ImageBinOp, MaskBinOp, VarianceBinOp> {
public:
    typedef typename exprTraits<ExprT1>::ImagePixelT ImagePixelT;
    typedef typename exprTraits<ExprT1>::MaskPixelT MaskPixelT;
    typedef typename exprTraits<ExprT1>::VariancePixelT VariancePixelT;
    BinaryExpr(ExprT1 e1, double e2,
               ImageBinOp imageOp=ImageBinOp(), MaskBinOp maskOp=MaskBinOp(), VarianceBinOp varOp=VarianceBinOp()) :
        _expr1(e1), _expr2(e2), _imageOp(imageOp), _maskOp(maskOp), _varOp(varOp) {}

    BinaryExpr(ExprT1 e1, double e2, double const alpha,
               ImageBinOp imageOp=ImageBinOp(), MaskBinOp maskOp=MaskBinOp(), VarianceBinOp=VarianceBinOp()) :
        _expr1(e1), _expr2(e2), _imageOp(imageOp), _maskOp(maskOp), _varOp(VarianceBinOp(alpha)) {}
    ImagePixelT image() const {
        return _imageOp(_expr1.image(), _expr2);
    }

    MaskPixelT mask() const {
        return _maskOp(_expr1.mask());
    }

    VariancePixelT variance() const {
        return _varOp(_expr1.image(), _expr2, _expr1.variance());
    }
private:
    typename exprTraits<ExprT1>::expr_type _expr1;
    double _expr2;
    ImageBinOp _imageOp;
    MaskBinOp _maskOp;
    VarianceBinOp _varOp;
};

/************************************************************************************************************/
template <typename ExprT1>
UnaryExpr<ExprT1,
          std::negate<typename exprTraits<ExprT1>::ImagePixelT>,
          noop<typename exprTraits<ExprT1>::MaskPixelT>,
          noop<typename exprTraits<ExprT1>::VariancePixelT> > operator-(ExprT1 e1) {
    return UnaryExpr<ExprT1,
                  std::negate<typename exprTraits<ExprT1>::ImagePixelT>,
                  noop<typename exprTraits<ExprT1>::MaskPixelT>,
                  noop<typename exprTraits<ExprT1>::VariancePixelT> >(e1);
}

//------------------------------------------
template <typename ExprT1,typename ExprT2>
BinaryExpr<ExprT1, ExprT2,
           std::plus<typename exprTraits<ExprT1>::ImagePixelT>,
           bitwise_or<typename exprTraits<ExprT1>::MaskPixelT>,
           variance_plus<typename exprTraits<ExprT1>::VariancePixelT> > operator+(ExprT1 e1, ExprT2 e2) {
    return BinaryExpr<ExprT1, ExprT2,
        std::plus<typename exprTraits<ExprT1>::ImagePixelT>,
        bitwise_or<typename exprTraits<ExprT1>::MaskPixelT>,
        variance_plus<typename exprTraits<ExprT1>::VariancePixelT> >(e1,e2);
}

template <typename ExprT1,typename ExprT2>
ExprT1 operator+=(ExprT1& e1, ExprT2 e2) {
    e1 = BinaryExpr<ExprT1, ExprT2,
        std::plus<typename exprTraits<ExprT1>::ImagePixelT>,
        bitwise_or<typename exprTraits<ExprT1>::MaskPixelT>,
        variance_plus<typename exprTraits<ExprT1>::VariancePixelT> >(e1,e2);
    return e1;
}

//
// Implementations of add that work for arithmetic or MaskedImage pixels
//
// The choice is made on the basis of std::is_arithmetic
namespace {
    template <typename ExprT1,typename ExprT2>
    ExprT1 doPlus(ExprT1 e1, ExprT2 e2,
                  double const,
                  boost::mpl::true_) {
        return e1 + e2;
    }

    template <typename ExprT1,typename ExprT2>
    BinaryExpr<ExprT1, ExprT2,
               std::plus<typename exprTraits<ExprT1>::ImagePixelT>,
               bitwise_or<typename exprTraits<ExprT1>::MaskPixelT>,
               variance_plus_covar<typename exprTraits<ExprT1>::VariancePixelT> > doPlus(ExprT1 e1, ExprT2 e2,
                                                                                         double const alpha,
                                                                                        boost::mpl::false_) {
        return BinaryExpr<ExprT1, ExprT2,
            std::plus<typename exprTraits<ExprT1>::ImagePixelT>,
            bitwise_or<typename exprTraits<ExprT1>::MaskPixelT>,
            variance_plus_covar<typename exprTraits<ExprT1>::VariancePixelT> >(e1, e2, alpha);
    }
}

template<typename ExprT1, typename ExprT2>
inline ExprT1 plus(ExprT1& lhs,          
                   ExprT2 const& rhs,    
                   float covariance      
                 ) {
    return doPlus(lhs, rhs, covariance, typename std::is_arithmetic<ExprT1>::type());
}

//------------------------------------------
template <typename ExprT1,typename ExprT2>
BinaryExpr<ExprT1, ExprT2,
           std::minus<typename exprTraits<ExprT1>::ImagePixelT>,
           bitwise_or<typename exprTraits<ExprT1>::MaskPixelT>,
           variance_plus<typename exprTraits<ExprT1>::VariancePixelT> > operator-(ExprT1 e1, ExprT2 e2) {
    return BinaryExpr<ExprT1, ExprT2,
        std::minus<typename exprTraits<ExprT1>::ImagePixelT>,
        bitwise_or<typename exprTraits<ExprT1>::MaskPixelT>,
        variance_plus<typename exprTraits<ExprT1>::VariancePixelT> >(e1,e2);
}

template <typename ExprT1,typename ExprT2>
ExprT1 operator-=(ExprT1& e1, ExprT2 e2) {
    e1 = BinaryExpr<ExprT1, ExprT2,
        std::minus<typename exprTraits<ExprT1>::ImagePixelT>,
        bitwise_or<typename exprTraits<ExprT1>::MaskPixelT>,
        variance_plus<typename exprTraits<ExprT1>::VariancePixelT> >(e1,e2);
    return e1;
}

//------------------------------------------
template <typename ExprT1,typename ExprT2>
BinaryExpr<ExprT1, ExprT2,
           std::multiplies<typename exprTraits<ExprT1>::ImagePixelT>,
           bitwise_or<typename exprTraits<ExprT1>::MaskPixelT>,
           variance_multiplies<typename exprTraits<ExprT1>::VariancePixelT> > operator*(ExprT1 e1, ExprT2 e2) {
    return BinaryExpr<ExprT1, ExprT2,
        std::multiplies<typename exprTraits<ExprT1>::ImagePixelT>,
        bitwise_or<typename exprTraits<ExprT1>::MaskPixelT>,
        variance_multiplies<typename exprTraits<ExprT1>::VariancePixelT> >(e1,e2);
}

template <typename ExprT1,typename ExprT2>
ExprT1 operator*=(ExprT1& e1, ExprT2 e2) {
    e1 = BinaryExpr<ExprT1, ExprT2,
        std::multiplies<typename exprTraits<ExprT1>::ImagePixelT>,
        bitwise_or<typename exprTraits<ExprT1>::MaskPixelT>,
        variance_multiplies<typename exprTraits<ExprT1>::VariancePixelT> >(e1,e2);
    return e1;
}

//------------------------------------------
template <typename ExprT1,typename ExprT2>
BinaryExpr<ExprT1, ExprT2,
           std::divides<typename exprTraits<ExprT1>::ImagePixelT>,
           bitwise_or<typename exprTraits<ExprT1>::MaskPixelT>,
           variance_divides<typename exprTraits<ExprT1>::VariancePixelT> > operator/(ExprT1 e1, ExprT2 e2) {
    return BinaryExpr<ExprT1, ExprT2,
        std::divides<typename exprTraits<ExprT1>::ImagePixelT>,
        bitwise_or<typename exprTraits<ExprT1>::MaskPixelT>,
        variance_divides<typename exprTraits<ExprT1>::VariancePixelT> >(e1,e2);
}

template <typename ExprT1,typename ExprT2>
ExprT1 operator/=(ExprT1& e1, ExprT2 e2) {
    e1 = BinaryExpr<ExprT1, ExprT2,
        std::divides<typename exprTraits<ExprT1>::ImagePixelT>,
        bitwise_or<typename exprTraits<ExprT1>::MaskPixelT>,
        variance_divides<typename exprTraits<ExprT1>::VariancePixelT> >(e1,e2);
    return e1;
}

/************************************************************************************************************/
template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
std::ostream& operator<<(std::ostream &os, SinglePixel<ImagePixelT, MaskPixelT, VariancePixelT> const& v) {
    return os << "(" << v.image() << ", " << v.mask() << ", " << v.variance() << ")";
}

template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
std::ostream& operator<<(std::ostream &os, Pixel<ImagePixelT, MaskPixelT, VariancePixelT> const& v) {
    return os << "(" << v.image() << ", " << v.mask() << ", " << v.variance() << ")";
}

template <typename ExprT1,typename ExprT2, typename BinOp, typename MaskBinOp, typename VarBinOp>
std::ostream& operator<<(std::ostream &os, BinaryExpr<ExprT1, ExprT2, BinOp, MaskBinOp, VarBinOp> const& v) {
    return os << "(" << v.image() << ", " << v.mask() << ", " << v.variance() << ")";
}

}}}}
#endif