图像变换(傅立叶变换), 图像增强, 边缘检测, 滤波, 图像压缩等.

实验工具:MATLAB软件

课程设计时间:2008年12月

实　验　部　分

1. 图像变换

程序代码及说明

clear all

N=100;

f=zeros(50,50); %产生一个50*50的全0数组

f(15:35,23:28)=1;

%定义图像数组，从15行到35行，23列到28列附值为1，为白色，其他区域为黑色

figure(1) %创建窗口的图形对象，句柄为1

imshow(f,'notruesize') %显示图像f

F=fft2(f,N,N); %在二维傅立叶变换前把f截断或者添加0，使其成为N*N的数组

F2=fftshift(abs(F)); %把傅立叶变换的零频率部分移到频谱的中间

figure(2)

x=1:N;y=1:N; %定义x和y的范围

mesh(x,y,F2(x,y));colormap(gray);colorbar

%绘制立体网状图，将图形对象的色度改为灰度图像，colorbar给坐标轴添加色彩条

%构建一个类似于figure(1)的矩形函数

N=200;

f=zeros(100,100);

f(30:70,45:55)=1;

%定义图像数组，从30行到70行，45列到55列附值为1，为白色，其他区域为黑色

imshow(f,'notruesize');

%然后对f进行二维快速傅立叶变换：以下列出你自己编写的代码…

N=200;

f=zeros(100,100);

f(30:70,45:55)=1;

imshow(f,'notruesize');

title('原始图像');

F=fft2(f,N,N); %对图像f进行二维快速傅立叶变换

grid

on %打开网格线

axis

on %打开坐标轴

imshow(F,[-1,5],'notruesize'); %显示傅立叶变换后的图像，图像数据的值域为[-1,5]

x=1:N;y=1:N;

title('二维快速傅立叶变换后的图像');

mesh(abs(F)); %绘制F的频谱图

title('傅立叶变换后的频谱图');

%然后对上述二维快速傅立叶变换提高分辨率：

要提高二维快速傅立叶变换的分辨率，在采样率一定的情况下，增大采

样点数N即可。对应的频谱图见测试结果。

N=300;

f=zeros(100,100);

f(30:70,45:55)=1;

figure(1)

imshow(f,'notruesize');

title('原始图像');

F=fft2(f,N,N);

axis on

figure(2)

imshow(F,[-1,5],'notruesize');

title('二维快速傅立叶变换后的图像');

x=1:N;y=1:N;

figure(3)

mesh(abs(F));

title('傅立叶变换后的频谱图');

N=400;

f=zeros(100,100);

f(30:70,45:55)=1;

figure(1)

imshow(f,'notruesize');

title('原始图像');

F=fft2(f,N,N);

axis on

figure(2)

imshow(F,[-1,5],'notruesize');

title('二维快速傅立叶变换后的图像');

x=1:N;y=1:N;

figure(3)

mesh(abs(F));

title('傅立叶变换后的频谱图');

�系数移动(使用函数fftshift)

N=200;

f=zeros(100,100);

f(30:70,45:55)=1;

figure(1)

imshow(f,'notruesize');

title('原始图像');

F=fft2(f,N,N);

axis on

figure(2)

imshow(F,[-1,5],'notruesize');

title('二维快速傅立叶变换后的图像');

x=1:N;y=1:N;

figure(3)

mesh(fftshift(abs(F)));

%把零频率部分移到频谱的中间

title('傅立叶变换后的频谱图');

%滤波器频率响应

x=1:N;y=1:N;

mesh(x,y,F2(x,y));colormap(gray);colorbar

%绘制立体网状图，将图形对象的色度改为灰度图像，colorbar给坐标轴添加色彩条

测试结果图：

1》对f进行二维快速傅立叶变换

(请自己运行查看)

2》对上述二维快速傅立叶变换提高分辨率

N=300时： N=400时：

N=800时： N=1000时：

(请自己运行查看)

从傅立叶变换的频谱图中可以看出，提高分辨率以后，其边缘更加平滑，锯齿状明显减弱。但其傅立叶变换后的图像没有明显改变。

3》DC系数移动

其系数移动以后，频谱分量都集中到了频谱的中间。

4》滤波器频率响应

(请自己运行查看)

2. 图像增强

图像增强是指按特定的需要突出一幅图像中的某些信息，同时，消弱或去除某些不需要的信息的处理方法。其主要目的是使处理后的图像对某些特定的应用比原来的图像更加有效。图像增强技术主要包含直方图修改处理、图像平滑化处理、图像尖锐化处理、和彩色处理技术等。图像增强有图像对比度增强、亮度增强，轮廓增强等等。

下面利用直方图统计算法对灰度图像进行增强：

程序代码：

I=imread('cameraman.tif');

subplot(121)

imshow(I);

title('原始图像');

subplot(122)

imhist(I,64)

%绘制图像的直方图，n=64为灰度图像灰度级，若I为灰度图像，默认n=256；

若I为二值图像，默认n=2。

title('图像的直方图');

(请自己运行查看)

n=256时：

(请自己运行查看)

下面利用直方图均衡化增强图像的对比度：

I=imread('cameraman.tif');

J=histeq(I);

%将灰度图像转换成具有64(默认)个离散灰度级的灰度图像

imshow(I)

title('原始图像')

figure,imshow(J)

title('直方图均衡化后的图像')

figure(1)

subplot(121);imhist(I,64)

title('原始图像的直方图')

subplot(122);imhist(J,64)

title('均衡化的直方图')

(请自己运行查看)

分析：从上图中可以看出，用直方图均衡化后，图像的直方图的灰度间隔被拉大了，均衡化的图像的一些细节显示了出来，这有利于图像的分析和识别。直方图均衡化就是通过变换函数histeq将原图的直方图调整为具有“平坦”倾向的直方图，然后用均衡直方图校正图像。

下面利用直方图规定化对图像进行增强：

I=imread('cameraman.tif');

figure,imshow(I);

title('原始图像');

hgram=50:2:250; %规定化函数

J=histeq(I,hgram);

figure,imshow(J);

title('直方图规定化后的图像');

figure,imhist(I,64);

title('原始图像的直方图');

figure,imhist(J,64);

title('直方图规定化后的直方图');

运行结果：

(请自己运行查看)

变换灰度间隔后的图像和直方图：

hgram=50:1:250; hgram=50:5:250;

3. 图像重建

图像重建的最典型的应用是医学上的计算机断层摄影技术(CT技术)。它用于人体头部、腹部等内部器官的无损伤诊断，其基本方法就是根据人体截面投影，经过计算机处理来重建截面图像。在人体中把需要扫描的部分取出一定厚度的断层面，再把断层面分成许多小的方块。当一束较窄的射线通过每个方块后强度就有一定程度的衰减，衰减的量由此方块的分子构成和组织密度决定。如果通过各种角度重复上述过程以获得一系列强度分布曲线，就有可能从这些数据中计算每一方块的衰减量。这样就能够重建断层或三维图像。

目前提出的图像重建方法有以下五种：

(1)联立方程法(也称矩阵法)；

(2)逆投影法；

(3)付立叶变换法；

(4)滤波──逆投影法(也称卷积法)；

(5)逐次逼近法。

滤波──逆投影法是当前用得较多的一种图像重建方法，在当代射线CT系统中几乎都用这种方法构成系统。它的特点是精度高，能快速实现。

在图像处理的工具箱中，MATLAB提供了一个计算图像沿着指定方向上的投影的函数——radon函数。

iradon函数可以实现radon逆变换，radon逆变换通常应用于X线断层摄影术中，可以从投影数据中重构图像。

下面利用radon函数和iradon函数计算图像的投影并从投影中重建图像，将Shepp-Logan的大脑图作为测试图。

函数radon和函数iradon的调用格式：

[R,xp]=radon(I,theta) 计算图像I在theta向量所指定的方向上的radon变换，I表示待处理的图像，theta表示radon变换的方向角度，可以是标量或向量值，返回值

R的每一列对应图像I在theta某一角度的radon变换值，xp向量表示沿着x'轴对应的坐标值。

IR=iradon(R,theta) 利用R各列中投影值来构造图像I的近似值。投影数越多，获得的图像越接近原始图像，角度theta必须是固定增量的均匀向量。

程序代码：

P=phantom(256); %用phantom函数产生Sheep-Logan的大脑图，n为图像p中的行列数，默认为256

imshow(P)

title('原始图像')

%以下为三种不同角度的投影模式

theta1=0:10:170;[R1,xp]=radon(P,theta1); %存在18个角度投影

theta2=0:5:175;[R2,xp]=radon(P,theta2); %存在36个角度投影

theta3=0:2:178;[R3,xp]=radon(P,theta3); %存在90个角度投影

figure,imagesc(theta3,xp,R3);colormap(hot);colorbar;

%显示图像Sheep-Logan的radon变换

title('经radon变换后的图像')

xlabel('\theta');ylabel('x\prime'); %定义坐标轴

%用三种情况的逆radon变换来重建图像

I1=iradon(R1,10);

I2=iradon(R2,5);

I3=iradon(R3,2);

figure,imshow(I1)

title('角度增值为10时的iradon变换图像')

figure,imshow(I2)

title('角度增值为5时的iradon变换图像')

figure,imshow(I3)

title('角度增值为2时的iradon变换图像')

运行结果图：

(请自己运行查看)

由上面重建的图像中可以看出，只用18个投影来重建图像效果很差，而36个投影来重建的图像要好的多，90个投影来重建的图像质量更好，失真也很小，由于R1重建图像的投影太少，所以存在许多虚假点，重建的效果与投影数目相关，投影数目越多图像重建的效果越好，所以要提高重建图像的质量，就需要增加投影角度的数目。

除此之外，还可以在Fan-Beam变换数据中用ifaanbeam函数重建图像。

4. 滤波

4.1 目的

运用中值滤波克服线性滤波器所带来的图像细节模糊。

4.2 使用设备

PC兼容机一台，操作系统为Windows2000(或Windows98，WindowsXP，以下默认为Windows2000)

4.3 使用滤波对图像进行增强

4.3.1

线性滤波(邻域平均)

线性低通滤波器最常用的是线性平滑滤波器，这种滤波器的所有系数都是正的，也称邻域平均。邻域平均减弱或消除了傅立叶变换的高频分量，对噪声的消除有所增强，但是由于平均而使图像变得更为模糊，细节的锐化程度逐渐减弱。

下面使用不同的平滑模板对图像进行滤波：(二维线性滤波fliter2)

程序代码：

I=imread('cameraman.tif');

imshow(I)

title('原始图像')

J=imnoise(I,'salt &

pepper'); %添加盐椒噪声，噪声密度为默认值0.05

figure,imshow(J)

title('添加盐椒噪声后的图像')

K1=filter2(fspecial('average',3),J)/255; %应用3×3邻域窗口法

figure,imshow(K1)

title('3×3窗的邻域平均滤波图像')

K2=filter2(fspecial('average',7),J)/255; %应用7×7邻域窗口法

figure,imshow(K2)

title('7×7窗的邻域平均滤波图像')

K3=filter2(fspecial('average',9),J)/255; %应用9×9邻域窗口法

figure,imshow(K3)

title('9×9窗的邻域平均滤波图像')

K4=filter2(fspecial('average',11),J)/255; %应用11×11邻域窗口法

figure,imshow(K4)

title('11×11窗的邻域平均滤波图像')

运行结果图：

(请自己运行查看)

4.3.2 中值滤波

中值滤波可以保留目标边缘，这是中值滤波器相对于均值滤波器的最大优势。中值滤波具有去噪的性能，可以消除孤立的噪声点，可以用来减弱随机干扰和脉冲干扰，但是边缘不模糊。

程序代码：

I=imread('cameraman.tif');

imshow(I)

title('原始图像')

J=imnoise(I,'salt &

pepper',0.02); %添加盐椒噪声，噪声密度为0.02

figure,imshow(J)

title('添加盐椒噪声后的图像')

K1=medfilt2(J); %在默认的3×3的邻域窗中进行中值滤波

figure,imshow(K1)

title('默认的3×3的邻域窗的中值滤波图像')

K2=medfilt2(J,[5

5]); %在5×5的邻域窗中进行中值滤波

figure,imshow(K2)

title('5×5的邻域窗的中值滤波图像')

运行结果图：

(请自己运行查看)

从上可见，中值滤波的效果要比邻域平均的低通滤波效果好，中值滤波以后的图像的轮廓比较清晰，而且使用较小的模板得到的视觉效果反而好一些。

4.3.3 锐化滤波

图像锐化处理的目的是使模糊图像变得清晰，锐化滤波器减弱或消除了傅立叶空间的低频分量，保留高频分量，从而加强了图像的轮廓，使图像看起来比较清晰。

下面应用Laplacian算子对图像进行锐化处理：

Laplacian算子是线性二次微分算子，其格式为：h = fspecial('laplacian',

alpha)，返回一个3×3的滤波器来近似二维Laplacian算子的形状，参数alpha决定了Laplacian算子的形状，alpha的取值范围为0.0~1.0，默认的值为0.2。

程序代码：

%应用Laplacian算子对图像进行锐化

I=imread('cameraman.tif');

imshow(I)

title('原始图像')

H=fspecial('laplacian');

%应用laplacian算子滤波锐化图像

laplacianH=filter2(H,I);

figure,imshow(laplacianH)

title('laplacian算子锐化后的图像')

运行结果图：

(请自己运行查看)

分析：由图可以看出，应用了Laplacian算子对图像锐化以后，将图像区域的边缘轮廓勾划了出来，因此Laplacian算子对于边缘检测也具有很好的功效。

5. 边缘检测

下面利用sobel算子对图像进行边缘检测：

使用edge函数实现图像的边缘检测，其调用格式为：

BW=edge(I,'sobel',thresh,direction) 根据指定的敏感阈值thresh用Sobel算子对图像进行边缘检测，edge函数忽略了所有小于阈值的边缘，如果没有指定阈值thresh或为空，函数自动选择参数值，direction指定Sobel算子边缘检测的方向，其参数值为'horizontal'，'vertical'或'both'(默认)。

程序代码：

I=imread('cameraman.tif');

imshow(I)

title('原始图像')

BW=edge(I,'sobel');

%以自动域值选择法对图像进行Sobel算子边缘检测

figure,imshow(BW);

title('自动域值的Sobel算子边缘检测')

[BW,thresh]=edge(I,'sobel');

%返回当前Sobel算子边缘检测的阈值

disp('sobel算子自动选择的阈值为：')

disp(thresh)

BW1=edge(I,'sobel',0.02,'horizontal');

%以域值为0.02水平方向对图像进行Sobel算子边缘检测

figure,imshow(BW1)

title('域值为0.02的水平方向的sobel算子检测')

BW2=edge(I,'sobel',0.02,'vertical');

%以域值为0.02垂直方向对图像进行Sobel算子边缘检测

figure,imshow(BW2)

title('域值为0.02的垂直方向的sobel算子检测')

BW3=edge(I,'sobel',0.05,'horizontal');

%以域值为0.05水平方向对图像进行Sobel算子边缘检测

figure,imshow(BW3)

title('域值为0.05的水平方向的sobel算子检测')

BW4=edge(I,'sobel',0.05,'vertical');

%以域值为0.05垂直方向对图像进行Sobel算子边缘检测

figure,imshow(BW4)

title('域值为0.05的垂直方向的sobel算子检测')

测试结果图：

(请自己运行查看)

sobel算子自动选择的阈值为：0.1433

由图可以看出，在采用水平和垂直方向的Sobel算子对图像进行边缘检测时，分别对应的水平和垂直方向上的边缘有较强的响应，阈值越小，检测的图像的边缘细节数越多，而增大阈值时，有些轮廓则未能检测出。

问题与思考：

根据Prewitt算子的定义设计实现用Prewitt算子进行图像的边缘检测。

其用法和Sobel算子类似。其调用格式为：

BW=edge(I,'prewitt',thresh,direction)

根据指定的敏感阈值thresh用Prewitt算子对图像进行边缘检测。

程序代码：

I=imread('cameraman.tif');

imshow(I)

title('原始图像')

BW=edge(I,'prewitt');

%以自动域值选择法对图像进行Prewitt算子边缘检测

figure,imshow(BW);

title('自动域值的prewitt算子边缘检测')

[BW,thresh]=edge(I,'prewitt');

%返回当前Prewitt算子边缘检测的阈值

disp('prewitt算子自动选择的阈值为：')

disp(thresh)

BW1=edge(I,'prewitt',0.02,'horizontal');

%以域值为0.02水平方向对图像进行Prewitt算子边缘检测

figure,imshow(BW1)

title('域值为0.02的水平方向的prewitt算子检测')

BW2=edge(I,'prewitt',0.02,'vertical');

%以域值为0.02垂直方向对图像进行Prewitt算子边缘检测

figure,imshow(BW2)

title('域值为0.02的垂直方向的prewitt算子检测')

BW3=edge(I,'prewitt',0.05,'horizontal');

%以域值为0.05水平方向对图像进行Prewitt算子边缘检测

figure,imshow(BW3)

title('域值为0.05的水平方向的prewitt算子检测')

BW4=edge(I,'prewitt',0.05,'vertical');

%以域值为0.05垂直方向对图像进行Prewitt算子边缘检测

figure,imshow(BW4)

title('域值为0.05的垂直方向的prewitt算子检测')

测试结果图：

(请自己运行查看)

prewitt算子自动选择的阈值为：0.1399

6. 图像压缩

图像压缩就是就是通过去除这些数据冗余来减少表示数据所需的比特数，去除多余数据。以数学的观点来看，这一过程实际上就是将二维像素阵列变换为一个在统计上无关联的数据集合。

图像压缩是指以较少的比特有损或无损地表示原来的像素矩阵的技术，也称图像编码。

图像数据之所以能被压缩，就是因为数据中存在着冗余。图像数据的冗余主要表现为：图像中相邻像素间的相关性引起的空间冗余；图像序列中不同帧之间存在相关性引起的时间冗余；不同彩色平面或频谱带的相关性引起的频谱冗余。

从压缩编码算法原理上可以分类为：

(1)无损压缩编码种类：哈夫曼编码、算术编码、行程编码、Lempel zev编码

(2)有损压缩编码种类：

预测编码：DPCM，运动补偿

频率域方法：正交变换编码(如DCT)，子带编码

空间域方法：统计分块编码

模型方法：分形编码，模型基编码

基于重要性：滤波，子采样，比特分配，矢量量化

(3)混合编码 ?JBIG，H261，JPEG，MPEG等技术标准

利用余弦变换实现图像压缩：

DCT先将整体图像分成N×N像素块(一般N=8

，即64个像素块)，再对N×N块像素逐一进行DCT变换。由于大多数图像高频分量较小，相应于图像高频成分的失真不太敏感，可以用更粗的量化，在保证所要求的图质下，舍弃某些次要信息。

程序代码：

I=imread('cameraman.tif');

imshow(I);

title('原始图像')

disp('原始图像大小：')

whos('I')

I=im2double(I);

%图像类型存储转换，将图像矩阵转换成双精度类型

T=dctmtx(8);

%离散余弦变换矩阵

B=blkproc(I,[8 8],'P1*x*P2',T,T');

mask=[1 1 1 1 0 0 0 0

1

1 1 0 0 0 0 0

1

1 0 0 0 0 0 0

1

0 0 0 0 0 0 0

0

0 0 0 0 0 0 0

0

0 0 0 0 0 0 0

0

0 0 0 0 0 0 0

0

0 0 0 0 0 0 0];

B2=blkproc(B,[8 8],'P1.*x',mask);

I2=blkproc(B2,[8 8],'P1*x*P2',T',T);

figure,imshow(I2);

title('压缩后的图像')

disp('压缩图像的大小：')

whos('I2')

运行结果：

(请自己运行查看)

原始图像大小：

Name Size Bytes Class Attributes

I 256x256 65536 uint8

压缩图像的大小：

Name Size Bytes Class Attributes

I2 256x256 524288 double

分析：由运行结果可以看出，经过DCT变换以后图像的大小几乎没有改变，我们知道DCT是一种空间变换，DCT变换的最大特点是对于一般的图像都能够将像块的能量集中于少数低频DCT系数上，这样就可能只编码和传输少数系数而不严重影响图像质量。DCT不能直接对图像产生压缩作用，但对图像的能量具有很好的集中效果，为压缩打下了基础。例如：一帧图像内容以不同的亮度和色度像素分布体现出来，而这些像素的分布依图像内容而变，毫无规律可言。但是通过离散余弦变换(DCT)，像素分布就有了规律。代表低频成份的量分布于左上角，而越高频率成份越向右下角分布。然后根据人眼视觉特性，去掉一些不影响图像基本内容的细节(高频分量)，从而达到压缩码率的目的。

利用小波变换实现图像压缩：

程序代码：

clear all

I=imread('cameraman.tif');

imshow(I); %显示图像

title('原始图像')

disp('原始图像I的大小:');

whos('I')

I=im2double(I);

[c,s]=wavedec2(I,2,'bior3.7'); %对图像用小波进行层分解

cal=appcoef2(c,s,'bior3.7',1); %提取小波分解结构中的一层的低频系数和高频系数

ch1=detcoef2('h',c,s,1); %提取二维水平方向细节系数

cv1=detcoef2('v',c,s,1); %提取二维垂直方向细节系数

cd1=detcoef2('d',c,s,1); %提取二维对角线方向细节系数

ca1=appcoef2(c,s,'bior3.7',1); %保留小波分解第一层低频信息

ca1=wcodemat(ca1,440,'mat',0); %首先对第一层信息进行量化编码

ca1=0.5*ca1; %改变图像高度

figure,image(ca1); %显示压缩后的图象

title('第一次压缩后的图像')

disp('第一次压缩图像的大小为：')

whos('ca1')

ca2=appcoef2(c,s,'bior3.7',2); %保留小波分解第二层低频信息进行压缩

ca2=wcodemat(ca2,440,'mat',0); %首先对第二层信息进行量化编码

ca2=0.25*ca2; %改变图像高度

figure,image(ca2); %显示压缩后的图象

title('第二次压缩后的图像')

disp('第二次压缩图像的大小为：')

whos('ca2')

运行结果：

(请自己运行查看)

原始图像I的大小:

Name Size Bytes Class

I 256x256 65536 uint8

array

Grand total is 65536 elements using 65536 bytes

第一次压缩图像的大小为：

Name Size Bytes Class

ca1 135x135 145800 double

array

Grand total is 18225 elements using 145800 bytes

第二次压缩图像的大小为：

Name Size Bytes Class

ca2 75x75 45000 double

array

Grand total is 5625 elements using 45000 bytes

分析：利用小波变换实现图像压缩，是利用小波分解图像的高频部分，只保留低频部分，从图中可以看出，第一次压缩是提取图像中小波分解的第一层低频信息，此时压缩效果较好，第二次压缩是提取第二层的低频部分，其压缩效果远不如第一次压