OpenCV 直方图计算函数 calcHist源码深度剖析
OpenCV中的直方图计算函数calcHist函数
·
OpenCV中的直方图计算函数calcHist函数可以计算给定的若干幅图像的指定的通道的统计直方图！
calcHist函数原型为
//!计算给定图像集合的联合密度直方图 （joint dense histogram）
CV_EXPORTS void calcHist( const Mat* images, int nimages, const int* channels, InputArray mask,
OutputArray hist, int dims, const int* histSize,
const float** ranges, bool uniform=true, bool accumulate=false );
该函数内部的实现方式是怎样的呢？
首先调用了函数 histPrepareImages(....）准备一些中间变量，
//该函数主要初始化了三个中间临时变量 ptrs，deltas，uniranges
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 )
接着根据输入图像的位深度调用不同的函数实现
针对8位深度的图像，调用了函数calcHist_8u(.....)
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 )
针对16位和32位深度图像，调用了模板类函数
template<typename T> static void
calcHist_( vector<uchar*>& _ptrs, const vector<int>& _deltas, Size imsize, Mat& hist, int dims, const float** _ranges, const double* _uniranges, bool uniform )
调用方式如下：
int depth = images[0].depth();//输入图像的深度
//根据输入图像的位深度调用不同的直方图计算函数
if( depth == CV_8U ) //8位深度
calcHist_8u(ptrs, deltas, imsize, ihist, dims, ranges, _uniranges, uniform );
else if( depth == CV_16U )//16位深度
calcHist_<ushort>(ptrs, deltas, imsize, ihist, dims, ranges, _uniranges, uniform );
else if( depth == CV_32F ) //32位深度
calcHist_<float>(ptrs, deltas, imsize, ihist, dims, ranges, _uniranges, uniform );
else
CV_Error(CV_StsUnsupportedFormat, "");
为了研究源码，我把位于...\OpenCV245\modules\imgproc\src\文件夹下的histogram.cpp文件下的相关函数加到自定义命名空间中，这么做就可以调试运行每一行代码啦
在命名空间MyHistFunc中定义了一系列与直方图的计算为主要目的之函数，所有在该Histogram.cpp文件中被定义为 static 类型的函数，在头文件Histogram.h中都没有进行声明，所以这些静态函数都只能在该文件中被该文件的定义在他之后的其他函数调用；而只有函数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 )在头文件中进行了声明，所以该函数是唯一暴露在外边，可以被其他文件中的函数（如Main()函数）调用的函数。 如果你想在别的文件的函数中调用该文件中的static类型的函数，那么只需要把该函数的声明复制到头文件Histogram.h中，就像calcHist(....)函数那样。
histogram.h 头文件的内容（这个头文件是在新建的工程中自己写的）：
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/imgproc/imgproc_c.h"
#include "opencv2/core/internal.hpp"
#include
#include
#include
#include
#include
#include
#include
using namespace std;
using namespace cv;
//自己定义了一个MyHistFunc的命名空间，然后把cv命名空间下的对应函数拷贝出来，
//以免编译器不知道该使用哪个命名空间中的calcHist(。。。）函数，在主函数调用时
//要加上命名空间的前缀，用法： MyHistFunc:: calcHist(....)
namespace MyHistFunc
{
//! 计算给定的一系列图像的联合密度分布直方图（joint dense histogram）
void 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 );
};
#include "stdafx.h"
#include "Histogram.h"
/**
在命名空间MyHistFunc中定义了一系列与直方图的计算为主要目的之函数，所有在该Histogram.cpp文件中
被定义为 static 类型的函数，在头文件Histogram.h中都没有进行声明，所以这些静态函数都只能在该文件中
被该文件的定义在他之后的其他函数调用；而只有函数
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 )
在头文件中进行了声明，所以该函数是唯一暴露在外边，可以被其他文件中的函数（如Main()函数）调用
的函数。 如果你想在别的文件的函数中调用该文件中的static类型的函数，那么只需要把该函数的声明
复制到头文件Histogram.h中，就像calcHist(....)函数那样。
*/
namespace MyHistFunc
{
//#define HAVE_TBB 1
////////////////// Helper functions //////////////////////
static const size_t OUT_OF_RANGE = (size_t)1 << (sizeof(size_t)*8 - 2);
//不管有多少通道，对于8位深度的图像，单一通道上的每个像素点的取值范围是从[0,256),
//所以，在计算8位深度的图像的直方图的时候，总是假定每一个维度上的取值范围都是[0,256)
//然后根据参数中给出的每一维度的范围制作查找表；在每一维上，查找表的长度都是256，
//查找表中第i维第j的位置上存放的就是灰度值等于j的像素点所对应的i维度上的bin的编号
static void
calcHistLookupTables_8u( const Mat& hist, const SparseMat& shist,
int dims, const float** ranges, const double* uniranges,
bool uniform, bool issparse, vector& _tab )
{
const int low = 0, high = 256;//位深度为8位的图像，一个基础元素的范围都在[0,255)内
int i, j;
//为查找表开辟空间，查找表长度 = 直方图维数*256
//每一个维度所占的存储位置都是256个
_tab.resize((high-low)*dims);
size_t* tab = &_tab[0];//指向查找表数组头的指针
if( uniform ) //均匀直方图，也就是同一维度内部的各个bin的区间宽度是一样的
{ //但是不同的维度上的bin的区间宽度不一定是一样的
for( i = 0; i < dims; i++ )
{
double a = uniranges[i*2];//第i维上单位区间内bin的上限
double b = uniranges[i*2+1];
int sz = !issparse ? hist.size[i] : shist.size(i);//第i维上可以容纳的bin数目
//step[i]是第i维的步长(跨度)，多维矩阵每一维度的步长不一样的
//对于单通道32位深度的浮点数矩阵，第0维的指针移动步长是sizeof(float)即==4字节长度
//第1维的指针移动步长是 cols*sizeof(float) = 4*hist.cols
size_t step = !issparse ? hist.step[i] : 1;
//下面这一循环就是把第i维上的每一个灰度值(0,255)压缩到参数规定的
//区间{ranges[i][low]到ranges[i][high]}内，然后在对应位置上保存bin索引值
for( j = low; j < high; j++ )// 灰度值j 从0到255 循环
{
//灰度值j在输入参数规定的第i维灰度区间内的索引值
int idx = cvFloor(j*a + b);//该索引值是数学上线性运算后得到的
size_t written_idx; //该索引是写入直方图矩阵时的索引值，
//由于查找表是单行矩阵，而直方图是多维矩阵，所以索引idx要乘以直方图
//每一维对应的跨度(step)才能与查找表的具体位置对应
if( (unsigned)idx < (unsigned)sz )
written_idx = idx*step;//索引idx要乘以直方图每一维对应的跨度(step)
else
written_idx = OUT_OF_RANGE;//对直方图的索引项不能超出直方图的总长度
//对查找表第i维第j个位置的索引是256*i+j,该位置保存的是当灰度值为j时，应该访问的
//直方图的对应位置的索引
tab[i*(high - low) + j - low] = written_idx; //相当于tab[i*256+j] = written_idx;
}
}
}
else //非均匀直方图的情形，同一维度内各个bin的区间宽度不一样，从最小值到最大值非均匀划分
{
for( i = 0; i < dims; i++ )
{
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;
for(j = low;;)
{
//低于第i维的最小区间的下界的那些灰度值
for( ; j < limit; j++ )
tab[i*(high - low) + j - low] = written_idx;
//在第i维指定区间内的灰度-索引值查找表
if( (unsigned)(++idx) < (unsigned)sz )
{
limit = std::min(cvCeil(ranges[i][idx+1]), high);
written_idx = idx*step;
}
else
{
//超出了第i维度规定的区间中的最大上界的那些灰度值 j
for( ; j < high; j++ )
tab[i*(high - low) + j - low] = OUT_OF_RANGE;
break;
}
}
}
}
}
//该函数主要初始化了三个中间临时变量 ptrs，deltas，uniranges
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& ptrs, vector& deltas,
Size& imsize, vector& uniranges )
{
int i, j, c;
CV_Assert( channels != 0 || nimages == dims );
imsize = images[0].size();//images是个图像数组，所有图像的尺寸必须一样的
//输入图像的位深度，elemSize1 表示的是矩阵中数据类型的大小，即 elemSize/channels的大小
int depth = images[0].depth(), esz1 = (int)images[0].elemSize1();//对于CV_8UCx类的图像，==1
bool isContinuous = true; //数据在内存中是否连续存放
//dims是要构建的直方图的维数,+1是因为ptrs中最后一个位置要留给掩膜矩阵mask呀
ptrs.resize(dims + 1);//指向对应的图像数组首地址，数组中最后一个指针指向mask
deltas.resize((dims + 1)*2);// +1 是因为deltas中最后两个位置要留给mask呀
for( i = 0; i < dims; i++ )
{
if(!channels)
{//channels为空的时候，该句将执行
c = 0;
CV_Assert( images[i].channels() == 1 );
CV_Assert( images[i].size() == imsize && images[i].depth() == depth );
if( !images[i].isContinuous() )
isContinuous = false;
ptrs[i] = images[i].data + c*esz1; //ptrs指针指向对应的图像数据数组首地址，这时c==0
deltas[i*2] = images[i].channels();//第j个图像的通道数目，指针ptrs[i]在图像中同一行内移动的步长
deltas[i*2+1] = (int)(images[i].step/esz1 - imsize.width*deltas[i*2]);
}
else
{
c = channels[i];
CV_Assert( c >= 0 );
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[i] = images[j].data + c*esz1; //ptrs指针指向对应的图像矩阵c通道的首地址
deltas[i*2] = images[j].channels();//第j个图像的通道数目
deltas[i*2+1] = (int)(images[j].step/esz1 - imsize.width*deltas[i*2]);
}
}
//如果mask矩阵不为空，则初始化指向mask矩阵的数据的指针和步长
if( mask.data )
{
CV_Assert( mask.size() == imsize && mask.channels() == 1 );
isContinuous = isContinuous && mask.isContinuous();
ptrs[dims] = mask.data;
deltas[dims*2] = 1; //列与列之间的步长，掩膜矩阵mask是单通道的，所以 ==1
//mask.step/mask.elemSize1()是mask矩阵一行的字节数，
deltas[dims*2 + 1] = (int)(mask.step/mask.elemSize1());//指针在行与行之间的步长
}
#ifndef HAVE_TBB
if( isContinuous )
{
imsize.width *= imsize.height;
imsize.height = 1;
}
#endif
//ranges[][]是个具有dims行2列的数组，其中第i行ranges[i]代表了直方图第i维的取值区间
//该区间包括下界ranges[i][0]和上界ranges[i][1]
if( !ranges ) //ranges数组为空的话，说明参数中没有指定每一维的取值区间，
{ //则默认所有维度的取值区间是（0,255），而且默认把区间[0,255]
//划分成histSize个宽度一样的 bin
CV_Assert( depth == CV_8U );
uniranges.resize( dims*2 );
for( i = 0; i < dims; i++ )
{/*第i维上单位区间内直方图bin的数目*/
uniranges[i*2] = histSize[i]/256.;
uniranges[i*2+1] = 0;
}
}
else if( uniform )//如果uniform为真，则说明直方图的每一维上的区间是均匀划分的
{ //也就是直方图每一维的所有bin的宽度是一样的
uniranges.resize( dims*2 );
for( i = 0; i < dims; i++ )
{
//保证第i维的取值范围区间的合理性（a,b）满足a static void
calcHist_( vector& _ptrs, const vector& _deltas,
Size imsize, Mat& hist, int dims, const float** _ranges,
const double* _uniranges, bool uniform )
{
T** ptrs = (T**)&_ptrs[0];
const int* deltas = &_deltas[0];
uchar* H = hist.data;
int i, x;
const uchar* mask = _ptrs[dims];
int mstep = _deltas[dims*2 + 1];
int size[CV_MAX_DIM];
size_t hstep[CV_MAX_DIM];
for( i = 0; i < dims; i++ )
{
size[i] = hist.size[i];
hstep[i] = hist.step[i];
}
if( uniform )
{
const double* uniranges = &_uniranges[0];
if( dims == 1 )
{
#ifdef HAVE_TBB
calcHist1D_Invoker body(_ptrs, _deltas, hist, _uniranges, size[0], dims, imsize);
parallel_for(BlockedRange(0, imsize.height), body);
return;
#endif
double a = uniranges[0], b = uniranges[1];
int sz = size[0], d0 = deltas[0], step0 = deltas[1];
const T* p0 = (const T*)ptrs[0];
for( ; imsize.height--; p0 += step0, mask += mstep )
{
if( !mask )
for( x = 0; x < imsize.width; x++, p0 += d0 )
{
int idx = cvFloor(*p0*a + b);
if( (unsigned)idx < (unsigned)sz )
((int*)H)[idx]++;
}
else
for( x = 0; x < imsize.width; x++, p0 += d0 )
if( mask[x] )
{
int idx = cvFloor(*p0*a + b);
if( (unsigned)idx < (unsigned)sz )
((int*)H)[idx]++;
}
}
}
else if( dims == 2 )
{
#ifdef HAVE_TBB
calcHist2D_Invoker body(_ptrs, _deltas, hist, _uniranges, size, dims, imsize, hstep);
parallel_for(BlockedRange(0, imsize.height), body);
return;
#endif
double a0 = uniranges[0], b0 = uniranges[1], a1 = uniranges[2], b1 = uniranges[3];
int sz0 = size[0], sz1 = size[1];
int d0 = deltas[0], step0 = deltas[1],
d1 = deltas[2], step1 = deltas[3];
size_t hstep0 = hstep[0];
const T* p0 = (const T*)ptrs[0];
const T* p1 = (const T*)ptrs[1];
for( ; imsize.height--; p0 += step0, p1 += step1, mask += mstep )
{
if( !mask )
for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1 )
{
int idx0 = cvFloor(*p0*a0 + b0);
int idx1 = cvFloor(*p1*a1 + b1);
if( (unsigned)idx0 < (unsigned)sz0 && (unsigned)idx1 < (unsigned)sz1 )
((int*)(H + hstep0*idx0))[idx1]++;
}
else
for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1 )
if( mask[x] )
{
int idx0 = cvFloor(*p0*a0 + b0);
int idx1 = cvFloor(*p1*a1 + b1);
if( (unsigned)idx0 < (unsigned)sz0 && (unsigned)idx1 < (unsigned)sz1 )
((int*)(H + hstep0*idx0))[idx1]++;
}
}
}
else if( dims == 3 )
{
#ifdef HAVE_TBB
if( calcHist3D_Invoker::isFit(hist, imsize) )
{
calcHist3D_Invoker body(_ptrs, _deltas, imsize, hist, uniranges, dims, hstep, size);
parallel_for(BlockedRange(0, imsize.height), body);
return;
}
#endif
double a0 = uniranges[0], b0 = uniranges[1],
a1 = uniranges[2], b1 = uniranges[3],
a2 = uniranges[4], b2 = uniranges[5];
int sz0 = size[0], sz1 = size[1], sz2 = size[2];
int d0 = deltas[0], step0 = deltas[1],
d1 = deltas[2], step1 = deltas[3],
d2 = deltas[4], step2 = deltas[5];
size_t hstep0 = hstep[0], hstep1 = hstep[1];
const T* p0 = (const T*)ptrs[0];
const T* p1 = (const T*)ptrs[1];
const T* p2 = (const T*)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 )
{
int idx0 = cvFloor(*p0*a0 + b0);
int idx1 = cvFloor(*p1*a1 + b1);
int idx2 = cvFloor(*p2*a2 + b2);
if( (unsigned)idx0 < (unsigned)sz0 &&
(unsigned)idx1 < (unsigned)sz1 &&
(unsigned)idx2 < (unsigned)sz2 )
((int*)(H + hstep0*idx0 + hstep1*idx1))[idx2]++;
}
else
for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1, p2 += d2 )
if( mask[x] )
{
int idx0 = cvFloor(*p0*a0 + b0);
int idx1 = cvFloor(*p1*a1 + b1);
int idx2 = cvFloor(*p2*a2 + b2);
if( (unsigned)idx0 < (unsigned)sz0 &&
(unsigned)idx1 < (unsigned)sz1 &&
(unsigned)idx2 < (unsigned)sz2 )
((int*)(H + hstep0*idx0 + hstep1*idx1))[idx2]++;
}
}
}
else
{
for( ; imsize.height--; mask += mstep )
{
if( !mask )
for( x = 0; x < imsize.width; x++ )
{
uchar* Hptr = H;
for( i = 0; i < dims; i++ )
{
int idx = cvFloor(*ptrs[i]*uniranges[i*2] + uniranges[i*2+1]);
if( (unsigned)idx >= (unsigned)size[i] )
break;
ptrs[i] += deltas[i*2];
Hptr += idx*hstep[i];
}
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;
i = 0;
if( mask[x] )
for( ; i < dims; i++ )
{
int idx = cvFloor(*ptrs[i]*uniranges[i*2] + uniranges[i*2+1]);
if( (unsigned)idx >= (unsigned)size[i] )
break;
ptrs[i] += deltas[i*2];
Hptr += idx*hstep[i];
}
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];
}
}
}
else
{
// non-uniform histogram
const float* ranges[CV_MAX_DIM];
for( i = 0; i < dims; i++ )
ranges[i] = &_ranges[i][0];
for( ; imsize.height--; mask += mstep )
{
for( x = 0; x < imsize.width; x++ )
{
uchar* Hptr = H;
i = 0;
if( !mask || mask[x] )
for( ; i < dims; i++ )
{
float v = (float)*ptrs[i];
const float* R = ranges[i];
int idx = -1, sz = size[i];
while( v >= R[idx+1] && ++idx < sz )
; // nop
if( (unsigned)idx >= (unsigned)sz )
break;
ptrs[i] += deltas[i*2];
Hptr += idx*hstep[i];
}
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];
}
}
}
static void
calcHist_8u( vector& _ptrs, const vector& _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; //指向直方图数据的指针
int x;
const uchar* mask = _ptrs[dims];//指向掩膜矩阵的指针
int mstep = _deltas[dims*2 + 1];//mask step ,掩膜矩阵行与行之间的步长
vector _tab; //计算直方图时的查找表
//计算直方图查找表
calcHistLookupTables_8u( hist, SparseMat(), dims, _ranges, _uniranges, uniform, false, _tab );
const size_t* tab = &_tab[0]; //指向直方图查找表数组头的指针
if( dims == 1 )//一维直方图的情形
{
#ifdef HAVE_TBB
if( CalcHist1D_8uInvoker::isFit(hist, imsize) )
{
int treadsNumber = tbb::task_scheduler_init::default_num_threads();
int grainSize = imsize.height/treadsNumber;
tbb::mutex histogramWriteLock;
CalcHist1D_8uInvoker body(_ptrs, _deltas, imsize, hist, dims, _tab, &histogramWriteLock);
parallel_for(BlockedRange(0, imsize.height, grainSize), body);
return;
}
#endif
//deltas[]应该是指针每次移动的步子
//deltas[0]，d0 是指针横向移动的步长，在一行内每次移动的步长等于通道数
//deltas[1],step0 是纵向移动的步长，在行与行之间移动
int d0 = deltas[0] ; //其实就是输入图像的通道数
int 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 )//对应于单通道图像
{
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 //对应于多通道图像
{
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;
}
}
//上面的if-else语句执行之后，x的值已经==imsize.width了啊
//所以下面的循环好像不会被执行的
for( ; x < imsize.width; x++, p0 += d0 )
{
matH[*p0]++;
}
}
else //如果该行掩膜矩阵的指针mask不为空，则表示该行内有被禁止统计的
{ //像素位置，则该行内的每一个位置在统计前要判断是否被掩盖了，即
//mask[x]等于0的位置不予统计
for( x = 0; x < imsize.width; x++, p0 += d0 )
if( mask[x] )
matH[*p0]++;
}
}
//将matH中的统计信息根据查找表放入直方图矩阵中
for(int i = 0; i < 256; i++ )
{
//从查找表取出灰度值为i的像素点对应的直方图存储位置索引
size_t hidx = tab[i];//一维直方图的情形，所以查找表中有256项，有可能多个i对应于一个hidx
if( hidx < OUT_OF_RANGE )
*(int*)(H + hidx) += matH[i];//(int*)(H + hidx)定位到直方图的第hidx个bin处
}
}
else if( dims == 2 )
{
#ifdef HAVE_TBB
if( CalcHist2D_8uInvoker::isFit(hist, imsize) )
{
callCalcHist2D_8u(_ptrs, _deltas, imsize, hist, dims, _tab);
return;
}
#endif
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 )
{
size_t idx = tab[*p0] + tab[*p1 + 256];
if( idx < OUT_OF_RANGE )
++*(int*)(H + idx); //二维直方图没有用matH做统计，而是在原来的H上直接累加
}
}
else
{ //有掩膜矩阵的情况
for( x = 0; x < imsize.width; x++, p0 += d0, p1 += d1 )
{
if(mask[x]) //mask[x]等于0的位置不予统计
{
//从查找表取出*p0和*p1对应的直方图存储位置索引
size_t idx = tab[*p0] + tab[*p1 + 256] ;//p0是第一维，p1第二维
if( idx < OUT_OF_RANGE )
++*(int*)(H + idx); //在对应位置上累加直方图
}
}
}
}
}
else if( dims == 3 )
{
#ifdef HAVE_TBB
if( CalcHist3D_8uInvoker::isFit(hist, imsize) )
{
callCalcHist3D_8u(_ptrs, _deltas, imsize, hist, dims, _tab);
return;
}
#endif
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 )
{
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
{
for( ; imsize.height--; mask += mstep )
{
if( !mask )
for( x = 0; x < imsize.width; x++ )
{
uchar* Hptr = H;
int i = 0;
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];
}
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(int i = 0; i < dims; i++ )
ptrs[i] += deltas[i*2 + 1];
}
}
}
}
//! 计算给定的一系列图像的联合密度分布直方图（joint dense histogram）
//该函数按照直方图的数据类型和直方图的维数分情况计算直方图
void MyHistFunc::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(); //把InputArray类型的参数转换为Mat类型
//CV_Assert(...)语句只有Debug模式下会执行，release模式下不执行
CV_Assert(dims > 0 && histSize);//参数检验，
//参数_hist是输出参数，保存着计算好的直方图
uchar* histdata = _hist.getMat().data; //指向直方图的数据部分的指针
_hist.create(dims, histSize, CV_32F); //为输出直方图开辟空间
Mat hist = _hist.getMat(), ihist = hist; //把OutputArray类型的参数转换为Mat类型
ihist.flags = (ihist.flags & ~CV_MAT_TYPE_MASK)|CV_32S;//直方图矩阵的类型信息
//如果参数_hist中本来就有之前计算的直方图，判断是否重新计算或累加上去
if( !accumulate || histdata != hist.data )
hist = Scalar(0.); //清空原来的数据，当前直方图不累加到原来的直方图上去
else
hist.convertTo(ihist, CV_32S);//保留原来的直方图，当前计算的结果累加到原来的直方图
//定义一些中间临时变量
vector ptrs; //指向图像数据和掩膜矩阵数据的指针数组，有dims+1个指针
vector deltas; //存放每一维上的指针移动步长
vector uniranges;//直方图每一个维度上单位区间内的直方图数目
Size imsize; //图像尺寸
//检查掩膜矩阵的有效性，必须是8位单通道
CV_Assert( !mask.data || mask.type() == CV_8UC1 );
//前期准备工作，主要是填充前面声明的辅助变量ptrs, deltas, imsize, uniranges
histPrepareImages( images, nimages, channels, mask, dims, hist.size, ranges,
uniform, ptrs, deltas, imsize, uniranges );
//判断在每一维度上的划分是否刻度均匀，
const double* _uniranges = uniform ? &uniranges[0] : 0;
int depth = images[0].depth();//输入图像的深度
//根据输入图像的位深度调用不同的直方图计算函数
if( depth == CV_8U ) //8位深度
calcHist_8u(ptrs, deltas, imsize, ihist, dims, ranges, _uniranges, uniform );
else if( depth == CV_16U )//16位深度
calcHist_(ptrs, deltas, imsize, ihist, dims, ranges, _uniranges, uniform );
else if( depth == CV_32F ) //32位深度
calcHist_(ptrs, deltas, imsize, ihist, dims, ranges, _uniranges, uniform );
else
CV_Error(CV_StsUnsupportedFormat, "");
ihist.convertTo(hist, CV_32F);
}
主函数的实现：（主程序中调用了命名空间MyHistFunc中的函数calcHist(....)）
//
#include "stdafx.h"
#include
#include
#include "Histogram.h"
using namespace std;
using namespace cv;
Mat image; //母图像
Mat roi_hist; //兴趣区域图像块的直方图
Mat roi_patch; //兴趣区域的图像块
bool selectObject = false; //selectObject的初始值为false,
int calculateHist = 0; //calculateHist的初值是0
bool showHist = false;
Point origin;
Rect selection;
//用鼠标选择目标图像区域
static void onMouse( int event, int x, int y, int, void* );
//该函数用于计算给定矩形区域的图像块的直方图
static void CalculatePatchHistogram(Rect& patchRect);
//绘制直方图
void DrawHistogramImage();
//设定直方图中bin的数目
int histSize = 256;
//设定取值范围
float range[] = {0,256};//灰度特征的范围都是0到255
const float* histRange = {range};
bool uniform = true;
bool accumu = false;//均匀分布，无累加
int _tmain(int argc, _TCHAR* argv[])
{
string filename = "E:\\素材\\flower1.jpg";
string winname = "MainWindow";
namedWindow(winname,1);
//为 Histogram Demo窗口添加鼠标响应函数，用于选择图像块
setMouseCallback( "MainWindow", onMouse, 0 );
Mat src, dst;//声明原始图像和目标图像
/// 装载图像
src = imread( filename, 1 );
if( !src.data ) //判断图像是否成功读取
{ return -1; }
//src.copyTo(image);
cvtColor(src,image,CV_BGR2GRAY);
for(;;)
{
///循环显示读入的图像
imshow(winname,image);
if(calculateHist)
{
cout<<"开始计算所选区域直方图.... ... ..."<(i-1))),Scalar(0,0,255),1,8,0);
*/
//折线表示
line( histImage, Point( bin_w*(i-1), hist_h - cvRound(roi_hist.at(i-1)) ) ,
Point( bin_w*(i), hist_h - cvRound(roi_hist.at(i)) ),
Scalar( 0, 0, 255), 1, 8, 0 );
}
/// 显示直方图
namedWindow("ROI Hist",1);
imshow("ROI Hist", histImage );
}
需要看客们注意的地方是：我并没有把TBB下的函数拷贝到自定义的命名空间钟来，因为那是并行运算的代码。
如下所示：
只要没有定义HAVE_TBB标志，自己的命名空间就可以正常运行了，否则会提示你上述类calcHist1D_Invoker<T>没有定义喔#ifdef HAVE_TBBcalcHist1D_Invoker<T> body(_ptrs, _deltas, hist, _uniranges, size[0], dims, imsize);parallel_for(BlockedRange(0, imsize.height), body);return;#endif
立足具身智能前沿赛道，致力于搭建全球化、开源化、全栈式技术交流与实践共创平台。
更多推荐
1
0- 0
已为社区贡献1条内容
所有评论(0)