根據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到255T=T*255; %T中儲存的是灰階級映射關係，比如50的值為44，則代表原圖中灰階級為49的像素新的灰階級為44

可出最後經長條圖規定化處理後的映像的MATLAB源碼如下：

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  完整的體現了我整個改寫的過程，我設定了很多中間變數和數組，可以觀察程式是否正確！

歡迎大家加入Image Recognition技術交流群：271891601，另外，特別歡迎成都從事Image Recognition工作的朋友交流，我的QQ號248787278

根據MATLAB的histeq函數改寫的運行在OpenCV下的長條圖規定化C源碼!