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

àì夳堔傛蜴生んèń 2022-08-22 01:23 31阅读 0赞

2016-6-8日用C++也实现了直方图规定化，只是结果上和这篇博文的源码有差异，详情[http://blog.csdn.net/wenhao\_ir/article/details/51605451][http_blog.csdn.net_wenhao_ir_article_details_51605451]

据说，图像的直方图规定化比直方图均衡化用得更多，但是很奇怪的是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  
  
![20160510150519935][]  
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满足以下要求：

①保证![20160510150519935][]是最小的，其中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;
    }

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

![20160510153345566][]

![20160510153358409][]

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

![20160510155618457][]![20160510155625669][]

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

图像处理开发资料、图像处理开发需求、图像处理接私活挣零花钱，可以搜索公众号"qxsf321"，并关注！  
图像处理开发资料、图像处理开发需求、图像处理接私活挣零花钱，可以搜索公众号"qxsf321"，并关注！  
图像处理开发资料、图像处理开发需求、图像处理接私活挣零花钱，可以搜索公众号"qxsf321"，并关注！

[http_blog.csdn.net_wenhao_ir_article_details_51605451]: http://blog.csdn.net/wenhao_ir/article/details/51605451
[20160510150519935]: /images/20220724/5a09f49e5ddd4558af90bcdeee4ba4b1.png
[20160510153345566]: /images/20220724/ffc6eb93f0824854b6333b8de14d1a22.png
[20160510153358409]: /images/20220724/59786d92ee8e44039011d52ca4fe38e9.png
[20160510155618457]: https://img-blog.csdn.net/20160510155618457
[20160510155625669]: /images/20220724/daaa4d0f33d14349a4149e4f76bf23bc.png