空间域图像处理直方图算法-直方图匹配，局部直方图与直方图统计算法

直方图均衡化算法可以自己得到一个转换函数，将输出图像具有近似的均匀分布。直方图均衡化是结果可预测且容易实现。但对于一些特殊的案例，直方图均衡化试图得到均匀直方图的方法并不能达到效果，这类情况下，往往需要指定输出图像直方图的具体分布，能够输出具有指定分布直方图的算法就是直方图匹配；

算法原理：

　　仍然用(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)) 并且

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

如果(G^{-1})存在，且满足条件：单值且单调递增，值域为[0,1],方程1~3就定义了一个输出指定直方图匹配的过程。从输入图像得到(s)，从指定的直方图分布得到(z)，再有方程3就有了从输入(r)到输出(s)的映射。

但问题是，在实际实践中很难获得(G^{-1})和(T(r))的解析式。比较幸运的是，在离散情况下，问题会变得简单，最终得到的是理想指定分布的近似。

离散环境下：

(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) ight] 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函数即可以做直方图均衡化，也可以做直方图匹配，只需将指定的概率密度函数作为参数传入即可。

输入图像：

>>f=imread('D:MATLABimageDIP3E_Original_Images_CH03Fig0323(a)(mars_moon_phobos).tif');
>> hig = imhist(f);
>> hig1 = hig(1:256);
>> horz = 1:256;
>> figure,stem(horz,hig1,'fill')
>> axis([0 255 0 60000])
>> set(gca,'xtick',[0:16:256])
>> set(gca,'ytick',[0:6000:60000])
>> figure,plot(horz,hig1,'color','b','linestyle','-','marker','.')
>> axis([0 255 0 60000])
>> set(gca,'xtick',[0:16:256])
>> set(gca,'ytick',[0:6000:60000])

图（2）

输入图像及其直方图，大部分像素点集中在较暗区域，直方图的左侧，如果对此做直方图均衡化，左边大量像素点将映射到较高的灰度区，并无法分散像素点，无法提高对比度。

采用直方图匹配，首先要针对图像拟合一个可以提高对比度的输出图像的灰度概率分布函数，原图有两个峰值，将峰值出比较陡的灰度分布变为高斯分布，也许可以离散密集的灰度值，平缓峰值。可以根据方程（4）求得(s)，s的函数图像为：

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　图（3）

 拟合一个具有两个高斯峰值的输出函数，MATLAB代码：

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

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

 输出高斯函数及其CDF（(G)）：

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　图（4）

输出图像要具有图（4）所示的灰度概率密度分布，需要完成直方图匹配过程，也就是方程（5）与步骤（4），MATLAB实现为histeq函数：

>> g = histeq(f,p);
>> imshow(g)
Warning: Image is too big to fit on screen; displaying at 67% 
> In imuitoolsprivateinitSize at 73
  In imshow at 262
>> hg = imhist(g);
>> hg1 = hg(1:256);
>> horz = 1:256;
>> plot(horz,hg1,'color','b','linestyle','-','marker','.');
>> axis([0 255 0 60000])
>> set(gca,'xtick',[0:10:256])
>> set(gca,'ytick',[0:3000:60000])

 输出效果图像：

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　图（5）

实现直方图匹配的C++算法：

#define A1 1
#define A2 0.07
#define U1 0.15
#define SIG1 0.05
#define U2 0.75
#define SIG2 0.05
#define pi 3.1415626535627
using namespace cimg_library;

inline float twoModeGauss(float r){
    //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)));
    float x = (A1/sqrt(2*pi)*SIG1);
    float y = (A2/sqrt(2*pi)*SIG2);
    float a1 = exp(-(pow((r-U1),2.0)/(2*pow(SIG1,2.0))));
    float a2 = exp(-(pow((r-U2),2.0)/(2*pow(SIG2,2.0))));
    float w =K1+x*a1+a2*y;
    return w;
}
void main(int argc,char ** argv)
{
    CImg<unsigned char> img("D:\MATLAB\image\DIP3E_Original_Images_CH03\Fig0323(a)(mars_moon_phobos).tif");
    int width = img.width();
    int height = img.height();
    int consus[256];
    for(size_t m = 0;m<256;m++){
        consus[m] = 0;
    }
    unsigned sub = 0;
    double sum = width * height;
    for(int i = 0;i<width;i++){
        for(int j = 0;j<height;j++){
            sub = img(i,j);
            consus[sub]++;
        }
    }
    double s[256];
    s[0] = consus[0]/sum;
    for(size_t n = 1;n<256;n++){
        s[n] = s[n-1]+(consus[n]/sum);
    }
    double count[256];
    sum = 0;
    for(size_t a = 0;a<256;a++){
        double r = (double)a/256.0;
        count[a] = twoModeGauss(r);
        sum+=count[a];
    }
    for(size_t c = 0;c<256;c++){
        count[c]=count[c]/sum;
    }
    double v = 0;
    unsigned char z[256];
    size_t p = 0;
    for(size_t b = 0;b<256;b++){
        while(s[b]>v){
            v+=count[p];
            p++;
        }
        z[b] = p;
    }
    CImg<unsigned char>dist(img);
    for(int i = 0;i<width;i++){
        for(int j = 0;j<height;j++){
            dist(i,j) = z[img(i,j)];
        }
    }
    CImgDisplay disp1(img,"origin");
    CImgDisplay disp2(dist,"distination");
    disp1.wait();
    disp2.wait();
}

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　图（6）

c++匹配后的效果在某些位置有些失真，也许代码有待调整。

改正代码：

#include "CImg.h"
#include <iostream>
#define K1 0.002
#define A1 1
#define A2 0.07
#define U1 0.15
#define SIG1 0.05
#define U2 0.75
#define SIG2 0.05
#define pi 3.1415626535627
using namespace cimg_library;
/**************直方图匹配*********************************
*匹配分布函数：双峰值高斯函数；                            *
*图像读取显示库：CImg；                                   *
*时间复杂度：O(N^2)，操作二维图像,平方级出现在逐像素统计；   *
*author：Bruce Mu                                       *
*基础测试版                                              *
********************************************************/
//双峰值高斯函数
inline double twoModeGauss(double r){
    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))));
    return w;
}
void main(int argc,char ** argv)
{
    CImg<unsigned char> img("D:\MATLAB\image\DIP3E_Original_Images_CH03\Fig0323(a)(mars_moon_phobos).tif");
    int width = img.width();
    int height = img.height();
    int consus[257];
    for(size_t m = 0;m<257;m++){
        consus[m] = 0;
    }
    unsigned sub = 0;
    double sum = width * height;
    //统计直方图
    for(int i = 0;i<width;i++){
        for(int j = 0;j<height;j++){
            sub = img(i,j);
            consus[sub+1]++;
        }
    }

    //频数转频率 并求CDF
    double s[257];
    s[0] = consus[0]/sum;
    for(size_t n = 1;n<257;n++){
        s[n] = s[n-1]+(consus[n]/sum);
    }
    //求取指定直方图的CDF
    double count[257];
    count[0] = 0;
    sum = 0;
    for(size_t a = 0;a<257;a++){
        double r = (double)a/256.0;
        count[a+1] = twoModeGauss(r);
        sum+=count[a+1];
    }
    for(size_t c = 1;c<257;c++){
        count[c]=count[c]/sum;
    }
    //根据两个CDF求取灰度映射
    double v = count[0];
    unsigned char z[256];
    size_t p = 1;
    for(size_t b = 0;b<256;b++){
        while(s[b]>v){
            v+=count[p];
            p++;
        }
        z[b] = p;
    }
    //按像素点的映射值赋值目标图
    CImg<unsigned char>dist(img);
    for(int i = 0;i<width;i++){
        for(int j = 0;j<height;j++){
            dist(i,j) = z[img(i,j)];
        }
    }
    //显示
    CImgDisplay disp1(img,"origin");
    CImgDisplay disp2(dist,"distination");
    disp1.wait();
    disp2.wait();
}

 

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　图（7）

局部直方图：

直方图均衡化和直方图匹配都是全局的增强方法，如果对于一些例子，图像的某些部分无需改动，而也存在部分需要增强效果，就需要用到局部直方图。局部直方图的思想在于并非对整个图像采用直方图算法，而是改变目标为针对像素点的邻域。邻域内直方图均衡化或直方图匹配，目的是增强小区域中的细节。这些小区域往往有这样的特点，区域内像素数对全局算法和变换只起到微不足道的作用，所以全局算法很难作用到这样的小区域，而反之小区域内的增强变换并不会影响全局。

实现：逐像素取邻域，应用直方图算法修改中心点像素值；

使用CImg的C局部直方图算法实现：

#include <iostream>
#include "CImg.h"
#include <math.h>
using namespace cimg_library;
//根据统计好的直方图求邻域内的CDF
void getCDF(float *s,const float *const c){
    for(size_t i = 1;i<256;i++){
        s[i] = s[i-1]+c[i];
    }
}
//局部变换和邻域中心移动后修改统计信息；
inline void refresh(float *c,unsigned char sub,BOOL flag){
    if(flag){
        c[sub]+=1.0/49.0;
    }else{
        c[sub]-=1.0/49.0;
    }
}
//应用局部直方图算法
void localhist(CImg<unsigned char> &img){
    int width = img.width();
    int height = img.height();
    float count[256];
    for(size_t t= 0;t<256;t++){
        count[t]= 0;
    }
    //统计第一邻域
    unsigned char sub = 0;
    for(int m = 0;m<7;m++){
        for(int n = 0;n<7;n++){
            sub = img(m,n);  
            count[sub]+=(1.0/49.0);
        }
    }
    float s[257];
    s[0] = count[0];
    getCDF(s,count);
    sub = img(3,3);
    unsigned char newvalue = (unsigned char)((s[sub])*255+0.5);
    img(3,3) = newvalue;
    refresh(count,sub,false);
    refresh(count,newvalue,true);
    //邻域位移，逐邻域修改直方图信息：
    for(int i = 3;i<height - 4;i++){
        for(int j = 3;j<width-4;j++){
            BOOL dirty = true;
            if(i == 3 && j == 3){
                continue;
            }
            for(int p = 0;p<7;p++){
                unsigned char old =img(j-4,i-3+p);
                unsigned char add = img(j+3,i-3+p);
                if(old != add){
                    refresh(count,old,false);
                    refresh(count,add,true);
                    dirty = false;
                }
            }
            if(!dirty){
                getCDF(s,count);
                sub = img(j,i);
                refresh(count,sub,false);
                img(j,i) = (unsigned char)(s[sub]*255+0.5);
                refresh(count,img(j,i),true);
            }
        }
    }
}

void main(int argc,char **argv){
    CImg<unsigned char>img("E:\convert.jpeg");
    CImg<unsigned char>origin(img);
    CImgDisplay disp1(origin,"origin");
    localhist(img);
    CImgDisplay disp2(img,"change");
    while (!disp1.is_closed() && !disp2.is_closed()) {
        disp1.wait();
        disp2.wait();
    }
}

算法效果：

邻域内直方图统计算法：

#include <iostream>
#include "CImg.h"
using namespace std;
using namespace cimg_library;

float getMean(const float *c){
    float m = 0;
    for(size_t x = 0;x<256;x++){
        m+=x*c[x];
    }
    return m;
}

float getVariance(const float *c,float mean){
    float delta = 0;
    for(size_t y = 0;y<256;y++){
        if(c[y] != 0){
            delta += (pow((y-mean),2)*c[y]);
        }
    }
    delta = sqrt(delta);
    return delta;
}

void getCount(float *c,int originx,int originy,int wid,int hei,CImg<unsigned char> &img){
    for(size_t t = 0;t<256;t++){
        c[t] = 0;
    }
    for(int i = originx;i<wid;i++){
        for(int j = originy;j<hei;j++){
            unsigned char sub = img(i,j);
            c[sub] += 1.0/256.0;
        }
    }
}

void init(float * count){
    for(size_t t = 0;t<256;t++){
        count[t] = 0;
    }
}
int histStatistic(CImg<unsigned char> &img,int dim,float k0,float k1,float k2,float E){
    if(dim%2 == 0){
        cout<<"please input a odd num!!!"<<endl;
        return -1;
    }
    int width = img.width();
    int height = img.height();
    int sum = width * height;
    float count[256];
    for(size_t t = 0;t<256;t++){
        count[t]=0;
    }
    unsigned char sub = 0;
    for(int i = 0;i<width;i++){
        for(int j = 0;j<height;j++){
            sub = img(i,j);
            count[sub]+=1.0/sum;
        }
    }
    float mean = getMean(count);
    float delta = getVariance(count,mean);
    float mxy = 0;
    float deltaxy = 0;
    float local = dim*dim;
    for(int i = dim/2;i<height - dim/2 -1;i++){
        for(int j = dim/2;j<width - dim/2 -1;j++){
            //重新初始化统计数组
            init(count);
            //统计局部直方图
            //getCount(count,j-dim/2,i-dim/2,j+dim/2+1,i+dim/2+1,img);
            for(int p = j-dim/2;p<j+dim/2+1;p++){
                for(int q = i - dim/2;q < i+dim/2+1;q++){
                    sub = img(p,q);
                    count[sub] += 1.0/local;
                }
            }
            mxy = getMean(count);
            deltaxy = getVariance(count,mxy);
        //    cout<<mxy<<"  :  "<<deltaxy<<endl;
            //增强判断
            if(mxy <= mean*k0 && deltaxy <=k2*delta && deltaxy >= k1*delta){
                img(j,i) = img(j,i)*E;
            //    cout<<"enforcement_________________"<<endl;
            }
        }
    }
    return 0;
}

int main(int argc,char ** argv){
    CImg<unsigned char> img("E:\statistic.jpeg");
    CImg<unsigned char> dist(img);
    histStatistic(dist,3,0.4,0.008,0.2,5.0);
    CImgDisplay disp1(img,"origin");
    CImgDisplay disp2(dist,"dist");
    while(!disp1.is_closed() && !disp2.is_closed()){
        disp2.wait();
    }
}

显示效果：

 调节参数后效果出中仍有白色高亮点，类似椒盐噪声，究其原因是由于本地原图操作造成，如果统计后在另一幅目标图上做修改，效果会比较好，以下是异地操作效果：

 最后，用邻域内直方图统计算法对局部直方图算法中的图像做变换，在参数设置为 邻域为7*7，k0=0.4,k1=0,k2 = 0.4,E =20时的增强效果要好于使用局部直方图,以下是两图的比较：