Opencv meanStdDev mean函数详解

阅读了一下Opencv求均值和标准差的源码，总结下来，愿渡有缘人。

meanStdDev 函数在mean.dispatch.cpp中实现，计算均值和方差的公式如下：

 

meanStdDev这个函数的输入是src和mask，mask矩阵里为0位置的元素不被求和if( mask[i] ){才求和}；

这个函数的输出也是两个矩阵，这个输出矩阵的维数应该是channels * 1，每个通道存储相应的和（sum）或平方和（sqsum），如果空间太大，后面的空间会被设置为0.

 下面先说一下要理解这个函数的一些先决条件：

（1）NAryMatIterator类：这是一个N元矩阵迭代器，这个迭代器可同步处理多个矩阵，常用在处理src和mask的情况。每个NAryMatIterator::pannels是一个连续的地址空间，pannels是整体矩阵的切片。不连续的地址空间出现在如Mat(src, Rect(1, 1, 3, 3))的情况。就是用Rect截取了src矩阵的一部分，这个部分矩阵的地址空间不连续。

也有很多人写了关于这个类的详解，但是最终还是看源码上的注释才懂了的。

（2）这个函数调用了一个求和、求平方和的函数

static int sqsum8u(const uchar* src, const uchar* mask, int* sum, int* sqsum, int len, int cn)

是通过getSumSqrFunc()获取一个函数数组得到的。

这个函数的输入是src,mask

输出是：没有被mask的元素的和（保存在sum所指的变量内）、元素的平方和（保存在sqsum）

返回值：处理元素的长度len

如果src是多个通道，sum[0],sum[1]分别存储相应通道的和

（3）AutoBuffer：这个类简单，就相当于new

（4）这个函数里OCL/OVX/IPP：和图形加速、GPU相关，Opencv编译的时候应该没有使能这些选项，我也没有细究这些，高人可以补充一下。

（5）CV_INSTRUMENT_REGION()：这个宏Opencv里随处可见，我也没有细究，应该是标识下面的代码在链接过程应该放在哪个段内。（放在.code的指令区段内）

正题，解释一下这个函数内使用的变量的含义和难懂的语句：

（1）total：每个pannels的尺寸大小，这个是Pannels的总大小。

（2）blockSize：每次处理的元素个数，被初始化为total

（3）_buf（cn*4）：如果有3个通道，则对于通道一，使用_buf[0,3,6,9],通道二使用_buf[1,4,7,10]

其中_buf[0]保存最终的sum，_buf[3]保存最终的sqsum。_buf[6]和_buf[9]保存临时的sum和sqsum。

（4）intSumBlockSize：求和的时候，处理的元素个数超过intSumBlockSize的大小，就将求和的结果加到_buf[0]，intSumBlockSize被初始化为1<<15。这个值不能太大，否则求得的sum会超过int类型表示的范围。

（5） bool blockSum = depth <= CV_16S：

          bool blockSqSum = depth <= CV_8S：

这两条语句主要是因为求均值，那么元素大小不能超过2个字节，否则图像加起来会超过double（8字节）所能表示的范围；要计算方差，元素不能超过1字节，否则不计算方差。

正常情况这两个变量的值都应该为true

（6） int nz = func( ptrs[0], ptrs[1], (uchar*)sbuf, (uchar*)sqbuf, bsz, cn );

sbuf存储cn个通道的sum，sqbuf存储cn个通道的sqsum，返回处理元素的个数: nz = bsz

（7）if( blockSum && (count + blockSize >= intSumBlockSize || (i+1 >= it.nplanes && j+bsz >= total)) ){}

count为已经处理了的元素个数：已经处理的元素个数不能超过intSumBlockSize，因为sqsum8u返回的和是int类型的。（func就是static int sqsum8u( const uchar* src, const uchar* mask, int* sum, int* sqsum, int len, int cn )）

这里count + blockSize的意思是，如果马上要满了，下次循环处理的元素个数就超了intSumBlockSize，就把求得的和保存一下，保存到_buf[0]里面，这个变量是double类型的，所以足够存储一副图像的sum。

（8)nz0：count的值在处理的元素个数马上超过intSumBlockSize的时候就会清零，而nz0不会，nz0保存所有的处理过的元素个数。

（9）ptrs[0] += bsz*esz;
            if( ptrs[1] )
                ptrs[1] += bsz;

最初ptrs[0]指向第一个Pannels数据的首地址，ptrs[0] += bsz*esz就是指针在本Pannel内向后移动bsz*esz个字节。（因为ptrs[0]的类型是char *型）

ptrs[1] += bsz;不乘esz是因为如果要计算方差，元素必须是占1个字节

（10）后面的代码就简单了，就是先求均值和方差，然后通过一个2遍的for循环将计算值保存到输出参数内，for循环第一遍保存均值；第二遍保存标准差。

附上源码：

void meanStdDev(InputArray _src, OutputArray _mean, OutputArray _sdv, InputArray _mask)
{
    CV_INSTRUMENT_REGION();

    CV_Assert(!_src.empty());
    CV_Assert( _mask.empty() || _mask.type() == CV_8UC1 );

    CV_OCL_RUN(OCL_PERFORMANCE_CHECK(_src.isUMat()) && _src.dims() <= 2,
               ocl_meanStdDev(_src, _mean, _sdv, _mask))

    Mat src = _src.getMat(), mask = _mask.getMat();

    CV_OVX_RUN(!ovx::skipSmallImages<VX_KERNEL_MEAN_STDDEV>(src.cols, src.rows),
               openvx_meanStdDev(src, _mean, _sdv, mask))

    CV_IPP_RUN(IPP_VERSION_X100 >= 700, ipp_meanStdDev(src, _mean, _sdv, mask));

    int k, cn = src.channels(), depth = src.depth();

    SumSqrFunc func = getSumSqrFunc(depth);

    CV_Assert( func != 0 );

    const Mat* arrays[] = {&src, &mask, 0};
    uchar* ptrs[2] = {};
    NAryMatIterator it(arrays, ptrs);
    int total = (int)it.size, blockSize = total, intSumBlockSize = 0;
    int j, count = 0, nz0 = 0;
    AutoBuffer<double> _buf(cn*4);
    double *s = (double*)_buf.data(), *sq = s + cn;
    int *sbuf = (int*)s, *sqbuf = (int*)sq;
    bool blockSum = depth <= CV_16S, blockSqSum = depth <= CV_8S;
    size_t esz = 0;

	for( k = 0; k < cn; k++ )
	    s[k] = sq[k] = 0;

    if( blockSum )
    {
        intSumBlockSize = 1 << 15;
        blockSize = std::min(blockSize, intSumBlockSize);
        sbuf = (int*)(sq + cn);
        if( blockSqSum )
            sqbuf = sbuf + cn;
        for( k = 0; k < cn; k++ )
            sbuf[k] = sqbuf[k] = 0;
        esz = src.elemSize();
    }

    for( size_t i = 0; i < it.nplanes; i++, ++it )
    {
        for( j = 0; j < total; j += blockSize )
        {
            int bsz = std::min(total - j, blockSize);
            int nz = func( ptrs[0], ptrs[1], (uchar*)sbuf, (uchar*)sqbuf, bsz, cn );
            count += nz;
            nz0 += nz;
            if( blockSum && (count + blockSize >= intSumBlockSize || (i+1 >= it.nplanes && j+bsz >= total)) )
            {
                for( k = 0; k < cn; k++ )
                {
                    s[k] += sbuf[k];
                    sbuf[k] = 0;
                }
                if( blockSqSum )
                {
                    for( k = 0; k < cn; k++ )
                    {
                        sq[k] += sqbuf[k];
                        sqbuf[k] = 0;
                    }
                }
                count = 0;
            }
            ptrs[0] += bsz*esz;
            if( ptrs[1] )
                ptrs[1] += bsz;
        }
    }

    double scale = nz0 ? 1./nz0 : 0.;
    for( k = 0; k < cn; k++ )
    {
        s[k] *= scale;
        sq[k] = std::sqrt(std::max(sq[k]*scale - s[k]*s[k], 0.));
    }

    for( j = 0; j < 2; j++ )
    {
        const double* sptr = j == 0 ? s : sq;
        _OutputArray _dst = j == 0 ? _mean : _sdv;
        if( !_dst.needed() )
            continue;

        if( !_dst.fixedSize() )
            _dst.create(cn, 1, CV_64F, -1, true);
        Mat dst = _dst.getMat();
        int dcn = (int)dst.total();
        CV_Assert( dst.type() == CV_64F && dst.isContinuous() &&
                   (dst.cols == 1 || dst.rows == 1) && dcn >= cn );
        double* dptr = dst.ptr<double>();
        for( k = 0; k < cn; k++ )
            dptr[k] = sptr[k];
        for( ; k < dcn; k++ )
            dptr[k] = 0;
    }
}