根据MATLAB的histeq函数改写的运行在OpenCV下的直方图规定化C源码!

据说，图像的直方图规定化比直方图均衡化用得更多，但是很奇怪的是OpenCV居然没有图像直方图规定化的源码！所以，我就有必要在OpenCV下写一个图像直方图规定化处理的函数，以方便将来使用。

我在网上找了几个直方图均稀化的源码，并基于OpenCV来改写这些源码，效果都不如MATLAB的histeq函数，这其中改写的艰辛与繁琐就不细说了。最后，没办法，只好学习MATALB的histeq函数源码，并对其进行基于OpenCV的改写。

虽然我最终改写成功了，但是对算法还是不太理解，只能按照MATLAB的帮助文档中对histeq的讲解，大致理解其原理，其解释如下：

When you supply a desired histogram hgram, histeq chooses the grayscale transformation T to minimize



where c0 is the cumulative histogram of A, c1 is the cumulative sum of hgram for all intensities k. This minimization is subject to the constraints that T must be monotonic(单调的) and c1(T(a)) cannot overshoot(超过) c0(a) by more than half the distance between
the histogram counts at a. histeq uses the transformation b = T(a) to map the gray levels in X (or the colormap) to their new values.

在改写成功的基础上，我来试着翻译一下上面这段话哈：

对于histeq函数如果你提供了第二个参数，并且这个第二个参数是一个直方图向量，那么histeq函数将出计算一个映射T，用这个映射来将原图像的像素值映射到直方图匹配后的值。

这个映射T满足以下要求：

①保证

是最小的，其中c0是要作直方图规定化的图像的直方图累积函数，c1则是标准直方图的直方图累积函数，

②T是单调递增的

③经过T映射之后，c1在a点的值不能直那家c0在a点值的一半

原理就讲上面这么多，下面给源码：

源码中所需图片的下载链接为：

coins.png http://pan.baidu.com/s/1slilbPF
rice.png http://pan.baidu.com/s/1cDNVx8
先给出从MATLAB的histeq函数提取出的可作直方图规定化的源码：

clear all;
close all;
clc;
I1=imread('coins.png');
I2=imread('rice.png');

a=I1;

hgram=imhist(I2);

n = 256;
hgram = hgram*(numel(a)/sum(hgram));       % Set sum = numel(a)
m = length(hgram);

% [nn,cum] = computeCumulativeHistogram(a,n);
nn = imhist(a,n)';
cum = cumsum(nn);

% T = createTransformationToIntensityImage(a,hgram,m,n,nn,cum);
% Create transformation to an intensity image by minimizing the error
% between desired and actual cumulative histogram.

cumd = cumsum(hgram*numel(a)/sum(hgram));

tol = ones(m,1)*min([nn(1:n-1),0;0,nn(2:n)])/2;
err = (cumd(:)*ones(1,n)-ones(m,1)*cum(:)')+tol;
d = find(err < -numel(a)*sqrt(eps));
if ~isempty(d) %如果d为空 那么isempty的值为1,~isempty的值为0,则不会执行这个if语句
err(d) = numel(a)*ones(size(d));
end
[dum,T] = min(err); %
T = (T-1)/(m-1);%这是把灰度级由1到256变为0到255
T=T*255; %T中存储的是灰度级映射关系，比如50的值为44，则代表原图中灰度级为49的像素新的灰度级为44

再给利用MATLAB的histeq函数作直方图规定化的源码：

clear all;
close all;
clc;
I1=imread('coins.png');
I2=imread('rice.png');
histv=imhist(I2);
histv=histv';
fid=fopen('C:\Users\Administrator\Documents\MATLAB\histv3.txt','wt');
fprintf(fid,'%g ',histv);
fclose(fid);
histv=histv';%再转置回来,下面要用
J=histeq(I1,histv);
ranghistv=histv./65536;

histI1=imhist(I1);

subplot(1,2,1);
imshow(I1);
subplot(1,2,2);
imshow(J);

再给OpenCV环境下的C代码源码：

#include <opencv2/opencv.hpp>
#include <opencv2/legacy/compat.hpp>
#include <fstream>
using namespace std;
#pragma comment(linker, "/subsystem:\"windows\" /entry:\"mainCRTStartup\"")

void mycvCalcHist(IplImage *img,double out_hist[256])//自己写的直方图计算函数
{
int i=0, j=0;
double temp1=0;
int temp2=0;
const int hist_sz = 256;//0到255，一共256个灰度值
double hist[hist_sz];
memset(hist, 0, sizeof(hist));
for( i = 0; i < img->height; i++ )
{
for( j = 0; j < img->width; j++ )
{	temp1=cvGet2D(img,i,j).val[0];
temp2=int(temp1);//作类型转换
hist[temp2]++; //这里实现了hist中存储各灰度值出现的次数
}
}

memcpy(out_hist,hist, sizeof(hist)); //肯定有人要问,为啥不用数组名作为参数传递从而改变实参数组的值
//这种方法一般情况下都可以,我也测试了,然而这里就是不行,我估计与
//memset(hist, 0, sizeof(hist));这句语句有关

}

void my_hist_specification(IplImage *a,IplImage *I2)//a是源图像的指针,I2是a作直方图规定化所需的标准直方图对应图像的指针
{

int i=0;
int j=0;
int n=256;

// 计算图像的灰度直方图
double hgram[256];
mycvCalcHist(I2,hgram);
double nn[256];
mycvCalcHist(a,nn);

//相当于M中的hgram = hgram*(numel(a)/sum(hgram));
double sum_hgram=0;
for(i = 0; i < 256; i++)
sum_hgram=sum_hgram+hgram[i];
double M_a,N_a;
M_a=a->height;
N_a=a->width;
double numel_a;
numel_a=M_a*N_a;
double numel_aDivsum_hgram;
numel_aDivsum_hgram=numel_a/sum_hgram;
for(i = 0; i < 256; i++)
hgram[i] = hgram[i]*numel_aDivsum_hgram ;

//相当于m = length(hgram);
int m;
m=256;

//相当于cum = cumsum(nn);
double cum[256];
double val_1=0;
int index;
for (index = 0; index<256; index++)
{
val_1 = val_1+nn[index];
cum[index] = val_1;
}

//相当于cumd = cumsum(hgram*numel(a)/sum(hgram));
double cumd[256];
double cumd_temp1[256];
double cumd_temp2;
sum_hgram=0;
for(i = 0; i < 256; i++)
sum_hgram=sum_hgram+hgram[i];
cumd_temp2=numel_a/sum_hgram;
for(i = 0; i < 256; i++)
cumd_temp1[i]=hgram[i]*cumd_temp2;
val_1=0;
for (index = 0; index<256; index++)
{
val_1 = val_1+cumd_temp1[index];
cumd[index] = val_1;
}

//相当于tol = ones(m,1)*min([nn(1:n-1),0;0,nn(2:n)])/2;
double tol_temp1[256];
for(index = 0; index<256; index++)
tol_temp1[index]=nn[index]/2;
tol_temp1[0]=0;
tol_temp1[255]=0;
static double tol[256][256];
for(index = 0; index<256; index++)
memcpy(tol[index],tol_temp1, sizeof(tol_temp1)); //这里值是正确的，相当于对tol_temp1作了行复制；

//相当于err = (cumd(:)*ones(1,n)-ones(m,1)*cum(:)')+tol;
static double err_1[256][256];
for(i=0;i<256;i++)
for(j=0;j<256;j++)
err_1[i][j]=cumd[i];   //这里值是正确的，对cumd作了列复制；二次检查也无误

static double err_2[256][256];
for(index = 0; index<256; index++)
memcpy(err_2[index],cum, sizeof(cum)); //这里值是正确的，对cum作了行复制
static double err[256][256];
for(i=0;i<256;i++)
for(j=0;j<256;j++)
err[i][j]=err_1[i][j]-err_2[i][j]+tol[i][j];

//相当于d = find(err < -numel(a)*sqrt(eps));
double matlab_eps= 1.4901e-008;
double find_err_temp1=-numel_a*matlab_eps;
int k=0;
int breakvari=0;
double watch_err;
for(i=0;i<256;i++)  //这里是按列遍历，不是按行遍历
{
for(j=0;j<256;j++)
if(err[j][i]<find_err_temp1)
{	watch_err=err[j][i];
k++;
}
}
double *d = (double *)malloc(k*sizeof(double));
k=0;
for(i=0;i<256;i++)  //这里是按列遍历，不是按行遍历
{
for(j=0;j<256;j++)
if(err[j][i]<find_err_temp1)
{	watch_err=err[j][i];
d[k]=i*256+j;
k++;
}
}

/*
相当于
if ~isempty(d) %如果d为空 那么isempty的值为1,~isempty的值为0,则不会执行这个if语句
err(d) = numel(a)*ones(size(d));
end
*/
int d_size;
d_size=k;
char isempty_flag=1;
for(i=0;i<d_size;i++)
{
if(d[i]!=0)
{
isempty_flag=0;
break;
}
}
double err_temp1=0;
int d_m,d_n;
if(!isempty_flag)
{
for(i=0;i<d_size;i++)
{
d_m=int(d[i])%256;
d_n=int(d[i])/256;
err[d_m][d_n]=numel_a;

}

}
//相当于[dum,T] = min(err);
int T[256];
double err_min=0;
for(j=0;j<256;j++)
{	err_min=err[0][j];
T[j]=0;
for(i=0;i<256;i++)
{   if(err[i][j]<err_min)
{
T[j]=i;
err_min=err[i][j];//T中就是经过直方图匹配之后的像素值映射关系了！
}
}
}

//下面的程序是把原图像中的每一个像素值用T中的映射关系来映射
T[0] = 0;//这是我通过前面的程序计算出的灰度值映射关系,不管你怎么映射,0肯定是0
int s1_temp1;
CvScalar s1;
CvScalar s2;
for( i = 0; i < M_a; i++ )
{
for( j = 0; j < N_a; j++ )
{
s1 = cvGet2D(a,i,j);
s1_temp1=int(s1.val[0]);
s2.val[0]=T[s1_temp1];
cvSet2D(a,i,j,s2);
}

}

}

int main()
{
// 从文件中加载原图
IplImage *pSrcImage_coins = cvLoadImage("coins.png", CV_LOAD_IMAGE_UNCHANGED);

IplImage *a = cvLoadImage("coins.png", CV_LOAD_IMAGE_UNCHANGED);
IplImage *I2 = cvLoadImage("rice.png", CV_LOAD_IMAGE_UNCHANGED);

my_hist_specification(a,I2);

const char *pstrWindowsATitle = "原图";
const char *pstrWindowsBTitle = "直方图规定化的源图";
//创建窗口
cvNamedWindow(pstrWindowsATitle, CV_WINDOW_AUTOSIZE);
cvNamedWindow(pstrWindowsBTitle, CV_WINDOW_AUTOSIZE);
//在指定窗口中显示图像
cvShowImage(pstrWindowsATitle, pSrcImage_coins);
cvShowImage(pstrWindowsBTitle, a);
//等待按键事件
cvWaitKey();
cvDestroyWindow(pstrWindowsATitle);
cvDestroyWindow(pstrWindowsBTitle);
cvReleaseImage(&a);
cvReleaseImage(&I2);

return 0;
}

最后的运行结果如下图所示：





我用另外一幅图测试，也没有问题，如下：





可见，结果是一致的，程序测试成功！另外，在我备份的源码Histogram_Specification_05.cpp，下载链接为：http://pan.baidu.com/s/1i5g0Kzr 完整的体现了我整个改写的过程，我设置了很多中间变量和数组，可以观察程序是否正确！

欢迎大家加入图像识别技术交流群：271891601，另外，特别欢迎成都从事图像识别工作的朋友交流，我的QQ号248787278