直方图绘图与直方图的均衡化、规定化

1.使用imhist()函数求图片直方图的时候，如果是RGB彩色图，则要首先调用rgb2gray()函数将其转化为灰度图。

eg:

ImageData=imread('baby.jpg');
I=rgb2gray(ImageData );
figure(1);
subplot(2,1,1);
imshow(ImageData);
subplot(2,1,2);
imshow(I);
figure(2);
imhist(uint8(I) );

执行结果如下：

2.bar可以用作绘制条状图、stem()函数绘制杆状图、plot()函数用来绘制图像最常见。

{

linspace()函数：

令A=linspace(1,100) 回车：然后会看到结果生成的是1到100之间的整数，一共100个数字，我们可以看到默认情况下linspace(a1,a2) 是生成包括a1 a2在内的等差数组。

linspace(a1,a2，N) 此函数是用来生成a1与a2之间等距的数组，间距d=(a2-a1)/(N-1)。例如B=linspace(0,9,9)我们可以看到结果如下： 

B =

         0    1.1250    2.2500    3.3750    4.5000    5.6250    6.7500    7.8750    9.0000

因为间距d = 1.125=(9-0)/(9-1)。

}

{

stem()函数：

stem(horz, z, 'fill')中，'fill'参数用来指示杆状图的标记点的颜色和线条的颜色，fill指示默认的‘蓝色，实线，圆’。如果换成‘r--p’就变成了'红色，虚线，五角星形'。详细的参数表格见《数字图像处理(MATLAB版)》P35。

}

ImageData=imread('baby.jpg');
f=rgb2gray(ImageData );
figure(1);
imhist(f);                  %最原始的结果

figure(2);
h = imhist(f,25);
horz = linspace(0,255,25);
bar(horz,h)                 %条形统计图
axis([0 255 0 60000] )
set(gca,'xtick',0:50:255)
set(gca,'ytick',0:20000:60000)

figure(3);
stem(horz,h,'fill');        %杆状图
axis([0 255 0 60000])
set(gca,'xtick',0:50:255)
set(gca,'ytick',0:20000:60000)

figure(4);
plot(imhist(f));            
axis([0 255 0 15000])
set(gca,'xtick',0:50:255)
set(gca,'ytick',0:20000:60000)

结果如下图：

 3.直方图均衡。g= histeq(f,nlev),f为输入图像，nlev为输出图像设定的灰度级数，默认为64。

%{
ImageData=imread('test.tif');
f=rgb2gray(ImageData );
%}
clear,clc;
close all;
f = imread('test.tif');  
figure,imhist(f);                  %最原始的结果
title('原图像直方图');
ylim('auto');
g = histeq(f,256);
figure,imshow(g) 
title('均衡化处理后的图像');
figure,imhist(g)
title('均衡化后的直方图');
ylim('auto');

hnorm = imhist(f)./numel(f);
cdf = cumsum(hnorm);
x = linspace(0,1,256);
figure,plot(x,cdf);
title('变换函数');
axis([0 1 0 1]);
set(gca,'xtick',0:.2:1)
set(gca,'ytick',0:.2:1)
xlabel('Input intensity values','fontsize',9);
ylabel('Output intensity values','fontsize',9);

执行后结果为：

 

4.直方图规定化（自定义规定函数）

f = imread('test.tif'); %读入原图像
p = manualhist;         %输入规定直方图
g= histeq(f,p);         %调用histeq函数得到得到规定化直方图
figure,imhist(g);       %显示规定化后的直方图

如代码中所示，manualhist函数用来得到规定化直方图，并画出图像。 twomodegauss用来计算一个归一化到单位区域的双模态高斯函数。

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).  

% 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

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

执行结果为：

5.函数adapthisteq()：自适应直方图均衡，这种方法用直方图匹配的方法来逐个处理图像中的较小区域。然后用双线性内插的方法来逐个将相邻的小片组合起来，从而消除人工引入的边界。特别是在均匀灰度区域，可以限制对比度来避免放大噪声。adapthisteq的语法是g = adapthisteq(f,param1,val1,param2,val2,...);

eg:

clear,clc;
close all;
f = imread('test.tif'); %读入原图像
imhist(f);

g1 = adapthisteq(f);
figure,imhist(g1);

g2 = adapthisteq(f,'NumTiles',[25,25], 'ClipLimit',0.05 );
figure,imhist(g2);

得到的结果如下：很明显，均衡化处理的结果越来越好：