OpenCV深入学习(8)--calcHist源码分析

      距离上篇深入学习又过了N久了，当时分析过calcHist的源码，不过有些地方不是很明白，后来一直忙着毕业，也没继续，但是心里总是觉得有个东西没有完成，这两天有点时间就把看的calcHist的源码整理了一下，贴出来，也算完了一个任务。

     calcHist的源码中根据计算的Mat的深度分了几种情况，分别调用不同的static函数实现的，其中8U的是用的普通函数，16U和32F的是用的模板，这里只对8U的分析了一下，就是在源码里加了一下注释。读读源码从中还真是学到不少东西。其中有些自己也不是很明白，错误之处还请指出，共同探讨。

      首先是一个总的函数calcHist

/*计算给定图像数组的直方图*/
void cv::calcHist( const Mat* images, int nimages, const int* channels,
InputArray _mask, OutputArray _hist, int dims, const int* histSize,
const float** ranges, bool uniform, bool accumulate )
{
/*从输入掩膜代理中获得掩膜矩阵*/
Mat mask = _mask.getMat();
/*确保直方图维数至少是1维，并且指定直方图bin大小的数组不是空的*/
CV_Assert(dims > 0 && histSize);

/*从输出直方图代理中获得指向存储输出直方图数据的数组的指针*/
uchar* histdata = _hist.getMat().data;
/*
创建输出直方图的存储空间，该代理根据其代理对象不同，调用被代理对象的创建函数，ex：是Mat的话，
调用Mat::create(int ndims,const int *sizes, int type)函数来分配存储空间
---分配的是32F类型的存储
*/
_hist.create(dims, histSize, CV_32F);
Mat hist = _hist.getMat(), ihist = hist;//从代理中取出或者将代理对象转换为Mat;

/*设置输出的flag标志位--注意最后的CV_32S，后面计算直方图使用的都是ihist临时矩阵，所以可以看出计算直方图时用的是32S深度的矩阵*/
ihist.flags = (ihist.flags & ~CV_MAT_TYPE_MASK)|CV_32S;
/*
非累加直方图或者开始获得的指向输出存储的指针与新分配的存储的指针不是同一个指针--？？，将输出都设为0;
否则，即,采用累加直方图，将hist转换到ihist--32S【用ihist作为临时存储？？】
*/
if( !accumulate || histdata != hist.data )
hist = Scalar(0.);
else
hist.convertTo(ihist, CV_32S);

/*为计算直方图时快速寻址分配的临时数组，保存一些参数*/
vector<uchar*> ptrs;
vector<int> deltas;
vector<double> uniranges;
Size imsize;

CV_Assert( !mask.data || mask.type() == CV_8UC1 );
/*调用本地静态函数计算后面所用的统计直方图时的参数---前8个参数为输入，后4个为输出，
第6个参数传入的是输出直方图的size数组，即矩阵各维度的大小
*/
histPrepareImages( images, nimages, channels, mask, dims, hist.size, ranges, uniform,
ptrs, deltas, imsize, uniranges );
/*是否均匀直方图，是的话用uniranges[0]指针初始化常量*/
const double* _uniranges = uniform ? &uniranges[0] : 0;

int depth = images[0].depth();//获取深度--所有输入图像的深度是一样的。

/*
根据深度不同，调用相应的函数计算直方图，深度为CV_8U时调用的是本地静态函数，
而其他两种情况下调用的是函数模板;
*/
if( depth == CV_8U )/*8U的图像的输出直方图值ihist是32S的*/
calcHist_8u(ptrs, deltas, imsize, ihist, dims, ranges, _uniranges, uniform );
else if( depth == CV_16U )
calcHist_<ushort>(ptrs, deltas, imsize, ihist, dims, ranges, _uniranges, uniform );
else if( depth == CV_32F )
calcHist_<float>(ptrs, deltas, imsize, ihist, dims, ranges, _uniranges, uniform );
else
CV_Error(CV_StsUnsupportedFormat, "");
/*将ihist在转换到hist中，并且将深度变换为32F--结合前面可见用于计算统计时使用的是32S深度的即整型格式，最后输出时用的是32F，应该是为了加快速度*/
ihist.convertTo(hist, CV_32F);
}

其中首先调用的是一个histPrepareImages函数：

/*
根据给定参数计算统计直方图所用的其他参数，即进行参数的转换；
前8个参数是输入的，后4个参数是输出；
---函数返回之后，ptrs中是计算输出直方图的dims从0~dims-1的每一维所用到的图像数据的首地址；
deltas中则是与之对应的图像的通道数以及字节偏移；imsize中是图像大小【如果是连续存储的则进行了顺序访问的加速处理】；
uniranges中是均匀直方图中的每个bin对应的灰度范围相关的东西---【没理解到底是怎么对应的？？】
*/
static void histPrepareImages( const Mat* images, int nimages, const int* channels,
const Mat& mask, int dims, const int* histSize,
const float** ranges, bool uniform,
vector<uchar*>& ptrs, vector<int>& deltas,
Size& imsize, vector<double>& uniranges )
{
int i, j, c;
/*
直方图的维数由channels或者nimages决定，如果channels为0，则图像数组的元素个数与dims必须一样，即数组中每一幅图像取一个通道；
如果channels不为0，则由channels和nimages共同决定dims：将images中的所有图像按照通道依次摊开排列，通道数从0开始，
channels数组指出有哪些通道用于计算直方图，而且channels指出的是按照直方图的第0维开始的，即直方图的第一维计算的是
channels[0]指定的通道图像。
*/
CV_Assert( channels != 0 || nimages == dims );
/*提取图像大小*/
imsize = images[0].size();
/*获取图像深度以及图像单通道的元素的大小*/
int depth = images[0].depth(), esz1 = (int)images[0].elemSize1();
bool isContinuous = true;
/*重置uchar数组的大小和int数组的大小*/
ptrs.resize(dims + 1);
deltas.resize((dims + 1)*2);

/*
对每一维进行处理---将图像数组按照通道摊开排列起来，将每一通道图像的数据的首地址放到ptrs[i]中；
deltas中保存的是每一维【或者说每一个通道图像】的映射到该通道所处的图像的通道数以及图像字节对齐产生的偏移量
*/
for( i = 0; i < dims; i++ )
{
/*通道参数为0---用维数指定通道--每一幅图像只取通道0*/
if(!channels)
{
j = i;
c = 0;//通道为0
CV_Assert( images[j].channels() == 1 );
}
else
{
/*通道数不是0--提取通道*/
c = channels[i];
CV_Assert( c >= 0 );
/*
按照channels参数的解释，将c设置为当前维【第i维】对应的通道图像
--这个通道图像或者是图像数组中的一幅单通道的图像，或者是数组中的一幅多通道图像的某个通道；
下面的for循环的目的是找出channels[i]对应的通道图像在输入数组的第几个中；循环结束之后
j是当前维即第i维对应的数组中的图像【单通道或者多通道】，c指向该图像中的一个通道；
*/
for( j = 0; j < nimages; c -= images[j].channels(), j++ )
if( c < images[j].channels() )
break;
CV_Assert( j < nimages );
}
/*确保每一幅图像大小深度都是一样的*/
CV_Assert( images[j].size() == imsize && images[j].depth() == depth );
/*图像是否是连续存储的--字节对齐*/
if( !images[j].isContinuous() )
isContinuous = false;

/*ptrs中存储每一个用于计算直方图的通道图像的数据首地址*/
ptrs[i] = images[j].data + c*esz1;
deltas[i*2] = images[j].channels();//每一幅用于计算直方图的图像的通道数；
/*
images[j].step/esz1---从每行占用的字节数步长计算图像每行占有的通道个数，
imsize.width*deltas[i*2]---直接通过宽度*每个元素的通道数计算得到的通道个数，
---所以这句应该是计算字节对齐导致的偏移--以深度类型为单位的【即8U,32F,32S这些】，对于没有字节对齐的是0---
*/
deltas[i*2+1] = (int)(images[j].step/esz1 - imsize.width*deltas[i*2]);
}
/*使用掩膜*/
if( mask.data )
{
CV_Assert( mask.size() == imsize && mask.channels() == 1 );
isContinuous = isContinuous && mask.isContinuous();
ptrs[dims] = mask.data; //最后一个指针放的是指向掩膜的数据；
deltas[dims*2] = 1;//通道数--掩膜为8UC1，所以是1；
//每一行的元素的个数
deltas[dims*2 + 1] = (int)(mask.step/mask.elemSize1());
}

/*连续存储，调整宽高，加快访问速度*/
if( isContinuous )
{
imsize.width *= imsize.height;
imsize.height = 1;
}

/*bin范围没有指定，默认对0-255范围的统计均匀直方图--必须是8U深度的图像*/
if( !ranges )
{
CV_Assert( depth == CV_8U );
/*调整bin均匀直方图的范围数组大小，并赋值*/
uniranges.resize( dims*2 );
for( i = 0; i < dims; i++ )
{
uniranges[i*2] = histSize[i]/256.;
uniranges[i*2+1] = 0;
}
}
else if( uniform )	/*采用均匀直方图，并且指定了每一维的bin的个数*/
{
uniranges.resize( dims*2 );
for( i = 0; i < dims; i++ )
{
CV_Assert( ranges[i] && ranges[i][0] < ranges[i][1] );
/*指定每一维的统计范围*/
double low = ranges[i][0], high = ranges[i][1];
double t = histSize[i]/(high - low);
uniranges[i*2] = t;
uniranges[i*2+1] = -t*low;
}
}
else
{
/*非均匀直方图，保证指定的范围是小->大*/
for( i = 0; i < dims; i++ )
{
size_t j, n = histSize[i];
for( j = 0; j < n; j++ )
CV_Assert( ranges[i][j] < ranges[i][j+1] );
}
}
}

然后calcHist中就调用8U深度图像的统计函数calcHist_8u进行具体的统计计算过程

/*
8U深度的图像的直方图统计函数
首先生成了映射查找表加快处理速度---理解映射的对应关系是关键--对应关系需要对几个输入参数ptr，deltas以及最重要的查找表tab的存储结构搞清楚；
后面具体的统计计算过程，对dims>=2的过程基本一样，对dims=1的稍有不同，需要先搞清楚映射对应关系之后，这个统计计算过程才可以比较清晰。
---参数ptrs; deltas; imsize;都是上面通过histPrepareImages函数计算获得的，_uniranges指向的是histPrepareImages计算获得的uniranges数组的第一个元素；
*/
static void
calcHist_8u( vector<uchar*>& _ptrs, const vector<int>& _deltas,
Size imsize, Mat& hist, int dims, const float** _ranges,
const double* _uniranges, bool uniform )
{
/*指向用于计算直方图的图像数据的指针的指针*/
uchar** ptrs = &_ptrs[0];
/*对应上面的计算图像通道数以及偏移的数组的指针*/
const int* deltas = &_deltas[0];
uchar* H = hist.data;//输出直方图的数据指针--32S
int i, x;
const uchar* mask = _ptrs[dims];//掩膜的数据指针
int mstep = _deltas[dims*2 + 1];//掩膜的数据指针移动步长
vector<size_t> _tab;

/*生成映射查找表-----要弄清楚映射规则需要弄清楚几个数组或者矩阵的存储结构：Mat即hist；_tab；ptr；delta；_uniranges等*/
calcHistLookupTables_8u( hist, SparseMat(), dims, _ranges, _uniranges, uniform, false, _tab );
const size_t* tab = &_tab[0];
/*分不同的维不同对待*/
if( dims == 1 )
{
/*d0--通道数，step0字节对齐偏移*/
int d0 = deltas[0], step0 = deltas[1];
int matH[256] = {0};//对各个灰度值的统计暂存
const uchar* p0 = (const uchar*)ptrs[0];//图像数据指针
/**/
for( ; imsize.height--; p0 += step0, mask += mstep )
{
/*不采用掩膜*/
if( !mask )
{
/*单通道图像*/
if( d0 == 1 )
{
/*单通道，4字节对齐，对每个字节都进行统计*/
for( x = 0; x <= imsize.width - 4; x += 4 )
{
int t0 = p0[x], t1 = p0[x+1];
matH[t0]++; matH[t1]++;
t0 = p0[x+2]; t1 = p0[x+3];
matH[t0]++; matH[t1]++;
}
p0 += x;
}
else   /*多通道图像*/
{
/*多通道，只对第0个通道进行统计，通道数d0做偏移*/
for( x = 0; x <= imsize.width - 4; x += 4 )
{
int t0 = p0[0], t1 = p0[d0];
matH[t0]++; matH[t1]++;
p0 += d0*2;
t0 = p0[0]; t1 = p0[d0];
matH[t0]++; matH[t1]++;
p0 += d0*2;
}
}
/*非4字节的对齐，自己单独处理*/
for( ; x < imsize.width; x++, p0 += d0 )
matH[*p0]++;
}
else /*采用掩膜时直接统计*/
for( x = 0; x < imsize.width; x++, p0 += d0 )
if( mask[x] )
matH[*p0]++;
}
/*将统计得到的数据存入直方图矩阵中，采用了灰度值到bin的映射查找到对应的bin位置*/
for( i = 0; i < 256; i++ )
{
size_t hidx = tab[i];//hidx--即直方图中的位置
//只统计不超限的数值
if( hidx < OUT_OF_RANGE )
*(int *)(H + hidx) += matH[i];
}
}
else if( dims == 2 )/*2维直方图*/
{
/*取出用于计算两维直方图的图像的通道数，字节对齐偏移，数据指针等*/
int d0 = deltas[0], step0 = deltas[1],
d1 = deltas[2], step1 = deltas[3];
const uchar* p0 = (const uchar*)ptrs[0];
const uchar* p1 = (const uchar*)ptrs[1];
/**/
for( ; imsize.height--; p0 += step0, p1 += step1, mask += mstep )
{
if( !mask )/*无掩膜*/
{
/*首先根据映射表计算出在直方图中的位置，然后累加*/
for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1 )
{
/* idx 为第0维的偏移【tab[*p0]】+第1维的偏移【tab[*p1 + 256]】---类似于图像的行偏移+列偏移*/
size_t idx = tab[*p0] + tab[*p1 + 256];
if( idx < OUT_OF_RANGE )
++*(int*)(H + idx);
}
}
else/*使用掩膜时*/
{
/*跟上面的一样只是多了掩膜判断*/
for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1 )
{
size_t idx;
if( mask[x] && (idx = tab[*p0] + tab[*p1 + 256]) < OUT_OF_RANGE )
++*(int*)(H + idx);
}
}
}
}
else if( dims == 3 )/*3维直方图--与2维的类似，只是多了一维，寻找直方图位置的idx的计算复杂了*/
{
int d0 = deltas[0], step0 = deltas[1],
d1 = deltas[2], step1 = deltas[3],
d2 = deltas[4], step2 = deltas[5];

const uchar* p0 = (const uchar*)ptrs[0];
const uchar* p1 = (const uchar*)ptrs[1];
const uchar* p2 = (const uchar*)ptrs[2];
/**/
for( ; imsize.height--; p0 += step0, p1 += step1, p2 += step2, mask += mstep )
{
if( !mask )/**/
for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1, p2 += d2 )
{
/*第0+1+2维的偏移获得总的偏移量*/
size_t idx = tab[*p0] + tab[*p1 + 256] + tab[*p2 + 512];
if( idx < OUT_OF_RANGE )
++*(int*)(H + idx);
}
else/**/
for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1, p2 += d2 )
{
size_t idx;
if( mask[x] && (idx = tab[*p0] + tab[*p1 + 256] + tab[*p2 + 512]) < OUT_OF_RANGE )
++*(int*)(H + idx);
}
}
}
else/*高维直方图(>3维)*/
{
for( ; imsize.height--; mask += mstep )
{
if( !mask )/*无掩膜图像计算*/
for( x = 0; x < imsize.width; x++ )
{
uchar* Hptr = H;
for( i = 0; i < dims; i++ )
{
size_t idx = tab[*ptrs[i] + i*256];
if( idx >= OUT_OF_RANGE )
break;
Hptr += idx;
ptrs[i] += deltas[i*2];
}
/*判断是否有超限的--即前面设置的在ranges范围外的设置为了OUT_OF_RANGE，统计数据时需要跨过这些值*/
if( i == dims )
++*((int*)Hptr);
else
for( ; i < dims; i++ )
ptrs[i] += deltas[i*2];
}
else/*使用掩模图像*/
for( x = 0; x < imsize.width; x++ )
{
uchar* Hptr = H;
int i = 0;
if( mask[x] )
for( ; i < dims; i++ )
{
size_t idx = tab[*ptrs[i] + i*256];
if( idx >= OUT_OF_RANGE )
break;
Hptr += idx;
ptrs[i] += deltas[i*2];
}

if( i == dims )
++*((int*)Hptr);
else
for( ; i < dims; i++ )
ptrs[i] += deltas[i*2];
}
/*字节对齐偏移*/
for( i = 0; i < dims; i++ )
ptrs[i] += deltas[i*2 + 1];
}
}
}

 

在calcHist_8u中使用了一个查找表用于加快处理速度，生成查找表的函数如下

static const size_t OUT_OF_RANGE = (size_t)1 << (sizeof(size_t)*8 - 2);

/*
建立查找表--即,将8u每个灰度值映射到相应的bin值上加快统计直方图的计算速度，
_tab中是从直方图的每一维0-255的值到输出直方图中的位置的映射，即tab中的index是第0维的0~255，第1维的0~255...第dims-1维的0~255；
而value则是每一个0~255的值映射到最终直方图中的第i维的偏移【即相对于该维数据首地址的偏移量】。不在要求范围之内的设为OUT_OF_RANGE
*/
static void
calcHistLookupTables_8u( const Mat& hist, const SparseMat& shist,
int dims, const float** ranges, const double* uniranges,
bool uniform, bool issparse, vector<size_t>& _tab )
{
/*上下限？*/
const int low = 0, high = 256;
int i, j;
/*调整大小*/
_tab.resize((high-low)*dims);
size_t* tab = &_tab[0];
/*均匀直方图*/
if( uniform )
{
/**/
for( i = 0; i < dims; i++ )
{
/*
uniranges的值在histPrepareImages函数中
double _low = ranges[i][0], _high = ranges[i][1];
double t = histSize[i]/(_high - _low);
uniranges[i*2] = t;
uniranges[i*2+1] = -t*_low;
*/
double a = uniranges[i*2];
double b = uniranges[i*2+1];
/*是否是稀疏直方图*/
int sz = !issparse ? hist.size[i] : shist.size(i);
size_t step = !issparse ? hist.step[i] : 1; //step中的是输出直方图的第i维的步长
/*给tab赋值*/
for( j = low; j < high; j++ )
{
/*
--计算low-high范围内的每个值被划分到bin中哪一个
idx=cvFloor(j*a+b)=cvFloor(j*t-t*_low)=cvFloor(t*(j-_low))=cvFloor(histSize[i]*(j-_low)/(_high-_low));
--将a，b的值进行了替换，a，b的值见上面；而histSize[i]中的是给定的第i维的bin的个数，low是8u的下限0，_high和_low则是用户给出的range范围。
*/
int idx = cvFloor(j*a + b);
size_t written_idx;
/*
超出大小范围
--没有超范围，要写入的索引是idx与步长的乘积？？
--超出范围，设置为次高位为1，其他位为0的size_t的数。
*/
if( (unsigned)idx < (unsigned)sz )
written_idx = idx*step;
else
written_idx = OUT_OF_RANGE;
/*给tab赋值--tab中的idx为每一维的灰度值，value为该灰度值映射到输出直方图中的第i维的偏移--*/
tab[i*(high - low) + j - low] = written_idx;
}
}
}
else
{
/*非均匀直方图*/
for( i = 0; i < dims; i++ )
{
/*取得ranges下限？*/
int limit = std::min(cvCeil(ranges[i][0]), high);
int idx = -1, sz = !issparse ? hist.size[i] : shist.size(i);
size_t written_idx = OUT_OF_RANGE;
size_t step = !issparse ? hist.step[i] : 1;

/*
根据指定的ranges各段的上下限设置非均匀直方图的映射表；
对于[low,ranges_low)和(ranges_high,high)范围的值映射为超限，
中间的值每段映射到相应的值。采用对整个范围按索引小->大循环，循环中
再采用对分段范围循环赋值，然后在大循环中修改小循环的循环跳出值的方法
------循环一遍即可遍历整个范围，方法很好！！
*/
for(j = low;;)
{
/*对于从low到limit范围的值映射为指定的written_idx*/
for( ; j < limit; j++ )
tab[i*(high - low) + j - low] = written_idx;
/*
取得ranges的下一个端点，并将其设为下一个的bin的上限，
并且计算该段数值对应的映射到输出直方图中的偏移位置。
*/
if( (unsigned)(++idx) < (unsigned)sz )
{
limit = std::min(cvCeil(ranges[i][idx+1]), high);
written_idx = idx*step;
}
else
{
/*从指定范围上限到high的值设为超限--退出循环*/
for( ; j < high; j++ )
tab[i*(high - low) + j - low] = OUT_OF_RANGE;
break;
}
}
}
}
}

整个对8U深度的图像的直方图统计就涉及到这几个函数，以上注释的错误之处还望不吝赐教！

calcHist的学习就到此为止。