数字图像处理第二章

灰度变换与空间滤波

本章主要讨论两种重要的空间域处理方法：亮度(或灰度)变换与空间滤波。空间域技术直接对图像的像素进行操作，高度直观。

（一）背景知识

1. 空间域处理由下列表达式表示：
2. g (x, y)= T [ f (x， y)]
3. 其中，f(x,y)为输入图像，g(x，y)为输出(处理后的)图像，T是对图像f的算子，作用于点(x,y)定义的邻域。

（二）灰度变换函数

1.由于输出值仅取决于点的灰度值，而不是取决于点的邻域，因此灰度变换函数通常写成如下形式：
s=T( r）
其中，r表示图像f中的灰度，s表示图像g中的灰度。两者在图像中处于相同的坐标(x,y)处。

2.1 imadjust函数

imadjust 函数是针对灰度图像进行灰度变换的基本图像处理工具箱函数，语法如下：
4. g = imadjust(f,[low_in high_in],[low—out high_out},gamma)

此函数将 f 的灰度值映像到 g 中的新值，也就是将 low_in 与high in 之间的值映射到 low out 与 high out 之间的值。low in 以下与 high in 以上的值可以被截去。也就是将 low_in 以下的值映射为 low_out; 将 high in 以上的值映射为high_out。
参数 gamma 指明了由 f 映射生成图像 g 时曲线的形状。如果 gamma 的值小于 1，映射被加权至较高(较亮)的输出值。如果 gamma 的值大于 1，映射加权至较低(较暗)的输出值。如果省略函数参量，gamma 默认为 1(线性映射)。
代码案例：

f = imread('dog.png')        %读取图片
g1 = imadjust(f,[0 1],[1 0]);&emsp;&emsp;&emsp;%high_out<low_out,灰度反转得到负片
imshow(g1);                      %显示图像g1
g2 = imadjust(f, [0.5 0.75], [0 1]);&emsp;%截取0.5到0.75之间的灰度。
g3 = imadjust (f , [0.5 0.75 ] , [1 0]) ; &emsp;%截取0.5到0.75之间的灰度并取片负片。
g4= imadjust(f,  [0 1]，[0.5 0.75]); &emsp;%将[0-1]压缩至[0.5-0.75]。
g5 = imadjust(f, [], [],2);         %非线性变化处理
 g6= imadjust(f, [], [],0.6);         %非线性变化处理

             原图

               g1 g2 g3
               g4 g5 g6
小结：
g1是负片变换，即取反变换。对一个灰度图像f 取反（负片）变换即为255（2^8）减去f 中的每个像素点的灰度值，得到的新图像g1 即为负片变换后的图像。负片图像的这一过程对于增强在一大片主要的黑色区域中嵌入白色及灰色细节是非常有用的。常用于X光片。

g2图是将[0.5 0.75]之间的灰度扩大到[0 1]整个范围，即仅保留[128-192]灰度的像素，这种类型的处理对于强调图像中感兴趣的灰度区非常有用，能提高图片的对比度。
g3是g2图的负片。
g4图是将[0 1]之间的灰度压缩到[0.5 0.75]范围，图像整体变灰暗了，这种类型的处理能降低图片的对比度。
g5图灰度图像的靠近low_in的灰度值较低像素点灰度值变低，其灰度变化范围被压缩，灰度值靠近high_in的一端的灰度变化范围被拉伸，映射加权至较暗的输出值。
g6图灰度图像靠近low_in的灰度值较低像素点灰度值变高，其灰度变化范围被拉伸，灰度值靠近high_in的一端灰度变化范围被压缩，图像整体变明亮。
如果 gamma 的值小于 1，映射被加权至较高(较亮)的输出值（g6）。如果 gamma 的值大于 1，映射加权至较低(较暗)的输出值（g5）。

imadjust函数中的参数[0-1],实际取决于图像的类型，如果是unit8的图像，则实际区间为2⁸，如果是unit16图像，则实际区间为2¹⁶。区间映射关系为直接映射。

2.2 scretchlim函数

当给imadjust函数传参时，使用scretchlim函数能够用 imadjust 函数自动地而不必关心上面讨论的低高参数的处理。
语法形式：

1. Low_High = stretchlim(f,tol);

示例代码：

 f = imread(‘dog.png’)    %读取图片
 g7 = imadjust(f, stretchlim(f), [ ]);
 imshow(g7)

               左f，右g7

2.3 对数及对比度扩展变换

对数及对比度扩展变换是动态范围处理的基本工具。对数变换通过下列表达式实现:

5. g = c*log(1+(f))

其中，c 是常数，f 是浮点数。在两个标量值中，低值置为 0, 高值置为 1。然而，注意 gamma 曲线是可变的，但是对数函数的形状是固定的。扩展变换实现图像灰度扩展和压缩功能，扩展低灰度值而压缩高灰度值。例如傅立叶频谱的范围在[0,10^6]或更高范围是常见的，而显示器是无法完整的显示如此大范围的灰度值的，从而导致频谱中低灰度值细节会被丢失掉。而将得到的频谱值进行对数变换，可以将其动态范围变换到一个合适区间，通常期望得到的是将压缩值返回至显示的全域。
示例代码：

1. f = imread(‘dog.png’)   %读取图片
2. p1=im2uint8(mat2gray(f));   % mat2gray 将值限定在[0 , 1]范围内，使用 im2uint8 将值限定在[0，255]范围内。把图像转换为 uint8 类。
3. p2 = im2uint8(mat2gray(log(1 + double(f))));
   
              从上至下f, p1, p2

小结：由自然对数的图像可知，横坐标为原图f值，纵坐标为处理后的值，那么低灰度值的值将会提高，图像更亮。其中对数函数中的+1是为了保证输出值不小于0。

2.4.指定任意灰度变换

假设T (1) 是由输入图像的 0 灰度值映射而来的值，T (2)是 由 1 亮度值映射而来的值，依此类推。T(256)是由灰度值 255 映射而来的值。如果我们用值化[0，1]范围内的浮点数来表示输入和输出图像，程序将得到极大简化。其中各灰度值将是浮点数而不是整数。实现灰度映射的简单方法是使用interpl 函数，用法如下：

1.g = interpl(z, T, f)
其中，f 是输入图像，g 是输出图像，T 是刚才说明的列向量，z 是与 T 等长的列向量,形式如下:

1. z = linspace(0, 1, numel(T));

对于 f 中的像素值，interpl 首先查找横坐标(z)的值，然后查找(内插)相应的 T 中的值，
并在相应的像素位置输出 g 的内插值。
示例代码：

 z = 0:10;
 T = sin(z);
y = 0:.25:10;
g = interp1(z,T,y); 
 plot(z,T,'.',y,g)

2.5 针对灰度变换的某些公用M-函数

2.5.1处理可变数目的输入或输出

nargin函数：检测输入到M-函数的参量的数目；

1. n=nargin      %返回输入到M-函数的参量的实际数量

nargout函数：M-函数的输出；

2. n=nargout      %指定调用函数所返回的参数的个数

nargchk函数：能够在M-函数体中检测传递过来的参量的数量是否正确。

3. msg=nargchk（low,high,number）

示例代码：

function spir_len = examp(d, n, lcolor)      %计算螺旋线的周长，并以lcolor颜色填充螺旋线
% 输入参数： 
%   d: 螺旋的旋距
%   n: 螺旋的圈数
%   lcolor：画图线的颜色
% 输出参数：
%   spir_len：螺旋的周长 
if nargin > 3          % 输入参数的实际数量大于3
    error('输入变量过多！');
elseif nargin == 2         % 输入参数的实际数量等于2
    lcolor = 'b';       % 默认情况下为蓝色
end
j = sqrt(-1); 
phi = 0 : pi/1000 : n*2*pi;       %角度
amp = 0 : d/2000 : n*d;        %模长
spir = amp .* exp(j*phi);        %用复数形式表示
 if nargout == 1
    spir_len = sum(abs(diff(spir)));     %diff一阶查分求解螺旋线的周长
    fill(real(spir), imag(spir), lcolor); 
elseif nargout == 0
    plot(spir, lcolor);     %以spir元素为横坐标值， lcolor元素为纵坐标值绘制曲线。
else
    error('输出变量过多！'); 
end
 
axis('square');     %将当前坐标系图形设置为方形。横轴及纵轴比例是1：1

测试：

examp(2,3)

examp(3,5)

examp(1,2,3,4)

nargchk函数示例代码：

function msg = nargchk(low, high, number)
msg = '输入正确';
if (number < low)     %当函数在number的值小于low返回下列消息
msg = '输入参数太少 ' ;
elseif (number > high)      %当函数在number的值大于high返回下列消息
msg = ' 输入参数太多' ;
end

测试：

nargchk(1,2,3)

ans =

' 输入参数太多'

nargchk(1,3,2)

ans =

'输入正确'

nargchk(3,2,1)

ans =

'输入参数太少 '

2.5.2 另一个针对灰度变换的M-函数
intrans函数：能执行负片变换，对数变换，gamma变换及对比度扩展。

• 待办

（三）直方图处理与函数绘图

以从图像灰度直方图中提取信息为基础的灰度变换函数在诸如增强、压缩、分割、描述等方面的图像处理中起重要作用。

3.1生成并绘制图像的直方图

在[0,G]范围内总共有 L 级可能灰度的一幅数字图像的直方图定义为下列离散函数：
h(rk)=nk
其中，$r_k$ 是[0，G]间隔内第k级灰度，$n_k$为图像中出现$r_k$这种灰度的像素数。对于 uint8 类， G的值为 255; 对于 uintl6 类，G的值为 65535;对于浮点图像类，G的值为 1.0。注意,对于 uint8 和 uintl6 类的图像，G =L-1 。有时，利用归一化的直方图是必要的。用 h($r_k$)的所有元素除以图像中的像素总数，就可以简单地得到归一化直方图：

p($r_k$) = h($r_k$)/n = $n_k$/n

其中，对整幅图像k=0, 1, 2, … , L-1。从基础概率论的角度讲，我们认可用p($r_k$)表示灰度级$r_k$出现的概率。

工具箱中用于处理图像直方图的核心函数是imhist,基本语法如下：

h = imhist ( f, b)  

其中的参数b默认值是256，参数b的含义是将灰度级数区间等分为b份，如b= 4，那么区间将分为四份0-63 64-127 128-191 192-255。

通过下列表达式就可以得到归一化的直方图：

p = imhist(f, b) /numel(f)      %numel (f) 函数可以给出数组f中元素的个数 (也就是图像中的像素数)

绘图

1 只有传入图像f一个参数的imhist函数。
示例代码：

f = imread(‘zc.jpg’);
imhist(f)

分析：可以观察到直方图的横坐标为256且等分为256份，因为没传参数b，所已默认值为256。下面对比传参数b= 64的图片。

2 用条形图绘制直方图

bar(horz, z, width)      

其中，z 为包含被绘制点的行向量；horz 为与 z 同维数的向量，包含了水平标度值的增量;width 是一个介于 0和 1 之间的数0 换句话说，horz 的值给出了水平增量，z 的值对应垂直量。 若 horz 被省略，水平轴会从0 至 length (z) 等分为若干个单位。当 width 的值为 1 时，竖条较明显；当 width 的值为 0 时，竖条是垂直线，默认值为 0.8。当绘制条形图时，通常会将水平轴等分为几段以降低分辨率。
示例代码：

>> f= imread('zc.jpg');
>> h = imhist(f,25);                          %等分为25份
>> horz = linspace(0,255,25);              %水平增量
>> bar(horz,h,0.5)                           %绘制直方图				
 >> axis([0 255 0 60000])				%直方图总大小							
>> set(gca,'xtick',0:50:255)				%设置坐标轴
>> set(gca,'ytick',0:20000:60000)

结果如下：

分析：imhist函数与bar函数绘图其实是相同的原理，imhist函数传入参数b为10（近似分为25份），那么数值上会与上图的结果相同。单bar函数绘图更直观。

3 用杆状图绘制直方图

 stem(horz,h, 'r--p') 

示例代码:

f = imread(‘zc.jpg’);
h = imhist(f,25);
horz = linspace(0,255,25);
stem(horz,h,‘fill’) %未设置坐标轴

结果如下：

4 直线绘制直方图

 plot(horz, z, 'LineSpec')     %horz由所有输入点坐标的x值组成，z是由与horz中包含的x对应的y所组成的向量。LineSpec是用户指定的绘图样式。

示例代码：

>> f = imread('zc.jpg');
>> hc = imhist(f);
>> plot(hc)
>> axis([0 255 0 60000])
>> set(gca,'xtick',0:50:255)
>> set(gca,'ytick',0:20000:60000)

3.2直方图均衡化

假设灰度级为归一化在[0,1]范围内的连续量，让$p_r$( r )代表一幅给定图像的灰度级的概率密度函数(PDF)，下标用来区分输入图像和输出图像的概率密度函数。假设我们对输入灰度进行下列变换，得到输出(处理后的)灰度级s:

式中的 w 是积分虚变量。可以看出，输出灰度级的概率密度函数是均匀的，也就是:

当灰度级为离散值时，我们利用直方图并采用前面介绍的直方图均衡化技术。一般来说， 虽然由于变量的离散特性，处理后的图像直方图也不会完全均匀。对于离散的灰度级，我们采用求和的方式，将均衡化变换成为下列形式：

式中，k=0,1 ， 2， …， L-1, 且$S_k$是输出(处理后的)图像的灰度值，对应输入图像的灰度值为$r_k$。

直方图均衡化由工具箱中的 histeq 函数实现，语法如下:

g = histeq(f,nlev)

其中，f 为输入图像，nlev 为输出图像设定的灰度级数。
示例代码：

>> f = imread('timg.jpg');
>> imhist(f)
>> ylim('auto')    %自动设定纵坐标轴的取值范围和刻度线
>> g = histeq(f,256);   %均衡化图像
>> imshow(g)
>> imhist(g)
>> ylim('auto')    
>> 

左原图，右均衡化之后的图。

上原图的直方图，下为处理后的。
分析：可以从直方图中很明显的看出处理之前的像素的灰度较集中，经过均衡化之后，像素的灰度较为分散，这样便使得图像的对比度增强。

3.3直方图匹配化（规定化）

直方图均衡化通过把输入图像的灰度级扩展到较宽灰度范围来实现图像增强，但这种方法有时并不总能导致成功的结果。当图像中接近0的像素过多，进行直方图均衡的时候，会导致统计概率变大，直接映射到高灰度，直接导致图像整体变量，达不到原来想要的结果，此时就需要使用直方图匹配，让变化后的图像具有特定的直方图。有时候，我们为了使两幅图像的色调保持一致，也可以采用该方法。

这种方法在原理上很简单。考虑归一化之后在[0， 1]区间内的连续灰度级，令 r和z分别表 示输入图像与输出图像的灰度级。输入图像的灰度级有概率密度函数$P_r$( r )，输出图像的灰度级具有规定的概率密度函数$P_z$( z )。变换为：

得到的灰度级s具有均勻的概率密度函数$P_s$( s )。假设定义的变量z具有下列特性:

我们要寻找的是灰度级为 Z 的图像，且具有特定的概率密度$P_z$( z )。由前面的两个等式可得:

可以由输入图像得到$T_r$(这是上面讨论的直方图均衡化变换)，由此得出结论：只要找到 H-1, 就能利用前面的等式得到变换后的灰度级z，概率密度函数(PDF)为指定的$P_z$( z )。
语法如下：

g = histeq(f,hspec)     

其中，f 为输入图像，hspec 为特定的直方图(某个特定值的行向量)，g 为输出图像，直方图近似于指定直方图 hspec。histeq 的特性是当 length ( hspec ) 比图像 f 中的灰度级小很多时，图像 g 的直方图通常会较好地匹配hspec。

示例代码：

>> f = imread('timg.jpg');
>> g = imhist(imread('atg.jpg'));
>> match = histeq(f,g);
>> imshow(match)
>> subplot(2,3,1),imshow(f),title('原图');
>> subplot(2,3,2),imshow(g),title('模板图');
>> subplot(2,3,3),imshow(match),title('处理后的图');
>> subplot(2,3,4),imhist(f),title('原图直方图');
>> g1 = imread('atg.jpg');
>> subplot(2,3,5),imhist(g1),title('模板图直方图');
>> subplot(2,3,6),imhist(match),title('匹配完成直方图');

分析：均衡化是将灰度集中的像素点分散覆盖到整个范围，而有时这种算法并不能处理到我们想要的图像，而匹配化是将待处理的图按照模板去扩展动态范围，可以根据主观去调整像素灰度级数的分布。从上面的实验可知，匹配之后的直方图是以模板图为模板动态扩展。

3.4 函数 adapthisteq

这个工具箱函数执行所谓的对比度受限的自适应直方图均衡。这个方法是由用直方图规定化方法处理图像的小区域(称为小片)组成。然后用双线性内插将相邻小片组合起来以消除人工引入的边界效应。特别是可以限制均匀亮度区域的对比度， 以免放大噪声。
语法如下：

g = adapthisteq(f, paraml, vail, param2, val2, ...)            %f是输入图像，g是输出图像

示例代码：

>> f= imread('timg.jpg');   
>> g = rgb2gray(f);
>> g = rgb2gray(f);     %g为原图f的灰度图像
>> g1 = adapthisteq(g);
>> g1 = adapthisteq(g);   %g1为默认值
>> g2 = adapthisteq(g,'NumTiles',[25 25]);   %g2将默认值改为25 25
>> g3 = adapthisteq(g,'NumTiles',[25 25],'clipLimit',0.05);   %g3片大小为25 25，对比度0.05
>> subplot(2,5,1),imshow(f),title('原图');
>> subplot(2,5,2),imshow(g),title('灰度图');
>> subplot(2,5,3),imshow(g1),title('默认值图');
>> subplot(2,5,4),imshow(g2),title('片大小25 25');
>> subplot(2,5,5),imshow(g3),title('片大小25 25，对比度0.05');
>> subplot(2,5,6),imhist(f);
>> subplot(2,5,7),imhist(g);
>> subplot(2,5,8),imhist(g1);
>> subplot(2,5,9),imhist(g2);
>> subplot(2,5,10),imhist(g3);

小结：im2double():将图象数组转换成double精度类型，im2uint8():将图象数组转换成unit8类型。线性映射imadjust。对数动态扩展，将低灰度级调高。均衡化，分散灰度集中的像素。匹配化，根据输入模板动态扩展。