OpenCV源码解析：协方差矩阵的计算--calcCovarMatrix
协方差矩阵
在统计学与概率论中，协方差是指两个向量元素之间的相关性。
设为n维随机变量
方差的定义为：
当存在两个随机变量X，Y时，其各个维度偏离其均值的程度就可以用协方差来定义：
在物理上的理解，你可以认为协方差是指两个向量之相互影响的程度，单从数值上来看，协方差的数值越大，表示两个变量对其均值的变化同向的程度越大。
当随机变量有多个的时候，一般不再使用X，Y这样的表述，而是使用X1,X2,…Xn来描述一个多维向量
协方差矩阵定义为，
展开后的形式就是
这很直观，运算起来也不难。
注意，这是一个对称矩阵，在计算机处理中，一般协方差矩阵的计算是这样的：先让样本矩阵中心化，即每一维度减去该维度的均值（这样一来，每一维度上的均值为0），然后直接使用新得到的样本矩阵乘以它的转置，最后除以(N-1)。OpenCV正是采用了这种算法。
OpenCV协方差矩阵的计算详解
下面为了方便，我们举一个4维向量，每个向量有5个样本的例子，
double data[5][4];在后面的讲解中，我们对这个矩阵进行协方差计算。理解了这个过程，其他类似的就都很容易了。相信你看过之后，也能自己独立写出协方差运算的C/C++实现。
cv::calcCovarMatrix源码如下，
void cv::calcCovarMatrix( InputArray _src, OutputArray _covar, InputOutputArray _mean, int flags, int ctype )
{
CV_INSTRUMENT_REGION()
if(_src.kind() == _InputArray::STD_VECTOR_MAT || _src.kind() == _InputArray::STD_ARRAY_MAT)
{
std::vector<cv::Mat> src;
_src.getMatVector(src);
CV_Assert( src.size() > 0 );
Size size = src[0].size();
int type = src[0].type();
ctype = std::max(std::max(CV_MAT_DEPTH(ctype >= 0 ? ctype : type), _mean.depth()), CV_32F);
Mat _data(static_cast<int>(src.size()), size.area(), type);
int i = 0;
for(std::vector<cv::Mat>::iterator each = src.begin(); each != src.end(); ++each, ++i )
{
CV_Assert( (*each).size() == size, (*each).type() == type );
Mat dataRow(size.height, size.width, type, _data.ptr(i));
(*each).copyTo(dataRow);
}
Mat mean;
if( (flags & CV_COVAR_USE_AVG) != 0 )
{
CV_Assert( _mean.size() == size );
if( mean.type() != ctype )
{
mean = _mean.getMat();
_mean.create(mean.size(), ctype);
Mat tmp = _mean.getMat();
mean.convertTo(tmp, ctype);
mean = tmp;
}
mean = _mean.getMat().reshape(1, 1);
}
calcCovarMatrix( _data, _covar, mean, (flags & ~(CV_COVAR_ROWS|CV_COVAR_COLS)) | CV_COVAR_ROWS, ctype );
if( (flags & CV_COVAR_USE_AVG) == 0 )
{
mean = mean.reshape(1, size.height);
mean.copyTo(_mean);
}
return;
}
Mat data = _src.getMat(), mean;
CV_Assert( ((flags & CV_COVAR_ROWS) != 0) ^ ((flags & CV_COVAR_COLS) != 0) );
bool takeRows = (flags & CV_COVAR_ROWS) != 0; // 向量样本数是否为行数
int type = data.type(); // 数据类型
int nsamples = takeRows ? data.rows : data.cols; //样本个数
CV_Assert( nsamples > 0 );
Size size = takeRows ? Size(data.cols, 1) : Size(1, data.rows); // 用于计算均值的向量的大小
if( (flags & CV_COVAR_USE_AVG) != 0 )
{
mean = _mean.getMat();
ctype = std::max(std::max(CV_MAT_DEPTH(ctype >= 0 ? ctype : type), mean.depth()), CV_32F);
CV_Assert( mean.size() == size );
if( mean.type() != ctype )
{
_mean.create(mean.size(), ctype);
Mat tmp = _mean.getMat();
mean.convertTo(tmp, ctype);
mean = tmp;
}
}
else
{
ctype = std::max(CV_MAT_DEPTH(ctype >= 0 ? ctype : type), CV_32F);
reduce( _src, _mean, takeRows ? 0 : 1, CV_REDUCE_AVG, ctype );
mean = _mean.getMat();
}
mulTransposed( data, _covar, ((flags & CV_COVAR_NORMAL) == 0) ^ takeRows,
mean, (flags & CV_COVAR_SCALE) != 0 ? 1./nsamples : 1, ctype );
}其过程是这样的：先判断是否是向量或数组，如果是，就按向量或数组处理。这里不符合条件，直接跳过。
真正的计算从reduce函数开始的，reduce函数计算了所有向量的平均值，得到平均值向量mean。最后调用mulTransposed完成协方差矩阵的计算。
先来看reduce的源代码，
void cv::reduce(InputArray _src, OutputArray _dst, int dim, int op, int dtype)
{
CV_INSTRUMENT_REGION()
CV_Assert( _src.dims() <= 2 );
int op0 = op;
int stype = _src.type(), sdepth = CV_MAT_DEPTH(stype), cn = CV_MAT_CN(stype);
if( dtype < 0 )
dtype = _dst.fixedType() ? _dst.type() : stype;
dtype = CV_MAKETYPE(dtype >= 0 ? dtype : stype, cn);
int ddepth = CV_MAT_DEPTH(dtype);
CV_Assert( cn == CV_MAT_CN(dtype) );
CV_Assert( op == CV_REDUCE_SUM || op == CV_REDUCE_MAX ||
op == CV_REDUCE_MIN || op == CV_REDUCE_AVG );
CV_OCL_RUN(_dst.isUMat(),
ocl_reduce(_src, _dst, dim, op, op0, stype, dtype))
// Fake reference to source. Resolves issue 8693 in case of src == dst.
UMat srcUMat;
if (_src.isUMat())
srcUMat = _src.getUMat();
Mat src = _src.getMat();
// 从calcCovarMatrix中，takeRows = 1, dim = 0；否则为 1
_dst.create(dim == 0 ? 1 : src.rows, dim == 0 ? src.cols : 1, dtype);
Mat dst = _dst.getMat(), temp = dst;
// 设置操作模式这些模式包括
/**
#define CV_REDUCE_SUM 0
#define CV_REDUCE_AVG 1
#define CV_REDUCE_MAX 2
#define CV_REDUCE_MIN 3
*/
if( op == CV_REDUCE_AVG )
{
op = CV_REDUCE_SUM;
if( sdepth < CV_32S && ddepth < CV_32S )
{
temp.create(dst.rows, dst.cols, CV_32SC(cn));
ddepth = CV_32S;
}
}
// 根据数据类型和操作模式，设置操作方程
ReduceFunc func = 0;
if( dim == 0 )
{
if( op == CV_REDUCE_SUM )
{
if(sdepth == CV_8U && ddepth == CV_32S)
func = GET_OPTIMIZED(reduceSumR8u32s);
else if(sdepth == CV_8U && ddepth == CV_32F)
func = reduceSumR8u32f;
else if(sdepth == CV_8U && ddepth == CV_64F)
func = reduceSumR8u64f;
else if(sdepth == CV_16U && ddepth == CV_32F)
func = reduceSumR16u32f;
else if(sdepth == CV_16U && ddepth == CV_64F)
func = reduceSumR16u64f;
else if(sdepth == CV_16S && ddepth == CV_32F)
func = reduceSumR16s32f;
else if(sdepth == CV_16S && ddepth == CV_64F)
func = reduceSumR16s64f;
else if(sdepth == CV_32F && ddepth == CV_32F)
func = GET_OPTIMIZED(reduceSumR32f32f);
else if(sdepth == CV_32F && ddepth == CV_64F)
func = reduceSumR32f64f;
else if(sdepth == CV_64F && ddepth == CV_64F)
func = reduceSumR64f64f;
}
else 。。。
}
else
{
。。。
}
if( !func )
CV_Error( CV_StsUnsupportedFormat,
"Unsupported combination of input and output array formats" );
func( src, temp ); // 执行ReduceR的操作
if( op0 == CV_REDUCE_AVG )
temp.convertTo(dst, dst.type(), 1./(dim == 0 ? src.rows : src.cols));
}
在reduce这个函数中，先通过func( src, temp )计算了所有向量的样本的和， 最后得到的结果就是这个_dst（=temp）向量，其中的每一个值是矩阵中对应向量所有样本的和，比如在CV_COVAR_ROWS模式下，data[5][4]矩阵有4个向量，每个向量的5个样本值相加后会得到一个值，整个5x4的矩阵就得到4个平均值，构成一个1x4的向量。
然后用convertTo求这些和的平均值。convertTo是一个比较简单的矩阵运算，主要负责数据类型转换（注意区别cvtColor()函数是负责转换不同通道的），这里大材小用了一下，用来对前面的1x4的向量和取平均值。
接下来，我们看一下cv::calcCovarMatrix中的cv::mulTransposed函数，它负责计算矩阵的转置并计算最后的乘积。其中参数ata 代表乘积的顺序，为1的时候，转置矩阵在前面：Calculates (A-delta)T*(A-delta) ；ata为0的时候，转置矩阵在后面，进行的计算是 (A-delta)*(A-delta)T，
void cv::mulTransposed( InputArray _src, OutputArray _dst, bool ata,
InputArray _delta, double scale, int dtype )
{
CV_INSTRUMENT_REGION()
Mat src = _src.getMat(), delta = _delta.getMat();
const int gemm_level = 100; // boundary above which GEMM is faster.
int stype = src.type();
dtype = std::max(std::max(CV_MAT_DEPTH(dtype >= 0 ? dtype : stype), delta.depth()), CV_32F);
CV_Assert( src.channels() == 1 );
if( !delta.empty() )
{
CV_Assert( delta.channels() == 1,
(delta.rows == src.rows || delta.rows == 1),
(delta.cols == src.cols || delta.cols == 1));
if( delta.type() != dtype )
delta.convertTo(delta, dtype);
}
int dsize = ata ? src.cols : src.rows;
_dst.create( dsize, dsize, dtype );
Mat dst = _dst.getMat();
if( src.data == dst.data || (stype == dtype &&
(dst.cols >= gemm_level && dst.rows >= gemm_level &&
src.cols >= gemm_level && src.rows >= gemm_level)))
{
Mat src2;
const Mat* tsrc = &src;
if( !delta.empty() )
{
if( delta.size() == src.size() )
subtract( src, delta, src2 );
else
{
repeat(delta, src.rows/delta.rows, src.cols/delta.cols, src2);
subtract( src, src2, src2 );
}
tsrc = &src2;
}
gemm( *tsrc, *tsrc, scale, Mat(), 0, dst, ata ? GEMM_1_T : GEMM_2_T );
}
else
{
MulTransposedFunc func = 0;
if。。。
else if(stype == CV_64F && dtype == CV_64F)
{
if(ata)
func = MulTransposedR<double,double>;
// 计算对角线及其右侧的元素, 先转置再计算，A^T*A
else
func = MulTransposedL<double,double>;
}
if( !func )
CV_Error( CV_StsUnsupportedFormat, "" );
func( src, dst, delta, scale ); // 调用上面得到的函数开始计算
completeSymm( dst, false ); // 完成对角线拷贝
}
}
这个函数本身没有太多内容，真正的计算是通过调用MulTransposedR实现的。MulTransposedR完成矩阵转置，并进行协方差矩阵中对角线及右上角元素的计算。最后通过
completeSymm( dst, false );完成对角线拷贝，得到一个完整的协方差矩阵。
下面重点看一下MulTransposedR，先列一下其计算过程，
因为矩阵是对称的，所以只计算了对角线右上侧的部分，例如，对于这个dst[4][4]的协方差矩阵，计算顺序为
第1轮：dst[0][0], dst[0][1], dst[0][2], dst[0][3]
第2轮： dst[1][1], dst[1][2], dst[1][3]
第3轮： dst[2][2], dst[2][3]
第4轮： dst[3][3]
在每一个协方差的计算中，都会把相应的样本值减去该向量的样本平均值再进行乘法计算，样本的平均值就是前面通过reduce计算出来的，在MulTransposedR中就是那个传入参数deltamat。 这部分是计算的核心代码，所以我贴出了全部的源码
template<typename sT, typename dT> static void
MulTransposedR( const Mat& srcmat, Mat& dstmat, const Mat& deltamat, double scale )
{
int i, j, k;
const sT* src = srcmat.ptr<sT>();
dT* dst = dstmat.ptr<dT>();
const dT* delta = deltamat.ptr<dT>();
size_t srcstep = srcmat.step/sizeof(src[0]);
size_t dststep = dstmat.step/sizeof(dst[0]);
size_t deltastep = deltamat.rows > 1 ? deltamat.step/sizeof(delta[0]) : 0;
int delta_cols = deltamat.cols;
Size size = srcmat.size();
dT* tdst = dst;
dT* col_buf = 0;
dT* delta_buf = 0;
int buf_size = size.height*sizeof(dT);
AutoBuffer<uchar> buf;
if( delta && delta_cols < size.width )
{
assert( delta_cols == 1 );
buf_size *= 5;
}
buf.allocate(buf_size);
col_buf = (dT*)(uchar*)buf;
if( delta && delta_cols < size.width )
{
delta_buf = col_buf + size.height;
for( i = 0; i < size.height; i++ )
delta_buf[i*4] = delta_buf[i*4+1] =
delta_buf[i*4+2] = delta_buf[i*4+3] = delta[i*deltastep];
delta = delta_buf;
deltastep = deltastep ? 4 : 0;
}
if( !delta )
for( i = 0; i < size.width; i++, tdst += dststep )
{
for( k = 0; k < size.height; k++ )
col_buf[k] = src[k*srcstep+i];
for( j = i; j <= size.width - 4; j += 4 )
{
double s0 = 0, s1 = 0, s2 = 0, s3 = 0;
const sT *tsrc = src + j;
for( k = 0; k < size.height; k++, tsrc += srcstep )
{
double a = col_buf[k];
s0 += a * tsrc[0]; // s0表示第0行
s1 += a * tsrc[1]; // s0表示第1行
s2 += a * tsrc[2]; // s0表示第2行
s3 += a * tsrc[3]; // s0表示第3行
}
tdst[j] = (dT)(s0*scale);
tdst[j+1] = (dT)(s1*scale);
tdst[j+2] = (dT)(s2*scale);
tdst[j+3] = (dT)(s3*scale);
}
for( ; j < size.width; j++ )
{
double s0 = 0;
const sT *tsrc = src + j;
for( k = 0; k < size.height; k++, tsrc += srcstep )
s0 += (double)col_buf[k] * tsrc[0];
tdst[j] = (dT)(s0*scale);
}
}
else
for( i = 0; i < size.width; i++, tdst += dststep )
{
// 计算每一列的(x[i] – avg[i])
if( !delta_buf )
for( k = 0; k < size.height; k++ )
col_buf[k] = src[k*srcstep+i] - delta[k*deltastep+i];
else
for( k = 0; k < size.height; k++ )
col_buf[k] = src[k*srcstep+i] - delta_buf[k*deltastep];
// 逐行处理，把前面得到的列向量转置后变成行（col_buf[k]）乘以列向量
// for循环处理过程中，每个循环处理4个列,直到处理完全部的列
for( j = i; j <= size.width - 4; j += 4 )
{
double s0 = 0, s1 = 0, s2 = 0, s3 = 0;
const sT *tsrc = src + j;
const dT *d = delta_buf ? delta_buf : delta + j;
for( k = 0; k < size.height; k++, tsrc+=srcstep, d+=deltastep )
{
double a = col_buf[k];
s0 += a * (tsrc[0] - d[0]); // s0表示第0行 s0 = C00*(S00-avg0)+ C01*(S10-avg0) + C02*(S20-avg0) + C03*(S30-avg0) + …
s1 += a * (tsrc[1] - d[1]); // s1表示第1行 …
s2 += a * (tsrc[2] - d[2]); // s2表示第2行 …
s3 += a * (tsrc[3] - d[3]); // s3表示第3行 s3 = C00*(S03-avg3)+ C01*(S13-avg3) + C02*(S23-avg3) + C03*(S33-avg3) + …
}
// 把计算结果存放在tdst中。
tdst[j] = (dT)(s0*scale);
tdst[j+1] = (dT)(s1*scale);
tdst[j+2] = (dT)(s2*scale);
tdst[j+3] = (dT)(s3*scale);
}
//计算上面没的计算到的位置（这里只计算对角线及右侧的元素的元素及），
// 从对角线的上元素开始，依左向右循环填充tdst
for( ; j < size.width; j++ )
{
double s0 = 0;
const sT *tsrc = src + j; // 第循环一次j+=1，相当于向右边移动一个向量（维度）
const dT *d = delta_buf ? delta_buf : delta + j;
for( k = 0; k < size.height; k++, tsrc+=srcstep, d+=deltastep )
s0 += (double)col_buf[k] * (tsrc[0] - d[0]);
tdst[j] = (dT)(s0*scale);
}
}
}
最后，要特别提示一下的是，如果你使用了IPP，那么其中有些运算就会使用Intel的IPP库函数，而不再是OpenCV中的源码。比如convertTo中，使用IPP，最终运行的是
CV_IPP_RUN_FAST(ipp_convertTo(src, dst, alpha, beta ));
如果不使用IPP，convertTo会调用
template<typename T, typename DT, typename WT> static void
cvtScale_( const T* src,size_t sstep, DT* dst, size_t dstep, Size size,WT scale, WT shift )并通过SIMD函数
cvtScale_SIMD<double, double, double>完成最后的运算。我这里使用的，是关闭IPP后的效果。关于如何关闭IPP，请参考
https://blog.csdn.net/tanmx219/article/details/81452343
参考资料
