直方图均衡化和规定化

图像的直方图事实上就是图像的亮度分布的概率密度函数，是一幅图像的所有象素集合的最基本的统计规律。直方图反映了图像的明暗分布规律，可以通过图像变换进行直方图调整，获得较好的视觉效果。直方图均衡化是通过灰度变换将一幅图像转换为另一幅具有均衡直方图的图像，即在每个灰度级上都具有相同的象素点数的过程。处理后的图像直方图分布更均匀了，图像在每个灰度级上都有像素点。从处理前后的图像可以看出，许多在原始图像中看不清楚的细节在直方图均衡化处理后所得到的图像中都变得十分清晰。所谓直方图规定化，就是通过一个灰度映像函数，将原灰度直方图改造成所希望的直方图。所以，直方图修正的关键就是灰度映像函数。

一、直方图均衡化核心思想

直方图均衡化的核心思想就是通过构造灰度级变化函数来改造原图像的直方图，使得变化后的图像的直方图达到某一要求，最终目的是为了提高图像灰度对比度，提高图像的可视化质量。

例如：原图像灰度直方图如下图所示：

显然该图像亮度偏暗，可视化质量不高，事实也正如此，如下是原图像：

现在需要做的就是增大图像的对比度，那就要增大最高灰度级和最低灰度级之间的差值，因此需要在保证图像可是质量的条件下增大最大亮度和减小最小亮度，因此均衡化成为了比较好的选择。

对比度=最大亮度/最小亮度

    直方图均衡化是通过灰度变换将一幅图像转换为另一幅具有均衡直方图的图像，即在每个灰度级上都具有相同的像素点数的过程。处理后的图像直方图分布更均匀了，图像在每个灰度级上都有像素点。即通过变化函数使得图像的灰度级按照灰度级出现的相对频数近似均匀分布。但是由于“简并现象”的存在，使得均衡化处理后的图像的灰度级不会完全符合均匀分布。

    设r代表要增强的图像中像素的灰度级，s代表增强后新图像中的灰度级，则变换函数就是r经过函数T()映射以后得到的s即为r变换后的灰度级。

累计求和变化函数公式如下：

数字图像君合均衡化处理算法步骤如下：

1. 由统计原始图像直方图；
2. 由式（1.1）计算新的灰度级；
3. 修正为合理的灰度级；
4. 计算新的直方图；
5. 用处理后的新灰度代替处理前的灰度，生成新图像。

经过以上变换之后，得到新的图像和新图像的灰度直方图：

 二、图像规定化算法核心思想

直方图规定化，就是通过一个灰度映像函数，将原灰度直方图改造成所希望的直方图。所以，直方图修正的关键就是灰度映像函数。灰度映像函数可以是一张给定的图片，也可以是某一种分布。

例如：如下是通过一张目标图规定化原图的灰度级分布。

也可以将上图中的匹配图像换成双模态高斯函数规定化原图灰度级分布。直方图规定化就是要使得原图的直方图更像目标图的直方图，如何才能做到更像，最简单的思想就是灰度级rk映射为其灰度级累积相对频数最接近目标图像灰度级累积相对频数所对应的灰度级sk,这种思想所对应的算法就是SML算法。

直方图规定化的算法步骤为：

1. 分别计算出原图像和规定图像的累积直方图
2. 利用SML算法映射灰度值
3. 利用更新后的映射表转换原图像的灰度值

例如：

分别计算出原图像和规定图像的累积直方图：

利用SML算法映射灰度值

0->3是因为原始累积直方图0.19最接近规定累积直方图中的0.15；同理0.44最接近0.35；0.65最接近0.65；0.81和0.89都最接近0.85；0.95，0.98，1都最接近1

因此，按照上图这种映射关系，就可以得到原图关于规定直方图的直方图。

三、直方图均衡化代码

用累积求和函数、sigmoid()函数求解直方图均衡化

%一、自编直方图均衡化代码
Image=imread('D:\课程\数字图像处理\课程设计\image\lena.bmp');
histgram=imhist(Image);%histgram中存储原图像各个灰度级的个数
[h,w,c]=size(Image);%读取图片的大小以及通道信息
NewImage=zeros(h,w);%设置新图片的大小

%1、用累积求和函数求解直方图均衡化
CSF=zeros(256);%定义累积求和函数
CSF(1)=histgram(1);
%模拟累积求和过程
for gray_scale=2:256
    CSF(gray_scale)=CSF(gray_scale-1)+histgram(gray_scale);
end
%得到累积求和函数之后，原灰度级r->新灰度级s的映射关系为：s=CSF(r)/(h*w)
%根据映射关系生成新图片
for x=1:h
    for y=1:w
        NewImage(x,y)=CSF(Image(x,y)+1)/(h*w);
    end
end

%2、除了累积求和函数之外，还可以用sigmoid()函数
img = imread('D:\课程\数字图像处理\课程设计\image\lena.bmp');
%统计原图各灰度级的个数
histgram = imhist(img);
%利用sigmoid将原灰度级映射成新灰度级
new_img_hist = 1 ./ (1 + exp(-(histgram - 0.5)));
%使用histeq函数按照new_img_hist直方图规则对原图进行直方图规定化
new_img = histeq(img, new_img_hist);

%显示均衡化图像
figure;
axis tight;%将坐标轴范围设置为等同于数据范围，使轴框紧密围绕数据
subplot(2,3,1);imshow(Image);title('原始图像');
subplot(2,3,2);imhist(Image);title('原始图像直方图');
subplot(2,3,3);imshow(NewImage);title('累积求和函数均衡化图像');
subplot(2,3,4);imhist(NewImage);title('累积求和函数均衡化图像直方图');
subplot(2,3,5);imshow(new_img);title('sigmoid()函数处理均衡化图像');
subplot(2,3,6);imhist(new_img);title('sigmoid()函数处理均衡化图像直方图');

运行结果

四、直方图规定化代码

%二、直方图规定化
origin=imread('D:\课程\数字图像处理\课程设计\image\lena.bmp');
p=manualhist;
% p=twomodegauss(0.15,0.05,0.75,0.05,1,0.07,0.002);
[m_o,n_o,c_o]=size(origin);

origin_hist=imhist(origin)/(m_o*n_o);
standard_hist=p;

%定义原图累积求和数组、目标分布累积求和数组
origin_sum=[];
standard_sum=[];
%累积求和原图各灰度级个数，累积求和目标分布各灰度级个数
for i=1:256
    origin_sum=[origin_sum sum(origin_hist(1:i))];
    standard_sum=[standard_sum sum(standard_hist(1:i))];
end

%遍历原图各灰度级累计求和数组中的每一个值，假设这个值为a
%将原灰度级映射到a距离standard_sum中最接近值所对应的灰度级，也就是
%standard_sum的下标，用index数组记录这个映射关系，index数组的索引是原灰度级
%index中索引对应的值是目标灰度级
for i=1:256
    value{i}=standard_sum-origin_sum(i);
    value{i}=abs(value{i});
    [temp index(i)]=min(value{i});
end
%根据index中记录的映射关系生成规定化图像
NewImage=zeros(m_o,n_o);
for i=1:m_o
    for j=1:n_o
        NewImage(i,j)=index(origin(i,j)+1)-1;%加1是为了取值，减1是为了表示灰度级
    end
end
NewImage=uint8(NewImage);

%显示规定化图像
figure;
axis tight;
subplot(2,3,1);imshow(origin);title('原始图像');
subplot(2,3,2);imhist(origin);title('原始图像直方图');
subplot(2,3,3);plot(p);title('双模态高斯函数图像');
subplot(2,3,4);imshow(NewImage);title('直方图规定化图像');
subplot(2,3,5);imhist(NewImage);title('直方图规定化图像直方图');

%双模态高斯函数
function p=manualhist
% p=twomodegauss(0.15,0.05,0.75,0.05,1,0.07,0.002);
p=twomodegauss(0.15,0.05,0.75,0.2,1,0.07,0.002);

function p=twomodegauss(m1,sig1,m2,sig2,A1,A2,k)
    c1=A1*(1/((2*pi)^0.5)*sig1);
    k1=2*(sig1^2);
    c2=A2*(1/((2*pi)^0.5)*sig2);
    k2=2*(sig2^2);
    z=linspace(0, 1, 256);
    p=k+c1*exp(-((z-m1).^2)./k1)+c2*exp(-((z-m2).^2)./k2);
    p=p./sum(p(:));
end
end

运行结果