diff --git a/include/lsst/afw/image/Pixel.h b/include/lsst/afw/image/Pixel.h
index 2ec0e6867..7e46e1f50 100644
--- a/include/lsst/afw/image/Pixel.h
+++ b/include/lsst/afw/image/Pixel.h
@@ -88,7 +88,7 @@ class SinglePixel : public detail::MaskedImagePixel_tag {
     SinglePixel(rhsExpr const& rhs,
                 // ensure this ctor isn't invoked for simple numeric types, which should use
                 // the overload above.
-                typename std::enable_if<!std::is_fundamental<rhsExpr>::value, void*>::type dummy = nullptr)
+                typename std::enable_if_t<!std::is_fundamental_v<rhsExpr>, void*> = nullptr)
             : _image(rhs.image()), _mask(rhs.mask()), _variance(rhs.variance()) {}
 
     ImagePixelT image() const { return _image; }
@@ -123,14 +123,7 @@ struct PixelTypeTraits {
     }
 };
 
-/// Specialization for a %pixel of a MaskedImage
-template <typename _ImagePixelT, typename _MaskPixelT, typename _VariancePixelT>
-struct PixelTypeTraits<SinglePixel<_ImagePixelT, _MaskPixelT, _VariancePixelT> > {
-    using PixelT = SinglePixel<_ImagePixelT, _MaskPixelT, _VariancePixelT>;
 
-    /// The quantity to use when a pixel value is undefined
-    static inline const PixelT padValue() { return PixelT(); }
-};
 
 /** Return a SinglePixel
  *
@@ -152,10 +145,6 @@ class Pixel : public detail::MaskedImagePixel_tag {
     using VariancePixelT = _VariancePixelT;
 
 /// Construct a Pixel from references to its image/mask/variance components
-#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.
     //
@@ -163,7 +152,6 @@ class Pixel : public detail::MaskedImagePixel_tag {
             : _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) {}
@@ -400,26 +388,6 @@ struct variance_plus {
     T1 operator()(T1 const&, T1 const&, T1 const& vx) const { return vx; }
 };
 
-/** The variance of the sum of a pair of correlated pixels
- *
- * The covariance is modelled as alpha*sqrt(var_x*var_y)
- *
- * @note We provide a single-operand version for when the right-hand-side of an expression is a scalar, not a
- * masked pixel,
- */
-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;
-};
-
 /// Class for representing Unary operations
 template <typename ExprT1, typename ImageBinOp, typename MaskBinOp, typename VarianceBinOp>
 class UnaryExpr {
@@ -523,9 +491,7 @@ class BinaryExpr<ExprT1, double, ImageBinOp, MaskBinOp, VarianceBinOp> final {
 
 /// Template for -e1
 template <typename ExprT1>
-UnaryExpr<ExprT1, std::negate<typename exprTraits<ExprT1>::ImagePixelT>,
-          noop<typename exprTraits<ExprT1>::MaskPixelT>, noop<typename exprTraits<ExprT1>::VariancePixelT> >
-operator-(ExprT1 e1) {
+auto operator-(ExprT1 e1) {
     return UnaryExpr<ExprT1, std::negate<typename exprTraits<ExprT1>::ImagePixelT>,
                      noop<typename exprTraits<ExprT1>::MaskPixelT>,
                      noop<typename exprTraits<ExprT1>::VariancePixelT> >(e1);
@@ -539,31 +505,28 @@ struct is_specialization : std::false_type {};
 template <template <class...> class Template, class... Args>
 struct is_specialization<Template<Args...>, Template> : std::true_type {};
 // avoid invocation with pybind11::detail::descr::operator+ for pybind11 >= 2.3.0
-template <typename ExprT1, typename ExprT2,
-   typename = std::enable_if_t <
-          ( std::is_integral<ExprT1>::value  ||
-            is_specialization<ExprT1,BinaryExpr>{} ||
-            is_specialization<ExprT1,SinglePixel>{} ||
-            is_specialization<ExprT1,Pixel>{})
-            &&
-          ( std::is_integral<ExprT2>::value  ||
-            is_specialization<ExprT2,BinaryExpr>{} ||
-            is_specialization<ExprT2,SinglePixel>{} ||
-            is_specialization<ExprT2,Pixel>{} )
-         >
->
-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) {
+template <typename T>
+inline constexpr bool is_valid_operand = std::is_integral<T>::value  ||
+            std::is_floating_point_v<T> ||
+            is_specialization<T, BinaryExpr>{} ||
+            is_specialization<T, SinglePixel>{} ||
+            is_specialization<T, Pixel>{};
+template <typename T1, typename T2>
+constexpr bool valid_binary_operands = is_valid_operand<T1> && is_valid_operand<T2>;
+
+template <typename T1, typename T2>
+using valid_operands = std::enable_if_t<valid_binary_operands<T1, T2>>;
+
+template <typename ExprT1, typename ExprT2, typename = valid_operands<ExprT1, ExprT2>>
+auto 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 for e1 += e2
-template <typename ExprT1, typename ExprT2>
-ExprT1 operator+=(ExprT1& e1, ExprT2 e2) {
+template <typename ExprT1, typename ExprT2, typename = valid_operands<ExprT1, ExprT2>>
+ExprT1 &operator+=(ExprT1& e1, ExprT2 const &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);
@@ -573,26 +536,16 @@ ExprT1 operator+=(ExprT1& e1, ExprT2 e2) {
 //
 // 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>
+template <typename ExprT1, typename ExprT2, typename = valid_operands<ExprT1, 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);
-}
 }  // namespace
 
 /// Like operator+(), but assume that covariance's 2*alpha*sqrt(vx*vy)
-template <typename ExprT1, typename ExprT2>
+template <typename ExprT1, typename ExprT2, typename = valid_operands<ExprT1, ExprT2>>
 inline ExprT1 plus(
         ExprT1& lhs,        ///< Left hand value
         ExprT2 const& rhs,  ///< Right hand value
@@ -603,7 +556,7 @@ inline ExprT1 plus(
 
 //------------------------------------------
 /// Template to evaluate (e1 - e2)
-template <typename ExprT1, typename ExprT2>
+template <typename ExprT1, typename ExprT2, typename = valid_operands<ExprT1, ExprT2>>
 BinaryExpr<ExprT1, ExprT2, std::minus<typename exprTraits<ExprT1>::ImagePixelT>,
            bitwise_or<typename exprTraits<ExprT1>::MaskPixelT>,
            variance_plus<typename exprTraits<ExprT1>::VariancePixelT> >
@@ -614,7 +567,7 @@ operator-(ExprT1 e1, ExprT2 e2) {
 }
 
 /// Template to evaluate e1 -= e2
-template <typename ExprT1, typename ExprT2>
+template <typename ExprT1, typename ExprT2, typename = valid_operands<ExprT1, ExprT2>>
 ExprT1 operator-=(ExprT1& e1, ExprT2 e2) {
     e1 = BinaryExpr<ExprT1, ExprT2, std::minus<typename exprTraits<ExprT1>::ImagePixelT>,
                     bitwise_or<typename exprTraits<ExprT1>::MaskPixelT>,
@@ -624,7 +577,7 @@ ExprT1 operator-=(ExprT1& e1, ExprT2 e2) {
 
 //------------------------------------------
 /// Template to evaluate (e1 * e2)
-template <typename ExprT1, typename ExprT2>
+template <typename ExprT1, typename ExprT2, typename = valid_operands<ExprT1, ExprT2>>
 BinaryExpr<ExprT1, ExprT2, std::multiplies<typename exprTraits<ExprT1>::ImagePixelT>,
            bitwise_or<typename exprTraits<ExprT1>::MaskPixelT>,
            variance_multiplies<typename exprTraits<ExprT1>::VariancePixelT> >
@@ -635,7 +588,7 @@ operator*(ExprT1 e1, ExprT2 e2) {
 }
 
 /// Template to evaluate e1 *= e2
-template <typename ExprT1, typename ExprT2>
+template <typename ExprT1, typename ExprT2, typename = valid_operands<ExprT1, ExprT2>>
 ExprT1 operator*=(ExprT1& e1, ExprT2 e2) {
     e1 = BinaryExpr<ExprT1, ExprT2, std::multiplies<typename exprTraits<ExprT1>::ImagePixelT>,
                     bitwise_or<typename exprTraits<ExprT1>::MaskPixelT>,
@@ -645,7 +598,7 @@ ExprT1 operator*=(ExprT1& e1, ExprT2 e2) {
 
 //------------------------------------------
 /// Template to evaluate (e1 / e2)
-template <typename ExprT1, typename ExprT2>
+template <typename ExprT1, typename ExprT2, typename = valid_operands<ExprT1, ExprT2>>
 BinaryExpr<ExprT1, ExprT2, std::divides<typename exprTraits<ExprT1>::ImagePixelT>,
            bitwise_or<typename exprTraits<ExprT1>::MaskPixelT>,
            variance_divides<typename exprTraits<ExprT1>::VariancePixelT> >
@@ -656,7 +609,7 @@ operator/(ExprT1 e1, ExprT2 e2) {
 }
 
 /// Template to evaluate e1 /= e2
-template <typename ExprT1, typename ExprT2>
+template <typename ExprT1, typename ExprT2, typename = valid_operands<ExprT1, ExprT2>>
 ExprT1 operator/=(ExprT1& e1, ExprT2 e2) {
     e1 = BinaryExpr<ExprT1, ExprT2, std::divides<typename exprTraits<ExprT1>::ImagePixelT>,
                     bitwise_or<typename exprTraits<ExprT1>::MaskPixelT>,