数字图像处理实验指导书2010.doc

数字图像处理实验指导书2010

———————————————————————————————— 作者：
———————————————————————————————— 日期：

2

个人收集整理 勿做商业用途

《数字图像处理》实验

指导书

成功 编

中南大学地球科学与信息物理学院GIS中心
2010年9月

20

目 录
实验一、数字图像读取 1
实验二、图像变换 2
实验三、图像增强 5
实验四、图像分割(选修） 9
实验五、图像压缩与编码 10
实验六、应用KL变换进行图像的特征提取（选修） 18
实验七、图像融合 19

实验一、数字图像读取
一、实验目的：
（1）熟悉MATLAB软件；
(2)用MATLAB软件读取各种图像
二、实验内容：
（1） 计算出一幅灰度图像的直方图
clear
close all
I=imread(’图像1.bmp’）;
imhist（I）
title（'实验一（1) 直方图')；
（2） 对灰度图像进行简单的灰度线性变换,
figure
subplot(2，2,1）
imshow（I);
title（'试验2—灰度线性变换'）；
subplot（2，2，2）
histeq（I);
(3） 看其直方图的对应变化和图像对比度的变化。
原图像 f（m,n） 的灰度范围 [a,b］ 线形变换为图像 g(m，n），灰度范围[a',b’］
公式: g(m，n）=a'+(b’-a')＊ f（m,n） /（b-a)
figure
subplot（2，2,1)
imshow(I)
J=imadjust（I,[0.3，0。7],［0,1］，1)；
title(' 实验一(3）用 g(m，n)=a’+（b’-a’)＊ f（m，n) /(b-a）进行变换 ’）；
subplot（2，2，2）
imshow(J）
subplot(2,2,3）
imshow(I)
J=imadjust(I,［0.5 0。8]，[0,1]，1);
subplot（2,2,4)
imshow(J)
（4） 图像二值化 （选取一个域值,将图像变为黑白图像）
figure
subplot（2，2，1)
imshow(I）
J=find（I〈150)；
I（J)=0;
J=find(I〉=150);
I（J)=255;
title（' 实验一(4）图像二值化 （ 域值为150 )’）；
subplot（2，2，2）
imshow（I)
clc;
I=imread（’图像2。jpg’);
bw=im2bw(I,0.5);%选取阈值为0。5
figure；
imshow(bw） %显示二值图像

实验二、图像变换
一、实验目的
1．了解图像变换的意义和手段;
2．熟悉傅里叶变换的孩本性质;
3．热练掌握FFT方法及应用;
4．通过实验了解二维频谱的分布特点;
5．通过本实验掌握利用MATLAB编程实现数字图像的傅立叶变换。

二、实验原理
1．应用傅立叶变换进行图像处理
傅里叶变换是线性系统分析的一个有力工具，它能够定量地分析诸如数字化系统、采样点、电子放大器、卷积滤波器、噪音和显示点等的作用。通过实验培养这项技能，将有助于解决大多数图像处理问题。对任何想在工作中有效应用数字图像处理技术的人来说，把时间用在学习和掌握博里叶变换上是很有必要的。
2．傅立叶（Fourier）变换的定义
对于二维信号,二维Fourier变换定义为：

　　二维离散傅立叶变换为：

图像的傅立叶变换与一维信号的傅立叶变换变换一样，有快速算法，具体参见参考书目,有关傅立叶变换的快速算法的程序不难找到。实际上,现在有实现傅立叶变换的芯片，可以实时实现傅立叶变换。
3．利用MATLAB软件实现数字图像傅立叶变换的程序：
I=imread（‘原图像名.gif'）； ％读入原图像文件
imshow（I); %显示原图像
fftI=fft2（I); %二维离散傅立叶变换
sfftI=fftshift（fftI）； %直流分量移到频谱中心
RR=real(sfftI)； ％取傅立叶变换的实部
II=imag(sfftI)； ％取傅立叶变换的虚部
A=sqrt（RR。^2+II.^2）；%计算频谱幅值
A=(A—min(min（A)））/(max（max(A)）-min（min（A）))＊225；
%归一化
figure; ％设定窗口
imshow（A); ％显示原图像的频谱

三、实验内容：
1．傅立叶变换
熟悉其概念和原理，实现对一幅灰度图像的快速傅立叶变换，并求其变换后的系数分布。
2．离散余弦变换
熟悉其概念和原理，实现对一幅灰度和彩色图像作的离散余弦变换，选择适当的DCT系数阈值对其进行DCT反变换。

％ 图像的FFT变换
clc；
I=imread（'图像1。bmp'）；
subplot(1,2，1)
imshow(I)；
title('原图');
subplot（1，2，2）
imhist(I)；
title(’直方图’）；
colorbar；
J=fft2（I);
figure；
subplot(1，2,1）
imshow(J）;
title('FFT变换结果’）;
subplot（1,2,2)
K=fftshift（J);
imshow（K);
title（'零点平移')；
figure;
imshow(log(abs(K））,［]），colormap(jet（64)），colorbar;
title（’系数分布图'）;
% 图像的DCT变换
RGB=imread（’图像1。bmp'）；
figure；
subplot(1，2,1）
imshow(RGB）；
title(’彩色原图'）；
a=rgb2gray(RGB);
subplot(1，2,2)
imshow(a)；
title（’灰度图’);
figure;
b=dct2（a）;
imshow(log（abs(b)）,［])，colormap（jet（64）），colorbar;
title('DCT变换结果’)；
figure;
b（abs（b）〈10)=0;
% idct
c=idct2（b）/255;
imshow(c)；
title（’IDCT变换结果’);
3．小波变换
实验内容： 熟悉小波变换的概念和原理,熟悉matlab小波工具箱主要函数的使用。利用二维小波分析对一幅图像作2层小波分解，并在此基础上提取各层的低频信息实现图像的压缩。
程序如下:
clc
close all
clear
a=imread('实验图像.tif’)；
subplot(1，2,1）；
imshow（a);
title（'原始图像')；
I=rgb2gray(a）；
subplot（1,2，2);
imshow（I);
title（'原始图像的灰度图');
％ 进行二维小波变换
[a，b］ = wavedec2(I， 2, 'bior3.7’);
% 提取各层低频信息
figure；
c = appcoef2( a, b， ’bior3。7’, 1 )；
subplot(1，2，1);
imshow（c， ［]）;
title('一层小波变换结果’）；
d = appcoef2（ a, b, 'bior3.7’, 2 )；
subplot(1,2,2);
imshow(d， [］）；
title（’二层小波变换结果’）；

四、实验步骤
1．打开计算机，安装和启动MATLAB程序；程序组中“work”文件夹中应有待处理的图像文件；
2．利用MatLab工具箱中的函数编制FFT频谱显示的函数；
3．a)调入、显示“实验一”获得的图像；图像存储格式应为“。gif";
b)对这三幅图像做FFT并利用自编的函数显示其频谱；
c)讨论不同的图像内容与FFT频谱之间的对应关系。
4．记录和整理实验报告.

五、实验仪器
1计算机；
2 MATLAB程序；
3移动式存储器(软盘、U盘等)。
4记录用的笔、纸.

六、实验报告内容
1叙述实验过程；
2提交实验的原始图像和结果图像。

七、思考题
1．傅里叶变换有哪些重要的性质?
2．图像的二维频谱在显示和处理时应注意什么？

实验三、图像增强
一、实验目的:
1掌握灰度直方图的概念及其计算方法；
2熟练掌握直方图均衡化和直方图规定化的计算过程;
3熟练掌握空域滤波中常用的平滑和锐化滤波器；
4掌握色彩直方图的概念和计算方法；
5利用MATLAB程序进行图像增强。
二、实验原理
图像增强是指按特定的需要突出一幅图像中的某些信息，同时，消弱或去除某些不需要的信息的处理方法。其主要目的是处理后的图像对某些特定的应用比原来的图像更加有效。图像增强技术主要有直方图修改处理、图像平滑化处理、图像尖锐化处理和彩色处理技术等。本实验以直方图均衡化增强图像对比度的方法为主要内容,其他方法同学们可以在课后自行联系.
1、直方图均衡化增强
直方图是多种空间域处理技术的基础.直方图操作能有效地用于图像增强.除了提供有用的图像统计资料外，直方图固有的信息在其他图像处理应用中也是非常有用的，如图像压缩与分割。直方图在软件中易于计算,也适用于商用硬件设备，因此，它们成为了实时图像处理的一个流行工具.
直方图是图像的最基本的统计特征，它反映的是图像的灰度值的分布情况。直方图均衡化的目的是使图像在整个灰度值动态变化范围内的分布均匀化,改善图像的亮度分布状态，增强图像的视觉效果。灰度直方图是图像预处理中涉及最广泛的基本概念之一。
图像的直方图事实上就是图像的亮度分布的概率密度函数，是一幅图像的所有象素集合的最基本的统计规律。直方图反映了图像的明暗分布规律，可以通过图像变换进行直方图调整，获得较好的视觉效果。
直方图均衡化是通过灰度变换将一幅图像转换为另一幅具有均衡直方图,即在每个灰度级上都具有相同的象素点数的过程。
下面给出直方图均衡化增强图像对比度的MATLAB程序:
I=imread（‘原图像名.gif）; ％ 读入原图像
J=histeq(I)； ％对原图像进行直方图均衡化处理
Imshow(I）； ％显示原图像
Title（‘原图像’)； ％给原图像加标题名
Figure；imshow(J）； ％对原图像进行屏幕控制；显示直方图均衡化后的图像
Title（‘直方图均衡化后的图像’） ； %给直方图均衡化后的图像加标题名
Figure； subplot(1，2,1） ；％对直方图均衡化后的图像进行屏幕控制；作一幅子图作为并排两幅图的第1幅图
Imhist(I,64); ％将原图像直方图显示为64级灰度
Title(‘原图像直方图’) ; %给原图像直方图加标题名
Subplot（1,2,2); ％作第2幅子图
Imhist（J,64) ； ％将均衡化后图像的直方图显示为64级灰度
Title(‘均衡变换后的直方图’） ; ％给均衡化后图像直方图加标题名
处理后的图像直方图分布更均匀了,图像在每个灰度级上都有像素点。从处理前后的图像可以看出，许多在原始图像中看不清楚的细节在直方图均衡化处理后所得到的图像中都变得十分清晰.
2、滤波增强（模板卷积增强）
（1）平滑：平滑的目的是模糊和消除噪声.平滑是用低通滤波器来完成,在空域中全是正值。
(2)锐化：锐化的目的是增强被模糊的细节。锐化是用高通滤波器来完成,在空域中，接近原点处为正，在远离原点处为负.
利用模板进行图像增强就是进行模板卷积。
1、 利用二个低通邻域平均模板（3×3和9×9)对一幅图像进行平滑，验证模板尺寸对图像的模糊效果的影响。
2、 利用一个低通模板对一幅有噪图像(GAUSS白噪声）进行滤波，检验两种滤波模板(分别使用一个5×5的线性邻域平均模板和一个非线性模板：3×5中值滤波器)对噪声的滤波效果.
3、 选择一个经过低通滤波器滤波的模糊图像,利用sobel和prewitt水平边缘增强高通滤波器（模板）对其进行高通滤波图像边缘增强，验证模板的滤波效果。
4、 选择一幅灰度图像分别利用一阶Sobel算子和二阶Laplacian算子对其进行边缘检测,验证检测效果。

三、实验内容
1、直方图均衡化增强
下面给出直方图均衡化增强图像对比度的MATLAB程序：
I=imread(‘图片1。gif'）; ％ 读入原图像
J=histeq（I); %对原图像进行直方图均衡化处理
Imshow(I)； %显示原图像
Title(‘原图像’）; ％给原图像加标题名
figure（0，0，400,400）； imshow（J）； ％对原图像进行屏幕控制;显示直方图均衡化后的图像
Title（‘直方图均衡化后的图像’) ； %给直方图均衡化后的图像加标题名
figure； subplot（1,2,1） ;%对直方图均衡化后的图像进行屏幕控制；作一幅子图作为并排两幅图的第1幅图
imhist(I,64）； ％将原图像直方图显示为64级灰度
Title（‘原图像直方图’) ; ％给原图像直方图加标题名
Subplot（1,2，2); %作第2幅子图
Imhist（J,64) ； %将均衡化后图像的直方图显示为64级灰度
Title（‘均衡变换后的直方图’） ; ％给均衡化后图像直方图加标题名
处理后的图像直方图分布更均匀了，图像在每个灰度级上都有像素点。从处理前后的图像可以看出，许多在原始图像中看不清楚的细节在直方图均衡化处理后所得到的图像中都变得十分清晰。

2、利用低通邻域平均模板进行平滑：
I=imread(’ Couple。bmp '）;
subplot（1,3，1)；
imshow(I);
title（’原图');
J=fspecial（'average’)；
J1=filter2（J，I)/255；
subplot(1,3,2）;
imshow(J1）;
title('3＊3滤波'）；
K=fspecial（'average’,9）；
K1=filter2（K，I)/255;
subplot（1,3，3)；
imshow（K1）；
title('9＊9滤波'）;
3、中值滤波和平均滤波
I=imread(’ Couple.bmp ’）;
J=imnoise(I，’gaussian’,0,0。01);
subplot（2,2，1)；
imshow（I)；
title('原图'）；
subplot(2,2,2）;
imshow（J）；
title（’noise’）；
K=fspecial（'average’，5);
K1=filter2(K，J）/255；
subplot（2,2，3)；
imshow（K1)；
title('average'）；
L=medfilt2（J，［3 5］);
subplot(2,2,4）;
imshow(L）;
title（’medium')；
4、高通滤波边缘增强
I=imread（' Couple.bmp ');
subplot(2，2，1);
imshow(I）；
title(’original pic');
J=fspecial（'average’，3);
J1=conv2（I，J)/255；
%J1=filter2(J,I)/255；
subplot(2，2，2);
imshow（J1)；
title('3＊3lowpass’）;
K=fspecial（’prewitt’);
K1=filter2(K,J1）＊5；
subplot(2,2，3）;
imshow(K1);
title（'prewitt');
L=fspecial（'sobel');
L1=filter2(L，J1）*5；
subplot（2，2，4)；
imshow（L1）；
title(’sibel'）;
5、边缘检测
分别用sobel和laplacian算子来进行，程序如下：
I=imread('Couple.bmp'）；
subplot(1,3,1）；
imshow(I）；
title（’original pic’);
K=fspecial（’laplacian',0.7）；
K1=filter2（K,I)/100；
subplot(1,3，2)；
imshow(K1）；
title（'laplacian'）；
L=fspecial(’sobel’);
L1=filter2（L,I）/200;
subplot（1，3，3）；
imshow(L1）；
title('sibel’）；

四、实验步骤
1打开计算机，启动MATLAB程序；程序组中“work”文件夹中应有待处理的图像文件;
2调入“实验一”中获取的数字图像，并进行计算机均衡化处理；
3显示原图像的直方图和经过均衡化处理过的图像直方图。
4记录和整理实验报告

五、实验仪器
1计算机；
2 MATLAB程序;
3移动式存储器（软盘、U盘等)。
4记录用的笔、纸.

六、实验报告内容
1叙述实验过程；
2提交实验的原始图像和结果图像。

七、思考题
1． 直方图是什么概念？它反映了图像的什么信息？
　　2． 直方图均衡化是什么意思？它的主要用途是什么？

实验四、图像分割（选修)
实验目的：（1) 学习边缘检测
（2） 学习灰度阀值分割
实验内容：
(1）分别用sobel、Laplacian—Gaussian方法对一幅灰度图像进行边缘提取，给出对比结果。
i=imread（’eight。tif’);
figure;
subplot(2,2，1）；
imshow（i）;
title（’原始图像’）；
subplot（2，2,3);
imshow（i)；
title（’原始图像'）;
i1=edge（i，'sobel’)；
subplot(2，2，2)；
imshow(i1)；
title('sober方法提取的边缘’）；
i2=edge(i,’log’)；
subplot(2，2,4);
imshow（i2）；
title（'Laplacian-Gaussian方法提取的边缘’）；
比较提取边缘的效果可以看出，sober算子是一种微分算子，对边缘的定位较精确，但是会漏去一些边缘细节.而Laplacian-Gaussian算子是一种二阶边缘检测方法,它通过寻找图像灰度值中二阶过零点来检测边缘并将边缘提取出来,边缘的细节比较丰富。通过比较可以看出Laplacian—Gaussian算子比sober算子边缘更完整，效果更好。
（2）利用双峰法对一幅灰度图像进行灰度分割处理
i=imread（’eight.tif'）;
subplot(1，2，1）;
imhist(i)；
title('原始图像直方图'）；
thread=130/255；
subplot（1，2,2）；
i3=im2bw(i,thread）；
imshow（i3);
title(’分割结果'）;
根据原图像的直方图，发现背景和目标的分割值大约在130左右，并将灰度图像转为二值图像，分割效果比较理想。
（实例：基于MATLAB 图像处理的汽车牌照识别系统）

实验五、图像压缩与编码
一、实验目的
1． 理解有损压缩和无损压缩的概念；
2． 理解图像压缩的主要原则和目的;　
3． 了解几种常用的图像压缩编码方式；
4． 利用MATLAB程序进行图像压缩.
二、实验原理
1。图像压缩原理
图像压缩主要目的是为了节省存储空间，增加传输速度。图像压缩的理想标准是信息丢失最少,压缩比例最大.不损失图像质量的压缩称为无损压缩，无损压缩不可能达到很高的压缩比；损失图像质量的压缩称为有损压缩，高的压缩比是以牺牲图像质量为代价的。压缩的实现方法是对图像重新进行编码，希望用更少的数据表示图像。
信息的冗余量有许多种，如空间冗余，时间冗余，结构冗余，知识冗余，视觉冗余等，数据压缩实质上是减少这些冗余量。高效编码的主要方法是尽可能去除图像中的冗余成分，从而以最小的码元包含最大的图像信息。
编码压缩方法有许多种，从不同的角度出发有不同的分类方法,从信息论角度出发可分为两大类。
（1）.冗余度压缩方法，也称无损压缩、信息保持编码或嫡编码。具体说就是解码图像和压缩编码前的图像严格相同，没有失真，从数学上讲是一种可逆运算。
（2）信息量压缩方法,也称有损压缩、失真度编码或烟压缩编码。也就是说解码图像和原始图像是有差别的，允许有一定的失真。
应用在多媒体中的图像压缩编码方法，从压缩编码算法原理上可以分为以下3类:
（1）无损压缩编码种类
哈夫曼(Huffman）编码，算术编码,行程（RLE）编码,Lempel zev编码。
（2）有损压缩编码种类
预测编码，DPCM，运动补偿;
频率域方法：正交变换编码(如DCT），子带编码；
空间域方法:统计分块编码；
模型方法：分形编码，模型基编码;
基于重要性:滤波，子采样，比特分配,向量量化；
（3）混合编码.
有JBIG,H261，JPEG,MPEG等技术标准。
本实验主要利用MATLAB程序进行离散余弦变换（DCT）压缩和行程编码（Run Length Encoding， RLE)。

1) 离散余弦变换（DCT）图像压缩原理
离散余弦变换DCT在图像压缩中具有广泛的应用，它是JPEG、MPEG等数据压缩标准的重要数学基础。
和相同图像质量的其他常用文件格式(如GIF(可交换的图像文件格式)，TIFF(标签图像文件格式)，PCX(图形文件格式))相比，JPEG是目前静态图像中压缩比最高的。JPEG比其他几种压缩比要高得多，而图像质量都差不多(JPEG处理的图像只有真彩图和灰度图）。正是由于其高压缩比，使得JPEG被广泛地应用于多媒体和网络程序中。JPEG有几种模式,其中最常用的是基于DCT变换的顺序型模式，又称为基本系统（Baseline）。
用DCT压缩图像的过程为：
(1)首先将输入图像分解为8×8或16×16的块，然后对每个子块进行二维DCT
变换。
（2)将变换后得到的量化的DCT系数进行编码和传送，形成压缩后的图像格
式.
用DCT解压的过程为：
（1)对每个8×8或16×16块进行二维DCT反变换。
（2)将反变换的矩阵的块合成一个单一的图像。
余弦变换具有把高度相关数据能量集中的趋势，DCT变换后矩阵的能量集中在矩阵的左上角,右下的大多数的DCT系数值非常接近于0.对于通常的图像来说，舍弃这些接近于0的DCT的系数值,并不会对重构图像的画面质量带来显著的下降。所以，利用DCT变换进行图像压缩可以节约大量的存储空间。压缩应该在最合理地近似原图像的情况下使用最少的系数。使用系数的多少也决定了压缩比的大小。
在压缩过程的第2步中，可以合理地舍弃一些系数,从而得到压缩的目的.在压缩过程的第2步，还可以采用RLE和Huffman编码来进一步压缩.

2）行程编码（RLE)原理：
　　例如如下这幅 的二值图像，

　　如果采用行程编码可以按如下格式保存

　　其中10和8表示图像的宽和高.在这个小例子中行程编码并没有起到压缩图像的作用。这是由于这个图的尺寸过小，当图像尺寸较大时行程编码还是不错的无损压缩方法。对于灰度图像和二值图像,用行程编码—般都有很高的压缩率。行程编码方法实现起来很容易,对于具有长重复值的串的压缩编码很有效，例如：对于有大面积的阴影或颜色相同的图像，使用这种方法压缩效果很好。很多位图文件格式都采用行程编码,如TIFF，PCX，GEM，BMP等。

3. 图像压缩编码的MATLAB程序语句
1） 利用DCT变换进行图像压缩的MATLAB程序
RGB = imread(’原图像名。tif')；
I = rgb2gray(RGB）；
J = dct2(I）;
imshow（log（abs（J）),[]), colormap（jet(64)）, colorbar
J(abs(J） 〈 10) = 0;
K = idct2(J);
figure,imshow(I）
figure，imshow（K,［0 255]）

2）利用离散余弦变换进行JPEG图像压缩
I=imread（‘灰度图。tif'); ％读入原图像；
I=im2double（I）; ％将原图像转为双精度数据类型;
T=dctmtx（8）; %产生二维DCT变换矩阵
B=blkproc(I,［8 8］，’P1＊x＊P2’，T,T’）; ％计算二维DCT，矩阵T及其转置T’是DCT函数P1*x*P2的参数
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］； ％二值掩膜，用来压缩DCT系数，只留下DCT系数中左上角的10个
B2=blkproc（B,［8 8］,’ P1。*x。',mask）； %只保留DCT变换的10个系数
I2= blkproc(B2，［8，8］,’P1*x*P2’，T’,T); %逆DCT，重构图像
Subplot(1，2,1）；
Imshow（I）;title(‘原图像'）； %显示原图像
Subplot(1，2,2)；
Imshow（I2）;title（‘压缩图像’);%显示压缩后的图像。对比原始图像和压缩后的图像，虽然舍弃了85%的DCT系数,但图像仍然清晰(当然有一些质量损失）

3）利用行程编码（RLE）进行图像压缩
I=checkerboard(10,2)； %调入原图像
[m n］=size（I）；
J=［]；
for i=1:m
value=I（i,1)；
num=1;
for j=2：n
if I（i,j)==value
num=num+1；
else
J=［J num value]；
num=1；
value=I(i,j);
end
end
I=[J num value 0 0]; ％添加的行判断位 0 0
end
disp(‘原图像大小：’)
whos(‘I'）;
disp（‘压缩图像大小：’）
whos(‘J'）；
disp（‘图像的压缩比：’)
disp（m*n/length（J)）

三、实验内容
（1）实现基本JPEG的压缩和编码分三个步骤:
1． 首先通过DCT变换去除数据冗余；
2． 使用量化表对DCT系数进行量化；
3． 对量化后的系数进行Huffman编码。
具体源程序由主程序及两个子程序（DCT量化、Huffman编码)组成：
1．主程序
I=imread（’autumn.tif’);
yiq=rgb2ntsc(I);
my=［16 11 10 16 24 40 51 61;12 12 14 19 26 58 60 55；
14 13 16 24 40 57 69 56;14 17 22 29 51 87 80 62；
18 22 37 56 68 109 103 77;24 35 55 64 81 104 113 92；
49 64 78 87 103 121 120 101；72 92 95 98 112 100 103 99］;
miq=［17 18 24 47 99 99 99 99;18 21 26 66 99 99 99 99;
24 26 56 99 99 99 99 99；47 66 99 99 99 99 99 99；
99 99 99 99 99 99 99 99；99 99 99 99 99 99 99 99;
99 99 99 99 99 99 99 99；99 99 99 99 99 99 99 99］；
I1=yiq(：,:,1);I2= yiq (:,:，2）;
[m n］=size（I1)；
t1=8；ti1=1；
while（t1〈m）
t1=t1+8;ti1=ti1+1；
end
t2=8;ti2=1;
while（t2〈n)
t2=t2+8；ti2=ti2+1;
end
times=0;
for k=0：ti1-2
for j=0:ti2—2
dct8x8(I1（k＊8+1:k*8+8,j＊8+1:j＊8+8),my,times*64+1)；
dct8x8（I2(k*8+1：k*8+8,j*8+1:j*8+8)，miq，times*64+1)；
times=times+1；
end
blkproc （I2(k＊8+1：k＊8+8，j*8+1:t2），［8 8]， T）；
end
for j=0:ti2-2
dct8x8(I1(k*8+1:t1,j＊8+1：j＊8+8）,times*64+1);
times=times+1；
end
dct8x8（I1(k＊8+1：t1，j＊8+1:t2）,times*64+1);
2．function dct8x8（I,m,s） ％定义DCT量化子程序
T=inline('dctmtx(8)')；
y=blkproc（I，［8 8]，T)；
y=round(y。/m);
p=1;te=1;
while(p〈=64)
for q=1:te
y1(s+p)=y（te-q+1，q）；p=p+1；
end
te=te+1;
for q=te:—1:1
y1(s+p)=y（te-q+1,q)；p=p+1;
end
te=te+1;
end
f=haffman（y1);
c（s:s+64，1)=f(：,1）;c(s：s+64,2）=f（：,2）;c（s：s+64,3)=f（：,3)
3．function c=haffman（I） ％定义Huffman编码子程序
[m，n］=size（I)；
p1=1;s=m＊n;
for k=1:m
for h=1：n
f=0；
for b=1:p1—1
if(c(b，1)==I(k，h）） f=1；break;end
end
if（f==0） c（p1，1)=I（k,h);p1=p1+1；end
end
end
for g=1:p1-1
p（g)=0;c（g,2)=0；
for k=1:m
for h=1：n
if（c（g,1）==I(k,h)） p(g)=p(g）+1；end
end
end
p(g)=p(g)/s;
end
pn=0;po=1；
while（1)
if(pn>=1。0) break；
else
[pm, p2］=min(p（1：p1—1))；p（p2)=1。1;
［pm2，p3]=min（p(1:p1-1))；p（p3）=1.1;
pn=pm+pm2;p（p1)=pn；
tree(po，1）=p2；tree(po，2）=p3;
po=po+1；p1=p1+1;
end
end
for k=1：po—1
tt=k；m1=1;
if(or（tree（k，1）〈g,tree（k，2)〈g）)
if（tree（k,1）〈g)
c（tree（k,1)，2）=c(tree（k,1),2）+m1;
m2=1;
while（tt〈 po-1)
m1=m1＊2;
for h=tt:po-1
if（tree（h，1)==tt+g）
c（tree（k，1),2）=c(tree(k,1），2）+m1；
m2=m2+1;tt=h；break;
else if(tree（h,2)==tt+g）
m2=m2+1；tt=h;break;
end
end
end
c（tree(k,1）,3)=m2；
end
tt=k;m1=1；
if(tree(k，2）〈g）
m2=1；
while（tt〈 po-1)
m1=m1*2；
for h=tt:po-1
if(tree（h,1)==tt+g)
c（tree（k，2），2)=c（tree(k，2)，2)+m1；
m2=m2+1；tt=h;break;
else if（tree(l，2）==tt+g)
m2=m2+1，tt=h;break;
end
end
end
c(tree（k,2)，3)=m2;
end
end
end
(2)JPEG2000采用小波变换编码,
小波变换压缩编码实现程序为
load wbarb;
subplot（2,2，1），image(X）;colormap(map)
title（'原始图像');
[c,s］=wavedec2（X,2, ’bior3。7’）；
thr=20;
ca1=appcoef2（c,s, 'bior3.7’,1）;
ch1=detcoef2（’h’,c,s，1）;
cv1=detcoef2（'v’，c，s，1);
cd1=detcoef2（’d'，c,s，1）;
a1=wrcoef2(’a'，c,s, ’bior3.7’,1);
h1=wrcoef2('h',c,s， 'bior3。7’，1)；
v1=wrcoef2('v'，c,s， 'bior3.7',1）；
d1=wrcoef2('d’,c,s， 'bior3.7'，1)；
c1=［a1,h1，v1，d1];
ca1=appcoef2(c，s, ’bior3。7’,1）；
ca1=wcodemat(ca1,440, ’mat'，0）；
ca1=0.5＊ca1
subplot(2，2,2),image(ca1)
title（