给出直方图均衡化的数学推导步骤并用程序实现-Read

给出直方图均衡化的数学推导、步骤，并用程序实现。

解： 总体思想：首先考虑连续函数并且让变量r 代表待增强图像的灰度级，假设r 被归一化到区间]1,0[，且0=r 表示黑色及1=r 表示白色。然后再考虑一个离散公式并允许像素值在区间]1,0[-L 内。

对于连续函数而言，假设其变换函数为 )(r T s = 10≤≤r （公式1） 在原始图像中，对于每一个像素值r 产生一个灰度值s 。其中，变换函数要满足以下条件:

（1）)(r T 在区间10≤≤r 中为单值且单调递增。这是为了保证其逆函数的存在，并且输出图像从黑到白顺序增加；

（2）当10≤≤r 时，1)(0≤≤r T 。这保证输出灰度级与输入有同样的范围。

把公式1的逆函数表示为 )(1s T r -= 10≤≤s （公式2）

一幅图像的灰度级可被视为区间]1,0[的随机变量。随机变量的一个最重要的基本描述是其概率密度函数。令)(r P r 和)(s P s 分别代表随机变量r 和s 的概率密度函数。此处带有下标的)(r P r 和)(s P s 用于表示不同的函数。由基本概率理论得到一个基本结果:如果)(r P r 和)(r T 已知，且)(1s T -满足条件(1)，那么变换变量s 的概率密度函数)(s P s 可由以下简单公式得到： ||)()(ds

dr r P s P r s = （公式3） 因此，变换变量s 的概率密度函数由输人图像的灰度级概率密度函数和所选择的变换函数所决定。

在图像处理中一个尤为重要的变换函数: ?

==r

r dw w P r T s 0)()( (公式4)

该被积函数其值为正，并且函数积分是一个函数曲线下的面积，其内含为随机变量r 的累积分布函数，所以它遵守该变换函数是单值单调增加的条件，因此满足条件（1）。同样地，区间]1,0[也满足条件(2)。其积分过程如下：

)()()(0r P dw w P dr d dr r dT dr ds r r r =??

????==? （式5） 用这个结果代替ds

dr ，代入式4,取概率为正，得到： 1)

(1)()()(===r P r P ds dr r P s P r r r s 10≤≤s （式6）

因为)(s P s 是概率密度函数，在这里可以得出区间]1,0[以外它的值为0,可见式6中给出的)(s P s 形式为均匀概率密度函数。换句话来说，公式4给出的变换函数会得到一个随机变量，其特征为一个均匀概率密度函数，与)(r P r 的函数形式是无关的。总述以上，可以看出

?==r

r dw w P r T s 0

)()(便是一个直方图均衡化的基本原理，

该等式右边的意义就是随机变量r 的累积分布函数。这样便转化为了求输入图像灰度级r 的累积分布函数。

下面开始讨论离散函数。对于离散值，处理的是它函数概率的和，而不是概率密度函数的积分。一幅图像中灰度级k r 出现的概率近似为： n

n r p k k k =)( 1,....1,0-=L k (式7) 其中， n 是图像中像素的总和，k n 是灰度级为k r 的像素个数，L 为图像中可能的灰度级总数。式4中变换函数的离散形式为:

∑∑=====k j j k j j r k k n n r p r T s 00)()( 1,....1,0-=L k （式8） 因此，已处理的图像(即输出图像)由通过式8，将输入图像中灰度级为k r 的各像素映射到输出图像中灰度级为k s 的对应像素得到。与连续形式不同，一般不能证明离散变换能产生均匀概率密度函数的离散值(为均匀直方图)。但是不论怎样，可以很容易地看出，式8的应用有展开输入图像直方图的一般趋势，以至于直方图均衡化过的图像灰度级能跨越更大的范围。至此，便给出了整个的证明过程。

下面列出了直方图均衡化的一般实现过程：

（1） 统计原始输入图像各灰度级的像素数目i n ，1...,1,0-=L i ,其中L 为灰度总级

数；

（2） 计算原始图像直方图，即各灰度级的概率密度，n n r P i i i =

)(，n 为原始图像的总像素数目；

（3） 计算累积分布函数

∑=≈k

i i

i k k r P r s 0)()( 1,....1,0-=L k （式9） （4） 计算最后的输出灰度级，

()[])1/(5.0)(min min max -++-=L g r s g g INT g k k k 1,....1,0-=L k (式10) 式中[]INT 是取整算符。令0min =g ，1max -=L g ，则计算式简化为

()[])1/(5.0)(1-+-=L r s L INT g k k k (式11)

（5） 用k f （原图像的灰度级函数）和k g 的映射关系，修改原图像的灰度级，获得

输出图像，其直方图为近似均匀分布。

下面将会给出相关的程序实现代码（使用MATLAB 来实现）。

%在操作目录下面有一个2.jpg 的图片

ima='2.jpg';

info=imfinfo(ima);

if info.ColorType=='truecolor'

source1=imread(ima);

source=rgb2gray(source1);

imwrite(source,'temp.jpg');

source=imread('temp.jpg');

info=imfinfo('temp.jpg');

elseif info.ColorType=='grayscale'

source=imread(ima);

end

figure(1),imshow(source)

figure(2),imhist(source)

%下面用来进行直方图均衡化处理

info=imfinfo(ima);

L=2^info.BitDepth;

recordin=[0,0]; %原图像的概率密度统计

sk=[0,0]; %原图像的累积函数

recordout=[0,0];%变换后图像的概率密度统计

change=[0,0]; %转换函数

for i=3:1:L %初始化空间

recordin(i)=0;

recordout(i)=0;

sk(i)=0;

change(i)=0;

end

for i=1:1:info.Height

for j=1:1:info.Width

recordin(source(i,j)+1)=recordin(source(i,j)+1)+1; %

统计原图灰度的概率分布 end

end

n=info.Height*info.Width; %计算像素总值

%下面的方法亦可以用来计算n

% n=0;

% for i=1:1:L

% n=n+recordin(i);

% end

for i=1:1:L

recordin(i)=recordin(i)/n; %归一化处理

end

sk(1)=recordin(1);

for i=2:1:L

sk(i)=sk(i-1)+recordin(i); %计算原图像的累积函数

end

%求出转换函数

for i=1:1:L

t=floor((L-1)*sk(i)+0.5);

change(i)=t; %记录转换函数的值

end

%求输出图像的距阵,进行原图像的直方图均匀化处理

output=zeros(size(source));

for i=1:1:info.Height

for j=1:1:info.Width

output(i,j)=change(source(i,j)+1);

end

end

imwrite(output,'out.jpg');

%求输出图像的概率密度统计

for i=1:1:info.Height

for j=1:1:info.Width

recordout(output(i,j)+1)=recordout(output(i,j)+1)+1; %输出图像的概率密度统计 end

end

%用于检查输出图像像素的正确性

%j=recordout(1);

%for i=2:1:L

% j=j+recordout(i);

%end

%for i=1:1:L

% recordout(i)=recordout(i)/n; %归一化处理

%end

I=double(output)/255;

figure(3),imagesc(I,[0 1])

colormap(gray);

figure(4),imhist(I)

t=1:1:L;

figure(5),plot(t,change,t,t,t,recordout，'black')

最后结果如下图：

图1 原先的图片图2 处理后的图片

图3 原图（见图1）的直方图

图4 变换后输出图片的直方图

比较图3和图4,可以看出图像变换前后的效果。为了能够更加清楚的看清楚，我们可以看下面的图并且加以说明。

图5 不同结果的比较

红色线代表未变化直线，蓝色线表示变换函数曲线，黑色线表示由程序统计出来的输出图像概率密度曲线。红色线与蓝色线偏离度的比较，可以通过不同图像的不同矩阵来观察出来；对于不同的图像的矩阵，其变化是不同的，差别相对而言也是较大的。黑色曲线主要是为了与图4进行验证性比较，可以看出用两种函数所做出来的直方图是相同的。下面是用同样的方法所做出来的另一组图像。

图6 变换前的图片图7 变换后的图片