直方图均衡化-优化

直方图匹配

直方图匹配和直方图均衡化类似，但直方图均衡试图使输出图像具有一个平坦的直方图，而直方图匹配是为了得到一个特定形状的直方图。

我们先看一个例子。

>> I=imread('E:\cv\冈萨雷斯数字图像处理MATLAB版图片\dipum_images_ch03\Fig0316(a)(moon).tif');

>> imshow(I)

 

>> figure,imhist(I)

 

>> I2=histeq(I,256);

>> figure,imshow(I2)

 

>> figure,imhist(I2)

 

我们可以看到，经过直方图均衡化效果并不好，变换后的直方图并不平坦，而是缺少了很多低灰度级，视觉上只是明显变亮，并没有达到预期的效果。（造成这种情况的原因是0灰度级太多，统计概率太大，所以直接映射到了高灰度级）。

在这种情况下，一种改善这种效果的方法是直方图匹配。接下来我们先看一下直方图匹配的基础知识。

考虑一下归一到[0,1]区间的连续灰度级，令r和z分别表示输入图像和输出图像的灰度级。输入图像的灰度级具有概率密度函数，输出图像的灰度级具有特定的概率密度函数。

我们在之前的讨论中，知道

将会产生一个理想的均衡化的直方图。

现在我们定义一个随机变量z，使得它具有特征

我们可以发现，

我们可以求得T(r)，现在只要我们可以求得，就可以求得具有概率密度函数的随机变量z。

为什么经过这种变换得到的随机变量z的概率密度函数是呢？我们可以证明之。

因为

所以

 

Matlab里的直方图匹配函数

g=histeq(f,hspec)

f表示输入图像，hspec是特定的直方图（一个具有特定值的行向量），g是输出图像，它的直方图近似特定的直方图hspec。histeq的一个特征是，当hspec的长度远小于图像f的灰度级个数，g的直方图可以更好的匹配。

我们已经知道了histeq的基本用法，接下来我们就用histeq函数来实现我们上述问题的直方图匹配问题。

我们渴望得到的直方图，可以在低灰度级端有较小的几种，又可以保持原始图像直方图的基本形状，我们观察到原始图像的直方图是双峰的，在原点有一个大的峰，在高灰度级端有一个小的峰，这种类型的直方图可以通过多模式高斯函数建模。

 

我们使用双峰高斯函数建模，

matlab函数twomodegauss为：

function p = twomodegauss(m1, sig1, m2, sig2, A1, A2, k)

%TWOMODEGAUSS Generates a two-mode Gaussian function.

%   P = TWOMODEGAUSS(M1, SIG1, M2, SIG2, A1, A2, K) generates a

%   two-mode, Gaussian-like function in the interval [0,1].  P is a

%   256-element vector normalized so that SUM(P) equals 1.  The mean

%   and standard deviation of the modes are (M1, SIG1) and (M2,

%   SIG2), respectively. A1 and A2 are the amplitude values of the

%   two modes.  Since the output is normalized, only the relative

%   values of A1 and A2 are important.  K is an offset value that

%   raises the "floor" of the function.  A good set of values to try

%   is M1=0.15, S1=0.05, M2=0.75, S2=0.05, A1=1, A2=0.07, and

%   K=0.002.

 

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins

%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004

%   $Revision: 1.6 $  $Date: 2003/10/13 00:54:47 $

 

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(:));

 

另外需要输入函数的参数，我们这里使用交互式输入函数manualhist：

function p = manualhist

%MANUALHIST Generates a two-mode histogram interactively.

%   P = MANUALHIST generates a two-mode histogram using

%   TWOMODEGAUSS(m1, sig1, m2, sig2, A1, A2, k).  m1 and m2 are the

%   means of the two modes and must be in the range [0,1].  sig1 and

%   sig2 are the standard deviations of the two modes.  A1 and A2 are

%   amplitude values, and k is an offset value that raised the

%   "floor" of histogram.  The number of elements in the histogram

%   vector P is 256 and sum(P) is normalized to 1.  MANUALHIST

%   repeatedly prompts for the parameters and plots the resulting

%   histogram until the user types an 'x' to quit, and then it returns

%   the last histogram computed.

%

%   A good set of starting values is: (0.15, 0.05, 0.75, 0.05, 1,

%   0.07, 0.002).

 

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins

%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004

%   $Revision: 1.7 $  $Date: 2003/10/13 00:49:57 $

 

% Initialize.

repeats = true;

quitnow = 'x';

 

% Compute a default histogram in case the user quits before

% estimating at least one histogram.

%p = twomodegauss(0.15, 0.05, 0.75, 0.05, 1, 0.07, 0.002);

 

% Cycle until an x is input.

while repeats

   s = input('Enter m1, sig1, m2, sig2, A1, A2, k OR x to quit:','s');

   if s == quitnow

      break

   end

 

   % Convert the input string to a vector of numerical values and

   % verify the number of inputs.

   v = str2num(s);

   if numel(v) ~= 7

      disp('Incorrect number of inputs')

      continue

   end

 

   p = twomodegauss(v(1), v(2), v(3), v(4), v(5), v(6), v(7));

   % Start a new figure and scale the axes. Specifying only xlim

   % leaves ylim on auto.

   figure, plot(p)

   xlim([0 255])

end

首先选择合适的直方图，

>> I=imread('E:\cv\冈萨雷斯数字图像处理MATLAB版图片\dipum_images_ch03\Fig0310(a)(Moon Phobos).tif');

>> g=histeq(I,manualhist);

Enter m1, sig1, m2, sig2, A1, A2, k OR x to quit:0.15 0.05 0.75 0.05 1 0.07 0.002

Enter m1, sig1, m2, sig2, A1, A2, k OR x to quit:x

 

>> figure,imhist(g)

 

>> figure,imhist(I)

 

>> figure,imshow(I)

>> figure,imshow(g)

在这里，利用直方图匹配改善了图像的视觉效果，相比直方图均衡，使得最低灰度级没有那么高。

最后，我们修改一下双模式高斯函数的参数，使得k=0.002变为0.02，结果如下：

>> I=imread('E:\cv\冈萨雷斯数字图像处理MATLAB版图片\dipum_images_ch03\Fig0310(a)(Moon Phobos).tif');

>> imshow(I)

 

>> figure,imhist(I)

 

>> I2=histeq(I,manualhist);

Enter m1, sig1, m2, sig2, A1, A2, k OR x to quit:0.15 0.05 0.75 0.05 1 0.07 0.02

Enter m1, sig1, m2, sig2, A1, A2, k OR x to quit:x

 

>> figure,imhist(g)

 

>> figure,imshow(g)

运行时可能出错的地方：

首先，对网上原始图片存储后，要将图片转换为二维图像，运行程序后，与文章中的会有些许差异；

其次，检查两个调用函数中是否存在小错误；

在命令窗口交互式输入命令语句时没有错误，若在编辑框中直接运行，Enter语句加注释即可，只需在命令行中输入相应参数。

我们发现，仅改变了参数k，使得匹配后的直方图有这么大的改变，导致发生这么大的变化，有待探究。但我们可以知道，双混合高斯的建模要注意。