空間域圖像處理直方圖算法-直方圖匹配，局部直方圖與直方圖統計算法

時間 2019-12-14

標籤空間圖像處理直方圖算法匹配局部統計简体版

原文原文鏈接

直方圖均衡化算法能夠本身獲得一個轉換函數，將輸出圖像具備近似的均勻分佈。直方圖均衡化是結果可預測且容易實現。但對於一些特殊的案例，直方圖均衡化試圖獲得均勻直方圖的方法並不能達到效果，這類狀況下，每每須要指定輸出圖像直方圖的具體分佈，可以輸出具備指定分佈直方圖的算法就是直方圖匹配；ios

算法原理：c++

　　仍然用\(r\)和\(z\)表示連續灰度變量，\(p_r(r)\)和\(p_z(z)\)表示連續機率密度函數，\(s\)是一個隨機變量，算法

\(s = T(r) = \int_{0}^{r} p_r(w)dw\) (1)數組

　　隨機變量\(z\)具備性質：函數

　　　　\(G(z) = \int_{0}^{z} p_z(t)dt = s\)　　(2)測試

　　所以就有了，\(G(z) = T(r)\) 而且ui

\(z = G^{-1}(s) = G^{-1} \left[T(r)\right]\)　　(3)this

若是\(G^{-1}\)存在，且知足條件：單值且單調遞增，值域爲[0,1],方程1~3就定義了一個輸出指定直方圖匹配的過程。從輸入圖像獲得\(s\)，從指定的直方圖分佈獲得\(z\)，再有方程3就有了從輸入\(r\)到輸出\(s\)的映射。spa

但問題是，在實際實踐中很難得到\(G^{-1}\)和\(T(r)\)的解析式。比較幸運的是，在離散狀況下，問題會變得簡單，最終獲得的是理想指定分佈的近似。3d

離散環境下：

\(s_k = T(r_k) = \sum_{j=0}^{k} P_r(r_j) = \sum_{j=0}^{k} \frac{n_j}{n} \qquad k =0,1,2,...L-1\)　　(4)

\(v_k = G(z_k) = \sum_{i=0}^{k} p_z(z_j) = s_k \qquad k = 0,1,2,3,...L-1.\)　　(5)

\(z_k = G^{-1}\left[T(r_k)\right] \qquad k = 0,1,2,3,...L-1\)　　(6)

\(z_k = G^{-1}(s_k)\)

方程(4)將像素點從\(r\)的灰度空間映射到了\(s\)的灰度空間，而方程(5)定義了一個直方圖隨機變量的轉換函數，從\(z\)的灰度空間到s的空間。方程(6)是(5)的逆變換，獲得了具備指定直方圖分佈的隨機變量。

其實，通俗的講，有三個灰度空間\({r_j}\),\({s_j}\),\({z_j}\), \(j = 0,1,2,...L-1\)的一維數組。映射就是完成這樣的過程：指定下標查找；下標構成了一個查找表，對於灰度級j，\(r_j\)映射到\(s_j\),\(s_j\)映射到\(z_j\)，\(r_j\)和\(z_j\)都是[0,L-]的灰度空間，\(s\)是[0,1]的空間。

圖（1）闡明瞭這樣的關係：

圖（1）

左上角的圖表示了一個假設的轉換函數，對應於每個\(z_q\)，輸出轉換函數值\(v_q\)，事實上咱們須要的是這樣的反過程，由於\(v = s \)，因此要找的就是爲每個\(s_k\)，找到對應的\(z_k\)。

正由於\(v = s\)，而且處理的是離散整型數據，因此尋找的\(z_k\)應該要知足條件\(G(z_k) - s_k= 0\)，而因爲\(z\)是一整數位步長的，因此要作的只是從小到大遍歷\(z_i\)，若是\(z_i\)是知足條件

\(G(z_i)-s_k)\ge 0\)　　(7)

的最小整數，那麼\(z_k = z_i\);

總結一下步驟：

（1）：根據輸入圖像獲得圖像的灰度直方圖（離散機率密度）；

（2）：根據方程（4）獲得\(s\)

（3）：由具備指定機率分佈的隨機變量\(z\)根據方程（5）獲得G；

（4）：以整數步長迭代\(z\)，找到知足方程（7）的最小整數\(z_i\),即爲對應\(s_k\)的\(z_k\)值；

實現：

MATLAB提供的histeq函數便可以作直方圖均衡化，也能夠作直方圖匹配，只需將指定的機率密度函數做爲參數傳入便可。

輸入圖像：

 1 >>f=imread('D:\MATLAB\image\DIP3E_Original_Images_CH03\Fig0323(a)(mars_moon_phobos).tif');
 2 >> hig = imhist(f);
 3 >> hig1 = hig(1:256);
 4 >> horz = 1:256;
 5 >> figure,stem(horz,hig1,'fill')
 6 >> axis([0 255 0 60000])
 7 >> set(gca,'xtick',[0:16:256])
 8 >> set(gca,'ytick',[0:6000:60000])
 9 >> figure,plot(horz,hig1,'color','b','linestyle','-','marker','.')
10 >> axis([0 255 0 60000])
11 >> set(gca,'xtick',[0:16:256])
12 >> set(gca,'ytick',[0:6000:60000])

圖（2）

輸入圖像及其直方圖，大部分像素點集中在較暗區域，直方圖的左側，若是對此作直方圖均衡化，左邊大量像素點將映射到較高的灰度區，並沒有法分散像素點，沒法提升對比度。

採用直方圖匹配，首先要針對圖像擬合一個能夠提升對比度的輸出圖像的灰度機率分佈函數，原圖有兩個峯值，將峯值出比較陡的灰度分佈變爲高斯分佈，也許能夠離散密集的灰度值，平緩峯值。能夠根據方程（4）求得\(s\)，s的函數圖像爲：

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　圖（3）

擬合一個具備兩個高斯峯值的輸出函數，MATLAB代碼：

 1 function p= twomodegauss( u1,sig1,u2,sig2,A1,A2,k)
 2 %bimodal gaussian function
 3 %   Detailed explanation goes here
 4 c1 = A1 *(1/((2*pi)^0.5)*sig1);
 5 c2 = A2 *(1/((2*pi)^0.5)*sig2);
 6 k1 = 2 * (sig1^2);
 7 k2 = 2 * (sig2^2);
 8 z = linspace(0,1,256);
 9 p = k + c1*exp(-((z-u1).^2)./k1)+c2 * exp(-((z-u2).^2)./k2);
10 p = p./sum(p(:));
11 end

 1 function p = manualhist
 2 %UNTITLED2 Summary of this function goes here
 3 %   Detailed explanation goes here
 4 flag = true;
 5 quitnow = 'x';
 6 p = twomodegauss(0.15,0.05,0.75,0.05,1,0.07,0.002);
 7 while flag
 8     s = input('enter m1,sig1,m2,sig2,A1,A2,k or x to quit:','s');
 9     if s== quitnow
10         break
11     end
12     v = str2num(s);
13     if numel(v) ~= 7
14         disp('incorret number of imputs')
15         continue
16     end
17     p  = twomodegauss(v(1),v(2),v(3),v(4),v(5),v(6),v(7));
18 end
19     figure,plot(p)
20     xlim([0 255])
21 end

輸出高斯函數及其CDF（\(G\)）：

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　圖（4）

輸出圖像要具備圖（4）所示的灰度機率密度分佈，須要完成直方圖匹配過程，也就是方程（5）與步驟（4），MATLAB實現爲histeq函數：

 1 >> g = histeq(f,p);
 2 >> imshow(g)
 3 Warning: Image is too big to fit on screen; displaying at 67% 
 4 > In imuitools\private\initSize at 73
 5   In imshow at 262
 6 >> hg = imhist(g);
 7 >> hg1 = hg(1:256);
 8 >> horz = 1:256;
 9 >> plot(horz,hg1,'color','b','linestyle','-','marker','.');
10 >> axis([0 255 0 60000])
11 >> set(gca,'xtick',[0:10:256])
12 >> set(gca,'ytick',[0:3000:60000])

輸出效果圖像：

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　圖（5）

實現直方圖匹配的C++算法：

 1 #define A1 1
 2 #define A2 0.07
 3 #define U1 0.15
 4 #define SIG1 0.05
 5 #define U2 0.75
 6 #define SIG2 0.05
 7 #define pi 3.1415626535627
 8 using namespace cimg_library;
 9 
10 inline float twoModeGauss(float r){
11     //double w = K1+(A1/sqrt(2*pi)*SIG1)*exp(-(pow((r-U1),2.0)/2*pow(SIG1,2.0)))+(A2/sqrt(2*pi)*SIG2)*exp(-(pow((r-U2),2.0)/2*pow(SIG2,2.0)));
12     float x = (A1/sqrt(2*pi)*SIG1);
13     float y = (A2/sqrt(2*pi)*SIG2);
14     float a1 = exp(-(pow((r-U1),2.0)/(2*pow(SIG1,2.0))));
15     float a2 = exp(-(pow((r-U2),2.0)/(2*pow(SIG2,2.0))));
16     float w =K1+x*a1+a2*y;
17     return w;
18 }
19 void main(int argc,char ** argv)
20 {
21     CImg<unsigned char> img("D:\\MATLAB\\image\\DIP3E_Original_Images_CH03\\Fig0323(a)(mars_moon_phobos).tif");
22     int width = img.width();
23     int height = img.height();
24     int consus[256];
25     for(size_t m = 0;m<256;m++){
26         consus[m] = 0;
27     }
28     unsigned sub = 0;
29     double sum = width * height;
30     for(int i = 0;i<width;i++){
31         for(int j = 0;j<height;j++){
32             sub = img(i,j);
33             consus[sub]++;
34         }
35     }
36     double s[256];
37     s[0] = consus[0]/sum;
38     for(size_t n = 1;n<256;n++){
39         s[n] = s[n-1]+(consus[n]/sum);
40     }
41     double count[256];
42     sum = 0;
43     for(size_t a = 0;a<256;a++){
44         double r = (double)a/256.0;
45         count[a] = twoModeGauss(r);
46         sum+=count[a];
47     }
48     for(size_t c = 0;c<256;c++){
49         count[c]=count[c]/sum;
50     }
51     double v = 0;
52     unsigned char z[256];
53     size_t p = 0;
54     for(size_t b = 0;b<256;b++){
55         while(s[b]>v){
56             v+=count[p];
57             p++;
58         }
59         z[b] = p;
60     }
61     CImg<unsigned char>dist(img);
62     for(int i = 0;i<width;i++){
63         for(int j = 0;j<height;j++){
64             dist(i,j) = z[img(i,j)];
65         }
66     }
67     CImgDisplay disp1(img,"origin");
68     CImgDisplay disp2(dist,"distination");
69     disp1.wait();
70     disp2.wait();
71 }

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　圖（6）

c++匹配後的效果在某些位置有些失真，也許代碼有待調整。

改正代碼：

 1 #include "CImg.h"
 2 #include <iostream>
 3 #define K1 0.002
 4 #define A1 1
 5 #define A2 0.07
 6 #define U1 0.15
 7 #define SIG1 0.05
 8 #define U2 0.75
 9 #define SIG2 0.05
10 #define pi 3.1415626535627
11 using namespace cimg_library;
12 /**************直方圖匹配*********************************
13 *匹配分佈函數：雙峯值高斯函數；                            *
14 *圖像讀取顯示庫：CImg；                                   *
15 *時間複雜度：O(N^2)，操做二維圖像,平方級出如今逐像素統計；   *
16 *author：Bruce Mu                                       *
17 *基礎測試版                                              *
18 ********************************************************/
19 //雙峯值高斯函數
20 inline double twoModeGauss(double r){
21     double w = K1+(A1/sqrt(2*pi)*SIG1)*exp(-(pow((r-U1),2.0)/(2*pow(SIG1,2.0))))+(A2/sqrt(2*pi)*SIG2)*exp(-(pow((r-U2),2.0)/(2*pow(SIG2,2.0))));
22     return w;
23 }
24 void main(int argc,char ** argv)
25 {
26     CImg<unsigned char> img("D:\\MATLAB\\image\\DIP3E_Original_Images_CH03\\Fig0323(a)(mars_moon_phobos).tif");
27     int width = img.width();
28     int height = img.height();
29     int consus[257];
30     for(size_t m = 0;m<257;m++){
31         consus[m] = 0;
32     }
33     unsigned sub = 0;
34     double sum = width * height;
35     //統計直方圖
36     for(int i = 0;i<width;i++){
37         for(int j = 0;j<height;j++){
38             sub = img(i,j);
39             consus[sub+1]++;
40         }
41     }
42 
43     //頻數轉頻率 並求CDF
44     double s[257];
45     s[0] = consus[0]/sum;
46     for(size_t n = 1;n<257;n++){
47         s[n] = s[n-1]+(consus[n]/sum);
48     }
49     //求取指定直方圖的CDF
50     double count[257];
51     count[0] = 0;
52     sum = 0;
53     for(size_t a = 0;a<257;a++){
54         double r = (double)a/256.0;
55         count[a+1] = twoModeGauss(r);
56         sum+=count[a+1];
57     }
58     for(size_t c = 1;c<257;c++){
59         count[c]=count[c]/sum;
60     }
61     //根據兩個CDF求取灰度映射
62     double v = count[0];
63     unsigned char z[256];
64     size_t p = 1;
65     for(size_t b = 0;b<256;b++){
66         while(s[b]>v){
67             v+=count[p];
68             p++;
69         }
70         z[b] = p;
71     }
72     //按像素點的映射值賦值目標圖
73     CImg<unsigned char>dist(img);
74     for(int i = 0;i<width;i++){
75         for(int j = 0;j<height;j++){
76             dist(i,j) = z[img(i,j)];
77         }
78     }
79     //顯示
80     CImgDisplay disp1(img,"origin");
81     CImgDisplay disp2(dist,"distination");
82     disp1.wait();
83     disp2.wait();
84 }

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　圖（7）

局部直方圖：

直方圖均衡化和直方圖匹配都是全局的加強方法，若是對於一些例子，圖像的某些部分無需改動，而也存在部分須要加強效果，就須要用到局部直方圖。局部直方圖的思想在於並不是對整個圖像採用直方圖算法，而是改變目標爲針對像素點的鄰域。鄰域內直方圖均衡化或直方圖匹配，目的是加強小區域中的細節。這些小區域每每有這樣的特色，區域內像素數對全局算法和變換隻起到微不足道的做用，因此全局算法很難做用到這樣的小區域，而反之小區域內的加強變換並不會影響全局。

實現：逐像素取鄰域，應用直方圖算法修改中心點像素值；

使用CImg的C局部直方圖算法實現：

 1 #include <iostream>
 2 #include "CImg.h"
 3 #include <math.h>
 4 using namespace cimg_library;
 5 //根據統計好的直方圖求鄰域內的CDF
 6 void getCDF(float *s,const float *const c){
 7     for(size_t i = 1;i<256;i++){
 8         s[i] = s[i-1]+c[i];
 9     }
10 }
11 //局部變換和鄰域中心移動後修改統計信息；
12 inline void refresh(float *c,unsigned char sub,BOOL flag){
13     if(flag){
14         c[sub]+=1.0/49.0;
15     }else{
16         c[sub]-=1.0/49.0;
17     }
18 }
19 //應用局部直方圖算法
20 void localhist(CImg<unsigned char> &img){
21     int width = img.width();
22     int height = img.height();
23     float count[256];
24     for(size_t t= 0;t<256;t++){
25         count[t]= 0;
26     }
27     //統計第一鄰域
28     unsigned char sub = 0;
29     for(int m = 0;m<7;m++){
30         for(int n = 0;n<7;n++){
31             sub = img(m,n);  
32             count[sub]+=(1.0/49.0);
33         }
34     }
35     float s[257];
36     s[0] = count[0];
37     getCDF(s,count);
38     sub = img(3,3);
39     unsigned char newvalue = (unsigned char)((s[sub])*255+0.5);
40     img(3,3) = newvalue;
41     refresh(count,sub,false);
42     refresh(count,newvalue,true);
43     //鄰域位移，逐鄰域修改直方圖信息：
44     for(int i = 3;i<height - 4;i++){
45         for(int j = 3;j<width-4;j++){
46             BOOL dirty = true;
47             if(i == 3 && j == 3){
48                 continue;
49             }
50             for(int p = 0;p<7;p++){
51                 unsigned char old =img(j-4,i-3+p);
52                 unsigned char add = img(j+3,i-3+p);
53                 if(old != add){
54                     refresh(count,old,false);
55                     refresh(count,add,true);
56                     dirty = false;
57                 }
58             }
59             if(!dirty){
60                 getCDF(s,count);
61                 sub = img(j,i);
62                 refresh(count,sub,false);
63                 img(j,i) = (unsigned char)(s[sub]*255+0.5);
64                 refresh(count,img(j,i),true);
65             }
66         }
67     }
68 }
69 
70 void main(int argc,char **argv){
71     CImg<unsigned char>img("E:\\convert.jpeg");
72     CImg<unsigned char>origin(img);
73     CImgDisplay disp1(origin,"origin");
74     localhist(img);
75     CImgDisplay disp2(img,"change");
76     while (!disp1.is_closed() && !disp2.is_closed()) {
77         disp1.wait();
78         disp2.wait();
79     }
80 }

算法效果：

鄰域內直方圖統計算法：

  1 #include <iostream>
  2 #include "CImg.h"
  3 using namespace std;
  4 using namespace cimg_library;
  5 
  6 float getMean(const float *c){
  7     float m = 0;
  8     for(size_t x = 0;x<256;x++){
  9         m+=x*c[x];
 10     }
 11     return m;
 12 }
 13 
 14 float getVariance(const float *c,float mean){
 15     float delta = 0;
 16     for(size_t y = 0;y<256;y++){
 17         if(c[y] != 0){
 18             delta += (pow((y-mean),2)*c[y]);
 19         }
 20     }
 21     delta = sqrt(delta);
 22     return delta;
 23 }
 24 
 25 void getCount(float *c,int originx,int originy,int wid,int hei,CImg<unsigned char> &img){
 26     for(size_t t = 0;t<256;t++){
 27         c[t] = 0;
 28     }
 29     for(int i = originx;i<wid;i++){
 30         for(int j = originy;j<hei;j++){
 31             unsigned char sub = img(i,j);
 32             c[sub] += 1.0/256.0;
 33         }
 34     }
 35 }
 36 
 37 void init(float * count){
 38     for(size_t t = 0;t<256;t++){
 39         count[t] = 0;
 40     }
 41 }
 42 int histStatistic(CImg<unsigned char> &img,int dim,float k0,float k1,float k2,float E){
 43     if(dim%2 == 0){
 44         cout<<"please input a odd num!!!"<<endl;
 45         return -1;
 46     }
 47     int width = img.width();
 48     int height = img.height();
 49     int sum = width * height;
 50     float count[256];
 51     for(size_t t = 0;t<256;t++){
 52         count[t]=0;
 53     }
 54     unsigned char sub = 0;
 55     for(int i = 0;i<width;i++){
 56         for(int j = 0;j<height;j++){
 57             sub = img(i,j);
 58             count[sub]+=1.0/sum;
 59         }
 60     }
 61     float mean = getMean(count);
 62     float delta = getVariance(count,mean);
 63     float mxy = 0;
 64     float deltaxy = 0;
 65     float local = dim*dim;
 66     for(int i = dim/2;i<height - dim/2 -1;i++){
 67         for(int j = dim/2;j<width - dim/2 -1;j++){
 68             //從新初始化統計數組
 69             init(count);
 70             //統計局部直方圖
 71             //getCount(count,j-dim/2,i-dim/2,j+dim/2+1,i+dim/2+1,img);
 72             for(int p = j-dim/2;p<j+dim/2+1;p++){
 73                 for(int q = i - dim/2;q < i+dim/2+1;q++){
 74                     sub = img(p,q);
 75                     count[sub] += 1.0/local;
 76                 }
 77             }
 78             mxy = getMean(count);
 79             deltaxy = getVariance(count,mxy);
 80         //    cout<<mxy<<"  :  "<<deltaxy<<endl;
 81             //加強判斷
 82             if(mxy <= mean*k0 && deltaxy <=k2*delta && deltaxy >= k1*delta){
 83                 img(j,i) = img(j,i)*E;
 84             //    cout<<"enforcement_________________"<<endl;
 85             }
 86         }
 87     }
 88     return 0;
 89 }
 90 
 91 int main(int argc,char ** argv){
 92     CImg<unsigned char> img("E:\\statistic.jpeg");
 93     CImg<unsigned char> dist(img);
 94     histStatistic(dist,3,0.4,0.008,0.2,5.0);
 95     CImgDisplay disp1(img,"origin");
 96     CImgDisplay disp2(dist,"dist");
 97     while(!disp1.is_closed() && !disp2.is_closed()){
 98         disp2.wait();
 99     }
100 }