Add multi-channel mask support to mean, meanStdDev and setTo

This adds the possibility to use multi-channel masks for the functions
cv::mean, cv::meanStdDev and the method Mat::setTo. The tests have now a
probability to use multi-channel masks for operations that support them.
This also includes Mat::copyTo, which supported multi-channel masks
before, but there was no test confirming this.

modules/core/include/opencv2/core.hpp

@@ -608,7 +608,7 @@ CV_EXPORTS_W void findNonZero( InputArray src, OutputArray idx );
 
 The function cv::mean calculates the mean value M of array elements,
 independently for each channel, and return it:
-\f[\begin{array}{l} N =  \sum _{I: \; \texttt{mask} (I) \ne 0} 1 \\ M_c =  \left ( \sum _{I: \; \texttt{mask} (I) \ne 0}{ \texttt{mtx} (I)_c} \right )/N \end{array}\f]
+\f[\begin{array}{l} N_c =  \sum _{I: \; {\texttt{mask} (I)_c} \ne 0} 1 \\ M_c =  \left ( \sum _{I: \; {\texttt{mask} (I)_c} \ne 0}{ \texttt{src} (I)_c} \right )/N_c \end{array}\f]
 When all the mask elements are 0's, the function returns Scalar::all(0)
 @param src input array that should have from 1 to 4 channels so that the result can be stored in
 Scalar_ .
@@ -622,7 +622,7 @@ CV_EXPORTS_W Scalar mean(InputArray src, InputArray mask = noArray());
 The function cv::meanStdDev calculates the mean and the standard deviation M
 of array elements independently for each channel and returns it via the
 output parameters:
-\f[\begin{array}{l} N =  \sum _{I, \texttt{mask} (I)  \ne 0} 1 \\ \texttt{mean} _c =  \frac{\sum_{ I: \; \texttt{mask}(I) \ne 0} \texttt{src} (I)_c}{N} \\ \texttt{stddev} _c =  \sqrt{\frac{\sum_{ I: \; \texttt{mask}(I) \ne 0} \left ( \texttt{src} (I)_c -  \texttt{mean} _c \right )^2}{N}} \end{array}\f]
+\f[\begin{array}{l} N_c =  \sum _{I, {\texttt{mask} (I)_c}  \ne 0} 1 \\ \texttt{mean} _c =  \frac{\sum_{ I: \; {\texttt{mask} (I)_c} \ne 0} \texttt{src} (I)_c}{N_c} \\ \texttt{stddev} _c =  \sqrt{\frac{\sum_{ I: \; {\texttt{mask} (I)_c} \ne 0} \left ( \texttt{src} (I)_c -  \texttt{mean} _c \right )^2}{N_c}} \end{array}\f]
 When all the mask elements are 0's, the function returns
 mean=stddev=Scalar::all(0).
 @note The calculated standard deviation is only the diagonal of the

modules/core/include/opencv2/core/mat.hpp

@@ -1192,8 +1192,8 @@ public:
     /** @overload
     @param m Destination matrix. If it does not have a proper size or type before the operation, it is
     reallocated.
-    @param mask Operation mask. Its non-zero elements indicate which matrix elements need to be copied.
-    The mask has to be of type CV_8U and can have 1 or multiple channels.
+    @param mask Operation mask of the same size as \*this. Its non-zero elements indicate which matrix
+    elements need to be copied. The mask has to be of type CV_8U and can have 1 or multiple channels.
     */
     void copyTo( OutputArray m, InputArray mask ) const;
 
@@ -1229,7 +1229,8 @@ public:
 
     This is an advanced variant of the Mat::operator=(const Scalar& s) operator.
     @param value Assigned scalar converted to the actual array type.
-    @param mask Operation mask of the same size as \*this.
+    @param mask Operation mask of the same size as \*this. Its non-zero elements indicate which matrix
+    elements need to be copied. The mask has to be of type CV_8U and can have 1 or multiple channels
      */
     Mat& setTo(InputArray value, InputArray mask=noArray());
 

modules/core/src/copy.cpp

@@ -334,7 +334,7 @@ static bool ipp_copyTo(const Mat &src, Mat &dst, const Mat &mask)
 #ifdef HAVE_IPP_IW
     CV_INSTRUMENT_REGION_IPP()
 
-    if(mask.channels() > 1 && mask.depth() != CV_8U)
+    if(mask.channels() > 1 || mask.depth() != CV_8U)
         return false;
 
     if (src.dims <= 2)
@@ -510,20 +510,23 @@ Mat& Mat::setTo(InputArray _value, InputArray _mask)
     Mat value = _value.getMat(), mask = _mask.getMat();
 
     CV_Assert( checkScalar(value, type(), _value.kind(), _InputArray::MAT ));
-    CV_Assert( mask.empty() || (mask.type() == CV_8U && size == mask.size) );
+    int cn = channels(), mcn = mask.channels();
+    CV_Assert( mask.empty() || (mask.depth() == CV_8U && (mcn == 1 || mcn == cn) && size == mask.size) );
 
     CV_IPP_RUN_FAST(ipp_Mat_setTo_Mat(*this, value, mask), *this)
 
-    size_t esz = elemSize();
+    size_t esz = mcn > 1 ? elemSize1() : elemSize();
     BinaryFunc copymask = getCopyMaskFunc(esz);
 
     const Mat* arrays[] = { this, !mask.empty() ? &mask : 0, 0 };
     uchar* ptrs[2]={0,0};
     NAryMatIterator it(arrays, ptrs);
-    int totalsz = (int)it.size, blockSize0 = std::min(totalsz, (int)((BLOCK_SIZE + esz-1)/esz));
+    int totalsz = (int)it.size*mcn;
+    int blockSize0 = std::min(totalsz, (int)((BLOCK_SIZE + esz-1)/esz));
+    blockSize0 -= blockSize0 % mcn;    // must be divisible without remainder for unrolling and advancing
     AutoBuffer<uchar> _scbuf(blockSize0*esz + 32);
     uchar* scbuf = alignPtr((uchar*)_scbuf, (int)sizeof(double));
-    convertAndUnrollScalar( value, type(), scbuf, blockSize0 );
+    convertAndUnrollScalar( value, type(), scbuf, blockSize0/mcn );
 
     for( size_t i = 0; i < it.nplanes; i++, ++it )
     {

modules/core/test/test_arithm.cpp

-#include "test_precomp.hpp"
+#include "test_precomp.hpp"
 #include <cmath>
 
 using namespace cv;
@@ -15,7 +15,7 @@ const int ARITHM_MAX_SIZE_LOG = 10;
 
 struct BaseElemWiseOp
 {
-    enum { FIX_ALPHA=1, FIX_BETA=2, FIX_GAMMA=4, REAL_GAMMA=8, SUPPORT_MASK=16, SCALAR_OUTPUT=32 };
+    enum { FIX_ALPHA=1, FIX_BETA=2, FIX_GAMMA=4, REAL_GAMMA=8, SUPPORT_MASK=16, SCALAR_OUTPUT=32, SUPPORT_MULTICHANNELMASK=64 };
     BaseElemWiseOp(int _ninputs, int _flags, double _alpha, double _beta,
                    Scalar _gamma=Scalar::all(0), int _context=1)
     : ninputs(_ninputs), flags(_flags), alpha(_alpha), beta(_beta), gamma(_gamma), context(_context) {}
@@ -467,7 +467,7 @@ struct CmpSOp : public BaseElemWiseOp
 
 struct CopyOp : public BaseElemWiseOp
 {
-    CopyOp() : BaseElemWiseOp(1, FIX_ALPHA+FIX_BETA+FIX_GAMMA+SUPPORT_MASK, 1, 1, Scalar::all(0)) {  }
+    CopyOp() : BaseElemWiseOp(1, FIX_ALPHA+FIX_BETA+FIX_GAMMA+SUPPORT_MASK+SUPPORT_MULTICHANNELMASK, 1, 1, Scalar::all(0)) {  }
     void op(const vector<Mat>& src, Mat& dst, const Mat& mask)
     {
         src[0].copyTo(dst, mask);
@@ -489,7 +489,7 @@ struct CopyOp : public BaseElemWiseOp
 
 struct SetOp : public BaseElemWiseOp
 {
-    SetOp() : BaseElemWiseOp(0, FIX_ALPHA+FIX_BETA+SUPPORT_MASK, 1, 1, Scalar::all(0)) {}
+    SetOp() : BaseElemWiseOp(0, FIX_ALPHA+FIX_BETA+SUPPORT_MASK+SUPPORT_MULTICHANNELMASK, 1, 1, Scalar::all(0)) {}
     void op(const vector<Mat>&, Mat& dst, const Mat& mask)
     {
         dst.setTo(gamma, mask);
@@ -1162,7 +1162,7 @@ struct CartToPolarToCartOp : public BaseElemWiseOp
 
 struct MeanOp : public BaseElemWiseOp
 {
-    MeanOp() : BaseElemWiseOp(1, FIX_ALPHA+FIX_BETA+FIX_GAMMA+SUPPORT_MASK+SCALAR_OUTPUT, 1, 1, Scalar::all(0))
+    MeanOp() : BaseElemWiseOp(1, FIX_ALPHA+FIX_BETA+FIX_GAMMA+SUPPORT_MASK+SUPPORT_MULTICHANNELMASK+SCALAR_OUTPUT, 1, 1, Scalar::all(0))
     {
         context = 3;
     };
@@ -1244,7 +1244,7 @@ struct MeanStdDevOp : public BaseElemWiseOp
     Scalar sqmeanRef;
     int cn;
 
-    MeanStdDevOp() : BaseElemWiseOp(1, FIX_ALPHA+FIX_BETA+FIX_GAMMA+SUPPORT_MASK+SCALAR_OUTPUT, 1, 1, Scalar::all(0))
+    MeanStdDevOp() : BaseElemWiseOp(1, FIX_ALPHA+FIX_BETA+FIX_GAMMA+SUPPORT_MASK+SUPPORT_MULTICHANNELMASK+SCALAR_OUTPUT, 1, 1, Scalar::all(0))
     {
         cn = 0;
         context = 7;
@@ -1394,7 +1394,8 @@ TEST_P(ElemWiseTest, accuracy)
         op->getRandomSize(rng, size);
         int type = op->getRandomType(rng);
         int depth = CV_MAT_DEPTH(type);
-        bool haveMask = (op->flags & cvtest::BaseElemWiseOp::SUPPORT_MASK) != 0 && rng.uniform(0, 4) == 0;
+        bool haveMask = ((op->flags & cvtest::BaseElemWiseOp::SUPPORT_MASK) != 0
+                || (op->flags & cvtest::BaseElemWiseOp::SUPPORT_MULTICHANNELMASK) != 0) && rng.uniform(0, 4) == 0;
 
         double minval=0, maxval=0;
         op->getValueRange(depth, minval, maxval);
@@ -1403,8 +1404,12 @@ TEST_P(ElemWiseTest, accuracy)
         for( i = 0; i < ninputs; i++ )
             src[i] = cvtest::randomMat(rng, size, type, minval, maxval, true);
         Mat dst0, dst, mask;
-        if( haveMask )
-            mask = cvtest::randomMat(rng, size, CV_8U, 0, 2, true);
+        if( haveMask ) {
+            bool multiChannelMask = (op->flags & cvtest::BaseElemWiseOp::SUPPORT_MULTICHANNELMASK) != 0
+                    && rng.uniform(0, 2) == 0;
+            int masktype = CV_8UC(multiChannelMask ? CV_MAT_CN(type) : 1);
+            mask = cvtest::randomMat(rng, size, masktype, 0, 2, true);
+        }
 
         if( (haveMask || ninputs == 0) && !(op->flags & cvtest::BaseElemWiseOp::SCALAR_OUTPUT))
         {

modules/ts/src/ts_func.cpp

@@ -353,26 +353,38 @@ void copy(const Mat& src, Mat& dst, const Mat& mask, bool invertMask)
         return;
     }
 
-    CV_Assert( src.size == mask.size && mask.type() == CV_8U );
+    int mcn = mask.channels();
+    CV_Assert( src.size == mask.size && mask.depth() == CV_8U
+               && (mcn == 1 || mcn == src.channels()) );
 
     const Mat *arrays[]={&src, &dst, &mask, 0};
     Mat planes[3];
 
     NAryMatIterator it(arrays, planes);
-    size_t j, k, elemSize = src.elemSize(), total = planes[0].total();
+    size_t j, k, elemSize = src.elemSize(), maskElemSize = mask.elemSize(), total = planes[0].total();
     size_t i, nplanes = it.nplanes;
+    size_t elemSize1 = src.elemSize1();
 
     for( i = 0; i < nplanes; i++, ++it)
     {
         const uchar* sptr = planes[0].ptr();
         uchar* dptr = planes[1].ptr();
         const uchar* mptr = planes[2].ptr();
-
-        for( j = 0; j < total; j++, sptr += elemSize, dptr += elemSize )
+        for( j = 0; j < total; j++, sptr += elemSize, dptr += elemSize, mptr += maskElemSize )
         {
-            if( (mptr[j] != 0) ^ invertMask )
-                for( k = 0; k < elemSize; k++ )
-                    dptr[k] = sptr[k];
+            if( mcn == 1)
+            {
+                if( (mptr[0] != 0) ^ invertMask )
+                    for( k = 0; k < elemSize; k++ )
+                        dptr[k] = sptr[k];
+            }
+            else
+            {
+                for( int c = 0; c < mcn; c++ )
+                    if( (mptr[c] != 0) ^ invertMask )
+                        for( k = 0; k < elemSize1; k++ )
+                            dptr[k + c * elemSize1] = sptr[k + c * elemSize1];
+            }
         }
     }
 }
@@ -414,25 +426,37 @@ void set(Mat& dst, const Scalar& gamma, const Mat& mask)
         return;
     }
 
-    CV_Assert( dst.size == mask.size && mask.type() == CV_8U );
+    int cn = dst.channels(), mcn = mask.channels();
+    CV_Assert( dst.size == mask.size && (mcn == 1 || mcn == cn) );
 
     const Mat *arrays[]={&dst, &mask, 0};
     Mat planes[2];
 
     NAryMatIterator it(arrays, planes);
-    size_t j, k, elemSize = dst.elemSize(), total = planes[0].total();
+    size_t j, k, elemSize = dst.elemSize(), maskElemSize = mask.elemSize(), total = planes[0].total();
     size_t i, nplanes = it.nplanes;
+    size_t elemSize1 = dst.elemSize1();
 
     for( i = 0; i < nplanes; i++, ++it)
     {
         uchar* dptr = planes[0].ptr();
         const uchar* mptr = planes[1].ptr();
 
-        for( j = 0; j < total; j++, dptr += elemSize )
+        for( j = 0; j < total; j++, dptr += elemSize, mptr += maskElemSize )
         {
-            if( mptr[j] )
-                for( k = 0; k < elemSize; k++ )
-                    dptr[k] = gptr[k];
+            if( mcn == 1)
+            {
+                if( mptr[0] )
+                    for( k = 0; k < elemSize; k++ )
+                        dptr[k] = gptr[k];
+            }
+            else
+            {
+                for( int c = 0; c < mcn; c++ )
+                    if( mptr[c] )
+                        for( k = 0; k < elemSize1; k++ )
+                            dptr[k + c * elemSize1] = gptr[k + c * elemSize1];
+            }
         }
     }
 }
@@ -2574,11 +2598,11 @@ void divide(const Mat& src1, const Mat& src2, Mat& dst, double scale)
 
 
 template<typename _Tp> static void
-mean_(const _Tp* src, const uchar* mask, size_t total, int cn, Scalar& sum, int& nz)
+mean_(const _Tp* src, const uchar* mask, size_t total, int cn, int mcn, Scalar& sum, Scalar_<int>& nz)
 {
     if( !mask )
     {
-        nz += (int)total;
+        nz += Scalar_<int>::all((int)total);
         total *= cn;
         for( size_t i = 0; i < total; i += cn )
         {
@@ -2586,23 +2610,41 @@ mean_(const _Tp* src, const uchar* mask, size_t total, int cn, Scalar& sum, int&
                 sum[c] += src[i + c];
         }
     }
-    else
+    else if( mcn == 1 )
     {
         for( size_t i = 0; i < total; i++ )
             if( mask[i] )
             {
-                nz++;
                 for( int c = 0; c < cn; c++ )
+                {
+                    nz[c]++;
                     sum[c] += src[i*cn + c];
+                }
+            }
+    }
+    else
+    {
+        total *= cn;
+        for( size_t i = 0; i < total; i += cn )
+        {
+            for( int c = 0; c < cn; c++ )
+            {
+                if( mask[i + c] )
+                {
+                    nz[c]++;
+                    sum[c] += src[i + c];
+                }
             }
+        }
     }
 }
 
 Scalar mean(const Mat& src, const Mat& mask)
 {
-    CV_Assert(mask.empty() || (mask.type() == CV_8U && mask.size == src.size));
+    CV_Assert(mask.empty() || (mask.depth() == CV_8U && mask.size == src.size &&
+                               (mask.channels() == 1 || mask.channels() == src.channels())));
     Scalar sum;
-    int nz = 0;
+    Scalar_<int> nz = Scalar_<int>::all(0);
 
     const Mat *arrays[]={&src, &mask, 0};
     Mat planes[2];
@@ -2610,7 +2652,7 @@ Scalar mean(const Mat& src, const Mat& mask)
     NAryMatIterator it(arrays, planes);
     size_t total = planes[0].total();
     size_t i, nplanes = it.nplanes;
-    int depth = src.depth(), cn = src.channels();
+    int c, depth = src.depth(), cn = src.channels(), mcn = mask.channels();
 
     for( i = 0; i < nplanes; i++, ++it )
     {
@@ -2620,32 +2662,34 @@ Scalar mean(const Mat& src, const Mat& mask)
         switch( depth )
         {
         case CV_8U:
-            mean_((const uchar*)sptr, mptr, total, cn, sum, nz);
+            mean_((const uchar*)sptr, mptr, total, cn, mcn, sum, nz);
             break;
         case CV_8S:
-            mean_((const schar*)sptr, mptr, total, cn, sum, nz);
+            mean_((const schar*)sptr, mptr, total, cn, mcn, sum, nz);
             break;
         case CV_16U:
-            mean_((const ushort*)sptr, mptr, total, cn, sum, nz);
+            mean_((const ushort*)sptr, mptr, total, cn, mcn, sum, nz);
             break;
         case CV_16S:
-            mean_((const short*)sptr, mptr, total, cn, sum, nz);
+            mean_((const short*)sptr, mptr, total, cn, mcn, sum, nz);
             break;
         case CV_32S:
-            mean_((const int*)sptr, mptr, total, cn, sum, nz);
+            mean_((const int*)sptr, mptr, total, cn, mcn, sum, nz);
             break;
         case CV_32F:
-            mean_((const float*)sptr, mptr, total, cn, sum, nz);
+            mean_((const float*)sptr, mptr, total, cn, mcn, sum, nz);
             break;
         case CV_64F:
-            mean_((const double*)sptr, mptr, total, cn, sum, nz);
+            mean_((const double*)sptr, mptr, total, cn, mcn, sum, nz);
             break;
         default:
             CV_Error(Error::StsUnsupportedFormat, "");
         }
     }
 
-    return sum * (1./std::max(nz, 1));
+    for( c = 0; c < cn; c++ )
+        sum[c] *= (1./std::max(nz[c], 1));
+    return sum;
 }