Matlab图像处理

原文：http://www.ilovematlab.cn/forum.php?mod=viewthread&s_tid=followedthread&tid=164340

基本概念

理论上图像是一种二维的连续函数，在计算机上对图像进行数字处理时，首先必须对其在空间和亮度上进行数字化，这就是图像的采样和量化的过程。空间坐标（x,y）的数字化称为图像采样，而幅值数字化称为灰度级量化。
对一幅图像采样时，若每行（横向）采样数为M，每列（纵向）采样数为N，则图像大小为M*N个像素， f(x,y)表示点(x,y) 处的灰度值，则F(x,y)构成一个M*N 实数矩阵
把采样后所得的各像素灰度值从模拟量到离散量的转换称为图像灰度的量化。量化是对图像幅度坐标的离散化，它决定了图像的幅度分辨率。

量化的方法包括：分层量化、均匀量化和非均匀量化。分层量化是把每一个离散样本的连续灰度值只分成有限多的层次；均匀量化是把原图像灰度层次从最暗至最亮均匀分为有限个层次，如果采用不均匀分层就称为非均匀量化。
当图像的采样点数一定时，采用不同量化级数的图像质量不一样。量化级数越多，图像质量越好；量化级数越少，图像质量越差。量化级数小的极端情况就是二值图像。
在C语言中，对M×N数字图像处理的核心代码如下：
for (j=1;j<N+1;j++)
   for(i=1;i<M+1;i++)
   ｛ 对I(i,j)的具体运算
｝；
在Matlab中，对M×N数字图像处理的核心代码如下：
for i=1:N
    for j=1:M
      对I(i,j)的具体运算
    end
end

图像变换的方法包括空间域变换、频率域变换、时频域变换、基于经典数学理论的变换、基于现代数学理论的变换。
以下程序段是用于图像放缩（空间变换）的MATLAB源程序：
*********************************************************************
function newImage=resample1(image,newRow,newCol)
% 功能：对图像进行缩放
% 输入：image-需要进行缩放的灰度图像；
% newRow-缩放后新图像的行数；
% newCol-缩放后新图像的列数；
% 输出：newImage-缩放后的图像。

[row,col]=size(image);
image_larger=zeros((row+1),(col+1));
image_larger(1:row,1:col)=image;
image=image_larger;
newImage=zeros(newRow,newCol);
for i=0:(newRow-1)
for j=0:(newCol-1)
x=j*(col-1)/(newCol-1);
y=i*(row-1)/(newRow-1);
fx=floor(x);
fy=floor(y);

area_ul=(x-fx)*(y-fy);
area_ur=(fx+1-x)*(y-fy);
area_bl=(x-fx)*(fy+1-y);
area_br=(fx+1-x)*(fy+1-y);

newImage(i+1,j+1)=...
image(fy+1,fx+1)*area_br+...
image(fy+1,fx+2)*area_bl+...
image(fy+2,fx+1)*area_ur+...
image(fy+2,fx+2)*area_ul;
end
end
*********************************************************************
在MATLAB中，图像的缩放也可以调用imresize函数来实现。imresize函数的调用格式如下：
B = imresize(A,m,method)
imrersize函数使用由参数method指定的插值运算来改变图像的大小。method的几种可选值： 'nearest'(默认值)最近邻插值； 'bilinear'双线性插值； 'bicubic'双三次插值； B = imresize(A,m)表示把图像A放大m倍。

（1）图像文件的读取
利用imread函数可以完成图像文件的读取操作。常用语法格式为：
I=imread(‘filename’,‘fmt’)或I=imread(‘filename.fmt’)；
其作用是将文件名用字符串filename表示的、扩展名用字符串fmt（表示图像文件格式）表示的图像文件中的数据读到矩阵I中。当filename中不包含任何路径信息时，imread会从当前工作目录中寻找并读取文件。要想读取指定路径中的图像，最简单的方法就是在filename中输入完整的或相对的地址。 MATLAB支持多种图像文件格式的读、写和显示。因此参数fmt常用的可能值有：
‘bmp’ Windows位图格式
‘jpg’or‘jpeg’ 联合图像专家组格式
‘tif’or‘tiff’ 标志图像文件格式
‘gif’ 图形交换格式
‘pcx’ Windows画刷格式
‘png’ 可移动网络图形格式
‘xwd’ X Window Dump格式
例如，命令行
>> I=imread(‘lena.jpg’);
将JPEG图像lena读入图像矩阵I中。

（2） 图像文件的写入（保存）
利用imwrite完成图像的输出和保存操作，也完全支持也完全支持上述各种
图像文件的格式。其语法格式为：
imwrite(I,‘filename’,‘fmt’)或imwrite(I,‘filename.fmt’)；
其中的I、filename和fmt的意义同上所述。
注意事项：当利用imwrite函数保存图像时，MATLAB默认的保存方式是将其简化为uint8的数据类型。与读取文件类型类似，MATLAB在文件保存时还支持16位的PNG和TIFF图像。所以，当用户保存这类文件时，MATLAB就将其存储在uint16中。

（3）图像文件的显示
图像的现实过程是将数字图像从一组离散数据还原为一幅可见图像的过程。
MATLAB的的图像处理工具箱提供了多种图像显示技术。例如imshow可以直接从文件显示多种图像；image函数可以将矩阵作为图像 ；colorbar函数可以用来显示颜色条；montage函数可以动态显示图像序列。这里仅对常用的显示函数进行介绍。
①图像的显示
imshow函数是最常用的显示各种图像的函数，其调用格式如下：
imshow(I,N)；
imshow(I,N)用于显示灰度图像，其中I为灰度图像的数据矩阵，N为灰度级数目，默认值为256。
例如下面的语句用于显示一幅灰度图像：
>> I=imread(‘lena.jpg’)；
>> imshow(I);
如果不希望在显示图像之前装载图像，那么可以使用以下格式直接进行图像文件的显示：
imshow filename
其中，filename为要显示的图像文件的文件名。
注意事项：该文件名必须带有合法的扩展名（指明文件格式），且该图像文件必须保存在当前目录下，或在MATLAB默认的目录下。
②添加色带
colorbar函数可以给一个坐标轴对象添加一条色带。如果该坐标轴对象包含一个图像对象，则添加的色带将指示出该图像中不同颜色的数据值。这对于了解被显示图像的灰度级特别有用。其调用格式为：
colorbar
③显示多幅图像
显示多幅图像最简单的方法就是在不同的图形窗口中显示它们。imshow总是在当前窗口中显示一幅图像，如果用户想连续显示两幅图像，那么第二幅图像就会替代第一幅图像。为了避免图像在当前窗口中的覆盖现象，在调用imshow函数显示下一幅图像之前可以使用figure命令来创建一个新的窗口。例如：
imshow(I1);
figure, imshow(I2);
figure, imshow(I3);
有时为了便于在多幅图像之间进行比较，需要将这些图像显示在一个图形窗口中。达到这一目的有两种方法：一种方法是联合使用imshow和subplot函数，但此方法在一个图形窗口只能有一个调色板；另一种方法是联合使用subimage和subplot函数，此方法可在一个图形窗口内使用多个调色板。
subplot函数将一个图形窗口划分为多个显示区域，其调用格式如下：
subplot(m,n,p)
subplot函数将图形窗口划分为m（行）×n（列）个显示区域，并选择第p个区域作为当前绘图区。

（4） 图像文件信息的查询
imfinfo函数用于查询图像文件的有关信息，详细地显示出图像文件的各种属性。其语法格式为：
info=imfinfo(‘filename’,‘fmt’)或info=imfinfo(‘filename.fmt’)
或imfinfo filename.fmt
imfinfo函数获取的图像文件信息依赖于文件类型的不同而不同，但至少应包
含以下内容：
文件名。如果该文件不在当前目录下，还包含该文件的完整路径。
文件格式。
文件格式的版本号。
文件最后一次修改的时间。
文件的大小。以字节为单位。
图像的宽度。
图像的高度。
每个像素所用的比特数。也叫像素深度。
图像类型。即该图像是真彩色图像、索引图像还是灰度图像。

数字图像处理的方法有两类：空间域处理法和频域法。频域法首先是要将图像从空间域变换到频率域，然后在频率域对图像进行各种处理，再将所得到的结果进行反变换，从而达到图像处理的目的，频域法具有很多明显的优点。目前，频域变换被广泛运用于图像特征提取、图像增强、图像复原、图像数据压缩和图像识别等领域。
（1）离散傅立叶变换（DFT）
离散傅立叶变换（DFT）在数字信号处理和数字图像处理中应用十分广泛。使用离散傅立叶变换的根本原因有二：一是DFT的输入、输出均为离散形式的，这使得计算机非常容易操作；二是因为计算DFT存在快速算法，即快速傅立叶变换（FFT），因而计算比较方便。
空间域是由f(x,y)所组成的坐标系，其中x和y用作（空间）变量。频率系统是由F(u,v)所组成的坐标系，其中u和v用作（频率）变量。由u=0,1,2,…M-1和v=0,1,2,…N-1定义的M×N矩形区域常称为频率矩形。显然，频率矩形的大小与输入图像的大小相同。为了与MATLAB的实现形式相一致，将1/MN项放置于逆变换公式的前面。由于MATLAB中的数组索引是以1而不是以0开头的，所以MATLAB中的F(1,1)和f(1,1)相应于正变换和逆变换中的数字量F(0,0)和f(0,0)。在频域原点处变换值[如F(0,0)]称为傅立叶变换的直流（dc）分量，F(0,0)等于f(x,y)的平均值的MN倍。
MATLAB图像处理工具箱提供了一些函数来进行傅立叶变换：
① 函数：fft2——用于计算二维快速傅立叶变换。调用格式为：
Y = fft2(X)；
Y = fft2(X,M,N)；
其中X是输入图像矩阵，Y是X进行二维傅立叶变换后的图像矩阵；X和Y大小相同。在Y=fft2(X,M,N)中，按照M、N指定的值对图像进行剪切或补0后进行傅立叶变换,返回变换矩阵的大小为M×N。变换结果的左上、右上、左下、右下四个角的周围对应于低频成分，中央部位对应于高频成分。
② 函数：fftn——用于计算n维傅立叶变换。调用格式：
Y = fftn(X)；
Y = fftn(X,SIZE)；
其中Y=fftn(X)计算图像的n维傅立叶变换，输出图像Y与X大小相同。在Y=fftn(X,SIZE)函数中，按照SIZE指定的值对图像X进行剪切或补0后进行傅立叶变换，返回变换矩阵的大小为SIZE。
③ 函数：fftshift——将变换后的图像频谱中心从矩阵的原点移到矩阵的中心。其调用格式为：
Y = fftshift(X)；
Y = fftshift(X,DIM)；
其中fftshift用于调整 fft、fft2和fftn的输出结果。对于向量X，将其左右两半交换位置，对于矩阵X，将其一、三象限和二、四象限进行互换，对于高维向量X，将矩阵各维的两半进行互换。利用这个函数可使变换结果的零频率分量位于中心。但应注意，在进行反变换时，必须使用四角代表低频成分，中央对应高频部分的变换结果。
④ 函数：ifft2——用于计算图像的二维傅立叶反变换。其调用格式为：
Y = ifft2(X)；
Y = ifft2(X,M,N)；
其中ifft2用于返回图像的二维傅立叶反变换矩阵，其参数定义同fft2。
⑤ 函数：ifftn——用于计算图像的n维傅立叶反变换。调用格式为：
Y = ifftn(X)；
Y = ifftn(X,SIZE)；
其参数定义同fftn。
实例：对某矩阵进行零填充后，进行快速傅立叶变换并显示其频谱。如图1.4.1所示。
f=zeros(30,30);
f(5:24,13:17)=1;
subplot(1,2,1),imshow(f);
F=fft2(f,256,256);
F2=fftshift(log(abs(F))); % 计算对数幅值并使零频率系数位于图形的中心
subplot(1,2,2),imshow(F2,[-1 5]);
colorbar;
（2）离散余弦变换（DCT）
离散余弦变换（DCT，Discrete Cosine Transform）的变换核为实数的余弦函数，因而DCT的计算速度要比变换核为复指数的DFT要快得多。离散余弦变换是仅次于K-L变换的次最佳正交变换，且有这样的性质：许多有关图像的重要可视信息都集中在DCT变换的一小部分系数中，因此已被广泛应用到图像压缩编码、语音信号处理等众多领域，并成为许多图像编码国际标准的核心。
MATLAB图像处理工具箱实现离散余弦变换有两种方法：其一是使用函数dct2，该函数用一个基于FFT的算法来提高当输入较大的输入方阵时的计算速度。其二是使用由dctmtx函数返回的DCT变换矩阵，这种方法较适合于较小的输入方阵（例如8×8或16×16）
提供DCT函数分别为：
① 函数：dct2——实现图像的二维离散余弦变换。调用格式为：
B = dct2(A)；
B = dct2(A,[M N])；
B = dct2(A,M,N)；
其中A表示要变换的图像，M和N是可选参数，表示填充后的图像矩阵大小，如果m和n比图像A小，则进行变换之前，将图像A进行剪切，B表示变换后得到的图像矩阵，各元素为离散余弦变换的系数B(k1,k2)。
② 函数：dctmtx——主要用于实现较小输入矩阵的离散余弦变换，调用格式为：D = dctmtx(N)，其中D是返回N×N的DCT变换矩阵，如果矩阵A是N×N方阵，则A的DCT变换可用D×A×D’来计算。这在有时比dct2计算快，特别是计算大量小的相同尺寸DCT时，矩阵D只需计算一次，因而速度快。
例如，在实现JPEG压缩时，要多次实现大小为8×8的图像块的DCT，为了实现这种变换，首先利用语句D=dctmtx(8)，然后，对每一个图像块执行运算B=D×A×D’。这种实现方法比调用函数dct2要快。
③ 函数：idct2——实现图像的二维离散余弦反变换。调用格式为：
B = idct2(A)；
B = idct2(A,[M N])；
B = idct2(A,M,N)；
其中参数同dct2。

灰度直方图
在数字图像处理中，灰度直方图是最简单且最有用的工具，可以说，对图像的分析与观察直到形成一个有效的处理方法，都离不开直方图。直方图表达的信息是每种亮度的像素点的个数。直方图是图像的一个重要特征，因为直方图用少量的数据表达图像的灰度统计特征。
根据图像直方图的定义编写的求灰度图像Matlab源程序。
*********************************************************************
%读入图像；
I=imread('taishan.jpg');
%将RGB图像转换为灰度图像；
B0=rgb2gray(I);
%将图像矩阵的类型转换成双精度型，便于后续的运算；
B=double(B0);
%求图像的行数与列数；
s=size(B);
%建立一个数组，用于存储1～256灰度级出现的个数；
h=zeros(1,256);
%根据定义，计算各像素灰度值出现的个数；
for i=1:s(1)
for j=1:s(2)
k=B(i,j);
k=floor(k);
h(k+1)=h(k+1)+1;
end
end
% 显示图像；
subplot(121),imshow(B0);
subplot(122),plot(h)
*********************************************************************

非线性灰度值变换
这种方法的目标与增强对比度相反。当原图的动态范围太大，超出了某些显示设备所允许的动态范围时，可采用对数形式的变换函数进行动态范围压缩：I=imread('yellowriver.jpg');
X1=rgb2gray(I);
>> figure,imshow(X1);
c=255/log(256);
x=0:1:255;
y=c*log(1+x);
figure,plot(x,y)
xlabel('f'),ylabel('g')
title('intensitytransformation')
%绘制变换曲线
[m,n]=size(X1);
X2=double(X1);
for i=1:m
for j=1:n
g(i,j)=0;
g(i,j)=c*log(X2(i,j)+1);
end
end
figure,imshow(mat2gray(g))

链码（又称为freeman码）是用曲线起始点的坐标和边界点方向代码来描述曲线或边界的方法，常被用来在图像处理、计算机图形学、模式识别等领域中表示曲线和区域边界。常用的链码按照中心像素点邻接方向个数的不同，分为4连通链码和8连通链码。4连通链码的邻接点有4个，分别在中心点的上、下、左和右。8连通链码比4连通链码增加了4个斜方向，因为任意一个像素周围均有8个邻接点，而8连通链码正好与像素点的实际情况相符，能够准确地描述中心像素点与其邻接点的信息。因此，8连通链码的使用相对较多。

8连通边界的链码生成程序：

function out=chaincode8(image)

%功能：实现8连通链码

%输入： 二值图像

%输出：链码的结果



n=[0 1;-1 1;-1 0;-1 -1;0 -1;1 -1;1 0;1 1];

%设置标志

flag=1;

%初始输出的链码串为空

cc=[];

%找到起始点

[x y]=find(image==1);

x=min(x);

imx=image(x,:);

y=min(find(imx==1));

first=[x y];



dir=7;

while flag==1

          tt=zeros(1,8);

          newdir=mod(dir+7-mod(dir,2),8);

          for i=0:7

              j=mod(newdir+i,8)+1;

              tt(i+1)=image(x+n(j,1),y+n(j,2));

          end

    d=min(find(tt==1));

          dir=mod(newdir+d-1,8);

          %找到下一个像素点的方向码后补充在链码的后面

    cc=[cc,dir];

    x=x+n(dir+1,1);y=y+n(dir+1,2);

    %判别链码的结束标志

    if x==first(1)&&y==first(2)

        flag=0;

    end

end

out=cc;

在Matlab图像处理工具箱中，提供了专门的bwlabel( )函数，对二值图像的连接部分进行标记。 其调用格式如下：
L = bwlabel(BW,n)

该函数返回一个和输入的二值图像BW大小相同的L矩阵，包含了标记了BW中每个连通区域的类别标签，这些标签的值为1、2、num(连通区域的个数)。n的值为4或8，表示是按4连通寻找区域还是8连通寻找，如果参数省略，则默认为8。

Harris角点检测的MATLAB 实现代码。
********************************************************************

function [posr, posc]=Harris1(in_image,a)

% 输入：in_image－待检测的RGB图像数组；

% 输入：a－角点参数响应，取值范围为：0.04～0.06；

% 输出：posr 所检测出角点的行坐标向量；

% 输出：posc所检测出角点的行坐标向量；

in_image=rgb2gray(in_image);        % 将RGB图像转化成灰度图像

ori_im=double(in_image);            % unit8型转化为双精度double64型

%%%=================计算图像在x、y 两个方向的梯度==============

fx = [-1 0 1];                       % x方向梯度算子模板

Ix = filter2(fx,ori_im);               % x方向滤波

fy = [-1;0;1];                       % y方向梯度算子

Iy = filter2(fy,ori_im);                % y方向滤波

%%%=================计算两个方向的梯度乘积====================

Ix2 = Ix.^2;

Iy2 = Iy.^2;

Ixy = Ix.*Iy;

%%%=================使用高斯函数对梯度乘积进行加权=============

h= fspecial('gaussian',[7 7],2);        % 产生7*7的高斯窗函数，sigma=2

Ix2 = filter2(h,Ix2);

Iy2 = filter2(h,Iy2);

Ixy = filter2(h,Ixy);

%%%=================计算每个像元的Harris响应值=================

[height,width]=size(ori_im);

R = zeros(height,width);                           

for i = 1:height

    for j = 1:width

        M = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)];            

        R(i,j) = det(M)-a*(trace(M))^2;          % 像素(i,j)处的Harris响应值

    end

end

%%%=================去掉小于阈值的Harris响应值=================

Rmax=max(max(R));

t=0.01* Rmax;                               %  阈值

for i = 1:height

for j = 1:width

if R(i,j)<t

  R(i,j) = 0;

  end

end

end

%%%=================进行3×3邻域非极大值抑制=================      

corner_peaks=imregionalmax(R);          % 进行非极大抑制，窗口大小3*3

countnum=sum(sum(corner_peaks));

%%%=================显示所提取的Harris角点=================

[posr, posc] = find(corner_peaks== 1);       % posr是用于存放行坐标的向

% pos是用于存放列坐标的向量；

figure

imshow(in_image)

hold on

for i = 1 : length(posr)

    plot(posc(i),posr(i),'r+');

end

********************************************************************

空域滤波
空域滤波是指借助模板对图像进行邻域操作，输出图像每一个像素的取值都是根据模板对输入像素相应邻域内的像素进行计算得到的。根据其特点，空域滤波一般可以分为线性的和非线性的两类处理方法；按照其功能，可以分为平滑滤波器和锐化滤波器。平滑滤波器主要实现模糊和去除噪声。模糊主要是在提取较大目标前，去除太小的细节或将目标内的小间断连接起来。锐化滤波器是为了增强被模糊的细节边缘。
例程1是实现邻域平均法的Matlab程序。
例程1
*****************************************************
x = imread('smile.jpg');
[width,height,dim]=size(x);
%均值滤波 模块12*12
k=12;
t=0;
%图像像素类型转换
x = double(x);
y=x;
%用邻域平均法进行滤波
for i=round(k/2):1:(width-round(k/2)+1)
for j=round(k/2):1:(height-round(k/2)+1)
for m=i-round(k/2)+1:1:i+round(k/2)-1
for n=j-round(k/2)+1:1:j+round(k/2)-1
s=x(m,n);
t=t+s;
end
end
y(i,j)=round(t/(k*k));
t=0;
end
end
subplot(1,2,1),imshow(uint8(x));
subplot(1,2,2),imshow(uint8(y));

在Matlab中，可以通过如下指令实现邻域均值滤波：
H=fspecial('average',n)；
Y=filter2(B,H)；
其中，n为邻域的大小，B为待滤波的二维灰度图像矩阵，H为所构造的邻域均值滤波器。
*********************************************************************
% 读入图像并进行数据转换
I=imread('junli.jpg');
I=rgb2gray(I);
% 添加椒盐噪声
J = imnoise(I,'salt & pepper',0.02);
% 采用不同大小邻域的均值滤波器进行滤波
K1=filter2(fspecial('average',3),J);
K2=filter2(fspecial('average',5),J);
K3=filter2(fspecial('average',7),J);
% 显示结果
subplot(221), imshow(J)
subplot(222),imshow(uint8(K1))
subplot(223),imshow(uint8(K2))
subplot(224),imshow(uint8(K3))
*********************************************************************

频域线性滤波和空间滤波一样，其基础都是卷积定理。
例程3对人为添加了少量椒盐噪声的图像进行巴特沃思低通滤波。
********************************************************************
I=imread('tejing.jpg');
I=rgb2gray(I);
J1=imnoise(I,'salt & pepper'); % 叠加椒盐噪声
subplot(1,2,1)
imshow(J1),title('原始图像')
f=double(J1); % 数据类型转换，MATLAB不支持图像的无符号整型的计算
g=fft2(f); % 傅立叶变换
g=fftshift(g); % 转换数据矩阵
[M,N]=size(g);
nn=2; % 二阶巴特沃斯(Butterworth)低通滤波器
d0=25;
m=fix(M/2); n=fix(N/2);
for i=1:M
for j=1:N
d=sqrt((i-m)^2+(j-n)^2);
h=1/(1+0.414*(d/d0)^(2*nn)); % 计算低通滤波器传递函数
result(i,j)=h*g(i,j);
end
end
result=ifftshift(result);
J2=ifft2(result);
J3=uint8(real(J2));
subplot(1,2,2) % 显示滤波处理后的图像
imshow(J3),title('滤波处理后的图像')
********************************************************************

分水岭分割方法，是一种基于拓扑理论的数学形态学的分割方法，其基本思想是把图像看作是测地学上的拓扑地貌，图像中每一点像素的灰度值表示该点的海拔高度，每一个局部极小值及其影响区域称为集水盆，而集水盆的边界则形成分水岭。分水岭的概念和形成可以通过模拟浸入过程来说明。在每一个局部极小值表面，刺穿一个小孔，然后把整个模型慢慢浸入水中，随着浸入的加深，每一个局部极小值的影响域慢慢向外扩展，在两个集水盆汇合处构筑大坝，即形成分水岭。

分水岭的计算过程是一个迭代标注过程。分水岭比较经典的计算方法是L. Vincent提出的。在该算法中，分水岭计算分两个步骤，一个是排序过程，一个是淹没过程。首先对每个像素的灰度级进行从低到高排序，然后在从低到高实现淹没过程中，对每一个局部极小值在h阶高度的影响域采用先进先出(FIFO)结构进行判断及标注。

分水岭变换得到的是输入图像的集水盆图像，集水盆之间的边界点，即为分水岭。显然，分水岭表示的是输入图像极大值点。因此，为得到图像的边缘信息，通常把梯度图像作为输入图像，即

g(x,y)=grad(f(x,y))={[f(x,y)-f(x-1,y)]2[f(x,y)-f(x,y-1)]2}0.5

式中，f(x,y)表示原始图像，grad{.}表示梯度运算。

分水岭算法对微弱边缘具有良好的响应，图像中的噪声、物体表面细微的灰度变化，都会产生过度分割的现象。但同时应当看出，分水岭算法对微弱边缘具有良好的响应，是得到封闭连续边缘的保证的。另外，分水岭算法所得到的封闭的集水盆，为分析图像的区域特征提供了可能。

为消除分水岭算法产生的过度分割，通常可以采用两种处理方法，一是利用先验知识去除无关边缘信息。二是修改梯度函数使得集水盆只响应想要探测的目标。

为降低分水岭算法产生的过度分割，通常要对梯度函数进行修改，一个简单的方法是对梯度图像进行阈值处理，以消除灰度的微小变化产生的过度分割。即

g(x,y)=max(grad(f(x,y)),gθ)

式中，gθ表示阈值。

程序可采用方法：用阈值限制梯度图像以达到消除灰度值的微小变化产生的过度分割，获得适量的区域，再对这些区域的边缘点的灰度级进行从低到高排序，然后在从低到高实现淹没的过程，梯度图像用Sobel算子计算获得。对梯度图像进行阈值处理时，选取合适的阈值对最终分割的图像有很大影响，因此阈值的选取是图像分割效果好坏的一个关键。缺点：实际图像中可能含有微弱的边缘，灰度变化的数值差别不是特别明显，选取阈值过大可能会消去这些微弱边缘。

分水岭变换是从局部极小点开始，即只能是在梯度图中用， 原始图是转换后才能用于分水岭变换的

一般图像中存在多个极小值点，通常会存在过分割现象，可以采用梯度阈值分割改进或者采用标记分水岭算法将多个极小值区域连在一起

插值法是根据两个自变量的已知函数值求这两个自变量之间各自变量对应函数值的近似计算方法。

% 本程序能够实现分水岭算法


% 分水岭算法
clear, close all;
clc;

PathName='t ';%t为自填内容,下面p类似
FileName=[PathName 'p'];
Image=imread(FileName);

subplot(2,2,1);subimage(Image);title('原图')
I=imread(t’tank.jpg’);
I=rgb2gray(I);
Image=I;
Image=double(Image);
B=[1,1,1;1,1,1;1,1,1];%方形结构元
E8=[-1,0;-1,1;0,1;1,1;1,0;1,-1;0,-1;-1,-1]; % 8－连通结构元坐标
maskLenth=length(E8); % 结构元点的个数

[X,Y]=size(Image);

%原始图像image 赋值给A1
n=1;
A(:,:,n)=Image;
M=zeros(X,Y);
Mark_Image=zeros(X,Y);

%产生距离图
while sum(sum(A(:,:,n)))~=0
A(:,:,n+1)= imerode(A(:,:,n),B);
U(:,:,n)= (A(:,:,n)-A(:,:,n+1))*n;
U(:,:,n)=double(U(:,:,n));
M=M+U(:,:,n);
n=n+1;
end
n=n-1;
subplot(2,2,2);imagesc(M,[0,n]);title('距离图');

% 搜寻局部最大值，将其放入Deal_Image
Deal_Image=zeros(X,Y);
while n>0
for high=1:X
for width=1:Y
%********************************************************************
Mark_Bool=0;
if M(high,width)==n
%______________________________________________________________
for dot=1:maskLenth
i=E8(dot,1); j=E8(dot,2);
if high+i>=1 & width+j>=1 & high+i<=X & width+j<=Y & M(high+i,width+j)>M(high,width);
Mark_Bool=1;break;
end % if_end
end % for dot_end
%______________________________________________________________
if Mark_Bool==0;
Deal_Image(high,width)=M(high,width);
end %if end
%______________________________________________________________
end %if end
%********************************************************************
end %for-end
end %for-end
n=n-1;
end % while n=0 end
Deal_Image =[Deal_Image>=1]
subplot(2,2,3);subimage(Deal_Image);title('输出图像');

Mark_Number=1;

while n>0
for high=1:X
for width=1:Y
Mark_Bool=0;
%********************************************************************
if M(high,width)==n
%______________________________________________________________
for dot=1:maskLenth
i=E8(dot,1); j=E8(dot,2);
if high+i>=1 & width+j>=1 & high+i<=X & width+j<=Y & Mark_Image(high+i,width+j)>0;
Mark_Image(high,width)=Mark_Image(high+i,width+j);
Mark_Bool=1;break;
end % if_end
end % for dot_end
%______________________________________________________________
if Mark_Bool==0;
Mark_Image(high,width)=Mark_Number;
Mark_Number=Mark_Number+1;
end %if end
%______________________________________________________________
pause;
subplot(2,2,2);imagesc(Mark_Image,[0,Mark_Number]);title('输出图像');
end %if end
%********************************************************************
end %for-end
end %for-end
n=n-1;
end % while n=0 end


subplot(2,2,3);imagesc(Mark_Image,[0,Mark_Number]);title('分割后的图像');

uicontrol('Style','edit','string',['分割出区域：',Num2str(Mark_Number-1),'个'],...
'Position', [400 0 150 18],'FontSize',12,'FontWeight','light');

基于Hough变换检测圆的基本思想是：对参数空间适当量化，得到一个三维的累加器阵列，并计算图像每点强度的梯度信息得到边缘，再计算与边缘上的每一个像素 距离半径为 的所有点，同时将相应立方小格的累加器加1。当检测完毕后，对三维阵列的所有累加器求峰值，大于所设定阈值的峰值小格的坐标就对应着图像空间圆形边界的圆心。

********************************************************************

function [hough_space,hough_circle,para] = hough_circle(BW,step_r,step_angle,r_min,r_max,p)

% 功能：基于Hough变换检测圆

% 输入：

% BW-二值图像；

% step_r-检测的圆半径步长

% step_angle-角度步长，单位为弧度

% r_min-最小圆半径

% r_max-最大圆半径

% p-以p*hough_space的最大值为阈值，p取0，1之间的数

% 输出：

% hough_space-参数空间，h(a,b,r)表示圆心在(a,b)半径为r的圆上的点数

% hough_circl-二值图像，检测到的圆

% para-检测到的圆的圆心、半径



[m,n] = size(BW);

size_r = round((r_max-r_min)/step_r)+1;

size_angle = round(2*pi/step_angle);



hough_space = zeros(m,n,size_r);



[rows,cols] = find(BW);

ecount = size(rows);



% Hough变换

% 将图像空间(x,y)对应到参数空间(a,b,r)

% a = x-r*cos(angle)

% b = y-r*sin(angle)

for i=1:ecount

    for r=1:size_r

        for k=1:size_angle

            a = round(rows(i)-(r_min+(r-1)*step_r)*cos(k*step_angle));

            b = round(cols(i)-(r_min+(r-1)*step_r)*sin(k*step_angle));

            if(a>0&a<=m&b>0&b<=n)

                hough_space(a,b,r) = hough_space(a,b,r)+1;

            end

        end

    end

end



% 搜索超过阈值的聚集点

max_para = max(max(max(hough_space)));

index = find(hough_space>=max_para*p);

length = size(index);

hough_circle=zeros(m,n);

for i=1:ecount

    for k=1:length

        par3 = floor(index(k)/(m*n))+1;

        par2 = floor((index(k)-(par3-1)*(m*n))/m)+1;

        par1 = index(k)-(par3-1)*(m*n)-(par2-1)*m;

        if((rows(i)-par1)^2+(cols(i)-par2)^2<(r_min+(par3-1)*step_r)^2+5&...

                (rows(i)-par1)^2+(cols(i)-par2)^2>(r_min+(par3-1)*step_r)^2-5)

            hough_circle(rows(i),cols(i)) = 1;

        end

    end

end



% 打印结果

for k=1:length

    par3 = floor(index(k)/(m*n))+1;

    par2 = floor((index(k)-(par3-1)*(m*n))/m)+1;

    par1 = index(k)-(par3-1)*(m*n)-(par2-1)*m;

    par3 = r_min+(par3-1)*step_r;

    fprintf(1,'Center %d %d radius %d\n',par1,par2,par3);

    para(:,k) = [par1,par2,par3];

end

图像融合就是将不同来源的同一对象的图像数据进行空间配准，然后采用一定的算法将各图像数据中所含的信息优势或互补性有机地结合起来产生新图像数据的信息技术。这种新数能够较全面地描述被研究对象，同单一信息源相比，能减少或抑制对被感知对象或环境解释中可能存在的多义性、不完全性、不确定性和误差，最大限度地利用各种信息源提供的信息。
一个简单的基于小波变换的图像融合程序，其基本原理是将两幅分别进行小波分解，将分解后的低频分量和高频分量分别叠加，然后将叠加后的高频分量和低频分量进行重构。
%导入待融合图像1
load bust
X1=X;
map1=map;
subplot(131);image(X1);
colormap(map1);title('原始图像1');
axis square
%导入待融合图像2
load mask
X2=X;
map2=map;
%对灰度值大于100的像素进行增强，小于100的像素进行减弱
for i=1:256
for j=1:256
if(X2(i,j)>100)
X2(i,j)=1.2*X2(i,j);
else
X2(i,j)=0.5*X2(i,j);
end
end
end
subplot(132)
image(X2);colormap(map2);title('原始图像2');
axis square
%对原始图像1进行小波分解
[c1,s1]=wavedec2(X1,2,'sym4');
%对分解后的低频部分进行增强
sizec1=size(c1);
for I=1:sizec1(2)
c1(I)=1.2*c1(I);
end
%对原始图像2进行分解
[c2,s2]=wavedec2(X2,2,'sym4');
%将分解后的低频分量和高频分量进行相加，并乘以权重系数0.5
c=c1+c2;
c=0.5*c;
s=s1+s2;
s=0.5*s;
%进行小波重构
xx=waverec2(c,s,'sym4');
subplot(133);image(xx);title('融合图像');
axis square

图像的去噪处理方法可分为空间域法和变换域法两大类。前者即是在原图像上直接进行数据运算，对像素的灰度值进行处理。变换域法是在图像的变换域上进行处理，对变换后的系数进行相应的处理，然后进行反变换达到图像去噪的目的。
（1） 基于离散余弦变换的图像去噪
基于离散余弦变换对图像的噪声抑制原理如下：一般而言，我们认为图像的噪声在离散余弦变换结果中处在其高频部分，而高频部分的幅值一般很小，利用这一性质，就很容易实现图像的噪声抑制。当然，这会同时失去图像的部分细节。例程1为基于离散余弦变换的图像去噪的MATLAB源程序。
例程1
*********************************************************************
%读取图像
X=imread('wangshi.jpg');
X=rgb2gray(X);
%读取图像尺寸
[m,n]=size(X);
%给图像加噪
Xnoised=imnoise(X,'speckle',0.01);
%输出加噪图像
figure(1);
imshow(Xnoised);
%DCT变换
Y=dct2(Xnoised);
I=zeros(m,n);
%高频屏蔽
I(1:m/3,1:n/3)=1;
Ydct=Y.*I;
%逆DCT变换
Y=uint8(idct2(Ydct));
%结果输出
figure(2);
imshow(Y);
*********************************************************************

（2） 基于小波变换的图像去噪
小波去噪是小波变换较为成功的一类应用，其去噪的基本思路:带噪信号经过预处理，然后利用小波变换把信号分解到各尺度中，在每一尺度下把属于噪声的小波系数去掉，保留并增强属于信号的小波系数，最后再经过小波逆变换恢复检测信号。

*********************************************************************
clear;
X=imread('life.jpg');
X=rgb2gray(X);
subplot(221);
imshow(X);
title('原始图像');
% 生成含噪图像并图示
init=2055615866;
randn('seed',init);
X=double(X);
% 添加随机噪声
XX=X+8*randn(size(X));
subplot(222);
imshow(uint8(XX));
title(' 含噪图像 ');
%用小波函数coif2对图像XX进行2层
% 分解
[c,l]=wavedec2(XX,2,'coif2');
% 设置尺度向量
n=[1,2];
% 设置阈值向量 , 对高频小波系数进行阈值处理
p=[10.28,24.08];
nc=wthcoef2('h',c,l,n,p,'s');
% 图像的二维小波重构
X1=waverec2(nc,l,'coif2');
subplot(223);
imshow(uint8(X1));
%colormap(map);
title(' 第一次消噪后的图像 ');
%再次对高频小波系数进行阈值处理
mc=wthcoef2('v',nc,l,n,p,'s');
% 图像的二维小波重构
X2=waverec2(mc,l,'coif2');
subplot(224);
imshow(uint8(X2));
title(' 第二次消噪后的图像 ');
*********************************************************************

不变矩(Invariant moments)是一种高度浓缩的图像特征，具有平移、灰度、尺度、旋转等多畸变不变性，因此矩和矩函数被广泛用于图像的模式识别、图像分类、目标识别和场景分析中。M．K．Hu在 1961年首先提出不变矩的概念，并将几何矩(Geometric moments，GMg)用于图像描述。
求解Hu的七个不变矩的Matlab源程序
function inv_m7 = invariable_moment(in_image)
%%%========================================================
format long
image=rgb2gray(in_image); %将输入的RGB图像转换为灰度图像
image=double(image); %将图像矩阵的数据类型转换成双精度型
%%%=================计算 、 、 =========================
m00=sum(sum(image)); %计算灰度图像的零阶几何矩
m10=0;
m01=0;
[row,col]=size(image);
for i=1:row
for j=1:col
m10=m10+i*image(i,j);
m01=m01+j*image(i,j);
end
end
%%%=================计算 、 ================================
u10=m10/m00;
u01=m01/m00;
%%%=================计算图像的二阶几何矩、三阶几何矩============
m20 = 0;m02 = 0;m11 = 0;m30 = 0;m12 = 0;m21 = 0;m03 = 0;
for i=1:row
for j=1:col
m20=m20+i^2*image(i,j);
m02=m02+j^2*image(i,j);
m11=m11+i*j*image(i,j);
m30=m30+i^3*image(i,j);
m03=m03+j^3*image(i,j);
m12=m12+i*j^2*image(i,j);
m21=m21+i^2*j*image(i,j);
end
end
%%%=================计算图像的二阶中心矩、三阶中心矩============
y00=m00;
y10=0;
y01=0;
y11=m11-u01*m10;
y20=m20-u10*m10;
y02=m02-u01*m01;
y30=m30-3*u10*m20+2*u10^2*m10;
y12=m12-2*u01*m11-u10*m02+2*u01^2*m10;
y21=m21-2*u10*m11-u01*m20+2*u10^2*m01;
y03=m03-3*u01*m02+2*u01^2*m01;
%%%=================计算图像的归格化中心矩====================
n20=y20/m00^2;
n02=y02/m00^2;
n11=y11/m00^2;
n30=y30/m00^2.5;
n03=y03/m00^2.5;
n12=y12/m00^2.5;
n21=y21/m00^2.5;
%%%=================计算图像的七个不变矩======================
h1 = n20 + n02; h2 = (n20-n02)^2 + 4*(n11)^2;
h3 = (n30-3*n12)^2 + (3*n21-n03)^2; h4 = (n30+n12)^2 + (n21+n03)^2;
h5 = (n30-3*n12)*(n30+n12)*((n30+n12)^2-3*(n21+n03)^2)+(3*n21-n03)*(n21+n03)*(3*(n30+n12)^2-(n21+n03)^2);
h6 = (n20-n02)*((n30+n12)^2-(n21+n03)^2)+4*n11*(n30+n12)*(n21+n03);
h7 = (3*n21-n03)*(n30+n12)*((n30+n12)^2-3*(n21+n03)^2)+(3*n12-n30)*(n21+n03)*(3*(n30+n12)^2-(n21+n03)^2);

inv_m7= [h1 h2 h3 h4 h5 h6 h7];
*********************************************************************

图像压缩的目的就是改变图像的表达方式，以尽量少的比特数表征图像，同时保持复原图像的质量。
于小波变换进行图像压缩的基本原理是：根据二维小波分解算法，一幅图像做小波分解后，可得到一些列不同分辨率的图像，而表现一幅图像最主要的部分是低频部分，如果去掉图像的高频部分而只保留低频部分，则可以达到图像压缩的目的。基于小波分析的图像压缩方法有很多，比较成功的有小波包最优基方法，小波域纹理模型方法，小波变换零树压缩，小波变换向量量化压缩等。
X=imread('robot.jpg');
X=rgb2gray(X);
X1=X;
% 分解图像，提取分解结构中的第一层系数
[c,l]=wavedec2(X,2,'bior3.7');
cA1=appcoef2(c,l,'bior3.7',1);
cH1=detcoef2('h',c,l,1);
cD1=detcoef2('d',c,l,1);
cV1=detcoef2('v',c,l,1);

% 重构第一层系数

A1=wrcoef2('a',c,l,'bior3.7',1);
H1=wrcoef2('h',c,l,'bior3.7',1);
D1=wrcoef2('d',c,l,'bior3.7',1);
V1=wrcoef2('v',c,l,'bior3.7',1);
c1=[A1 H1;V1 D1];
subplot(221),imshow(X1),title('原始图像'); axis square;
subplot(222),image(c1);title('分解后的高频和低频信息');
axis square

% 对图像进行压缩，保留第一层低频信息并对其进行量化编码
ca1=wcodemat(cA1,440,'mat',0);
ca1=0.5*ca1;
subplot(223);image(ca1);
axis square;
title('第一次压缩图像的大小:');
% 压缩图像，保留第二层低频信息并对其进行量化编码
cA2=appcoef2(c,l,'bior3.7',2);
ca2=wcodemat(cA2,440,'mat',0);
ca2=0.5*ca2;
subplot(224);
image(ca2);
title('第二次压缩图像');
axis square;

图像边缘是图像的重要特征，是图像中特性（如像素灰度、纹理等）分布的不连续处，图像周围特性有阶跃变化或屋脊状变化的那些像素集合。图像的边缘部分集中了图像的大部分信息，一幅图像的边缘结构与特点往往是决定图像特质的重要部分。图像边缘的另一个定义是指其周围像素灰度变化不连续的那些像素的集合。边缘广泛存在于物体与背景之间、物体与物体之间，因此，边缘是图像分割、图像理解及图像识别的重要特征。 图像边缘检测主要用于增强图像中的轮廓边缘、细节以及灰度跳变部分，形成完整的物体边界，达到将物体从图像中分离出来或将表示同一物体表面的区域检测出来的目的。目前为止最通用的方法是检测亮度值的不连续性，用一阶和二阶导数检测的。 （1）微分法 微分法的目的是利用微分运算求信号的变化率，加强高频分量的作用，从而使轮廓清晰。遵循如下两个基本准则之一：找到亮度的一阶导数在幅度上比指定的阈值大的地方；找到亮度的二阶导数有零交叉的地方。 （2）差分边缘检测方法 利用像素灰度的一阶导数算子在灰度迅速变化处得到高值来进行奇异点的检测。它在某一点的值就代表该点的边缘强度，通过对这些值设置阈值来进一步得到边缘图像。差分边缘检测方法是最原始、最基本的方法。但要求差分方向与边缘方向垂直，这就需要对图像的不同方向（一般为垂直方向、水平方向和对角线方向）都进行差分运算，增加了实际运算的繁琐性，目前很少采用。 （3）Roberts边缘检测算子 Roberts边缘检测算子根据任意一对互相垂直方向上的差分可用来计算梯度的原理，采用对角线方向相邻两像素之差然后计算出Roberts的梯度幅度值. Roberts检测器较为简单，但具有一些功能上的限制，例如，它是非对称的，而且不能检测诸如45°倍数的边缘。然而，它还是经常用于硬件实现中，因为它既简单又快速。 （4）Sobel边缘检测算子 对数字图像的每个像素，考察它上下左右邻点灰度的加权差，与之接近的邻点的权大。 Sobel算子很容易在空间上实现，边缘检测效果较好，且受噪声的影响也较小。邻域增大抗噪性会更好，但计算量也会增大，得出的边缘也会相应变粗。Sobel算子会检测出许多伪边缘，边缘定位精度不够高，在精度要求不高时是一种较常用的边缘检测方法。 （5）Prewitt边缘检测算子 图像中的每个像素都用这两个核作卷积，一个核对垂直边缘影响最大，另一个核对水平边缘影响最大。两个卷积的绝对值的最大值作为该点的输出值。不能简单的将小于0的值处理为0，这样会丢失信息。它比Sobel检测器在计算上要简单一些，但比较容易产生一些噪声。 （6）拉普拉斯边缘检测算子 拉普拉斯边缘检测算子是一种二阶微分算子，与其它边缘检测方法的不同之处在于，该方法是一种各向同性的检测方法，即其边缘的增强程度与边缘的方向无关，从而可以满足不同走向的边缘锐化的要求。 拉普拉斯算子自身很少被直接用作边缘检测，因为二阶导数对噪声具有无法接受的敏感性，它的幅度会产生双边缘，而且它不能检测边缘的方向。然而，当与其它边缘检测技术组合使用时，拉普拉斯算子是一种有效的补充方法。例如，虽然它的双边缘使得它不适合直接用于边缘检测，但该性质可用于边缘定位。
在MATLAB中可以由edge函数实现各算子对边缘的检测，其调用格式为：         BW = edge(I,’METHOD’)—自动选择阈值用指定的算子进行边缘检测； BW = edge(I,’METHOD’,THRESH)—根据给定的阈值THRESH用指定的算子进行边缘检测，忽略所有小于阈值的边缘。当THRESH为空时，自动选择阈值； [BW, THRESH] = edge(I,’METHOD’,…)—返回edge使用的阈值，以确定哪个梯度值足够大到可以称为边缘点。 其中,I为输入图像。返回图像BW为与I同大的二值图像，1表示边缘，0表示非边缘。I是unit8型、unit16型，或者是double型，BW是unit8型。         METHOD：使用检测算子的类型，经常使用的有： sobel：缺省值，用导数的sobel近似值检测边缘，那些梯度最大的点返回边缘。         roberts：用导数的roberts近似值检测边缘，那些梯度最大的点返回边缘。         prewitt：用导数的prewitt近似值检测边缘，那些梯度最大的点返回边缘。         LoG：用LoG算子检测边缘。         Canny：用Canny算子检测边缘。 THRESH：指定的阈值，所有不强于thresh的边缘都被忽略。
由edge函数实现各算子对图像的边缘检测 clear all; I = imread('d:\office.bmp'); I=rgb2gray(I); BW1 = edge(I,'sobel');   %利用Sobel算子进行边缘检测 BW2 = edge(I,'roberts');  %利用roberts算子进行边缘检测 BW3 = edge(I,'prewitt');  %利用prewitt算子进行边缘检测 BW4 = edge(I,'log');     %利用log算子进行边缘检测 BW5 = edge(I,'canny');   %利用canny算子进行边缘检测 subplot(2,3,1),imshow(I) subplot(2,3,2),imshow(BW1) subplot(2,3,3),imshow(BW2) subplot(2,3,4),imshow(BW3) subplot(2,3,5),imshow(BW4) subplot(2,3,6),imshow(BW5)	回复此楼 支持

图像增强是对图像进行加工，以得到视觉效果更好或更有用的新图像。图像增强在处理方法上可分为基于空域的图像增强和基于频域的图像增强。
空域图像增强常用方法：
（1）直方图均衡化
   直方图是一个离散函数，它表示数字图像每一灰度级与该灰度级出现频率的对应关系。有的图像的灰度直方图在低值灰度区间上频率较大，使得图像中较暗区域中的细节常常看不清楚。为了使图像清晰，可将图像灰度范围拉开，并且让灰度频率较小的灰度级变大，即让灰度直方图在较大的动态范围内趋于一致。
（2）灰度调整
  灰度调整时按一定的规则修改输入图像每一个像素的灰度，从而改变图像灰度的动态范围。它可对灰度范围进行扩展或压缩，也可对图像进行分段处理，根据图像的特点和要求在某段区间中进行压缩而在另外区间中进行扩展。
（3）空域滤波
    空域滤波可用来对图像的某些特征进行增强，而去除其他特征。空域滤波是一种模板操作，即输出图像的某一像素值由输入图像相应的像素及其邻域的像素值决定。模板操作可看成图像和模板的卷积。而不同的滤波操作，其模板就不同。

计算和显示直方图的MATLAB函数为imhist，其语法格式如下：
imhist(I,n);
其中参数I是待处理的灰度图像，由imread函数读入。参数n是指定的灰度级数目，默认值为256。
计算和显示索引色图像直方图：
imhist(I,map);
其中参数I是待处理的索引色图像，由imread函数读入。参数map是调色板。
     图像的直方图均衡化可用MATLAB中的histeq函数实现，其语法格式如下：
将原始图像的直方图变成用户指定的向量：
J=histeq(I,hgram);
其中参数I是待处理的灰度图像，由imread函数读入。参数hgram是用户指定的向量，向量中的每个元素要在0和1之间。
指定直方图均衡化后的灰度级数：
J=histeq(I,n);
其中参数I是待处理的灰度图像，由imread函数读入。参数n是用户指定的均衡化后的灰度级数，默认值为64。

图像的灰度调整可用MATLAB中的imadjust函数实现，其语法格式如下：
J=imadjust (I,[low  high],[bottom  top],gamma);
其中参数I是待处理的灰度图像，由imread函数读入。参数[low  high]是原图像要变换的灰度范围。向量中的每个元素要在0和1之间。参数[bottom  top]是变换后的灰度范围。向量中的每个元素要在0和1之间。参数gamma是矫正量。

图像的空域滤波可用MATLAB中的filter2函数实现，其语法格式如下：
Im=filter2(h,I,shape);
其中参数I是待处理的灰度图像，由imread函数读入。参数h是滤波算子（模板）。参数shape是指定滤波的范围，shape=’full’是做图像边界补零；shape=’same’是输出的图像Im和输入的图像I大小相同；shape=’valid’是不考虑边界补零，只计算有效输出部分。该参数可省略不写。
滤波模板h可以自己定义，也可调用MATLAB的fspecial函数获得，该函数语法格式如下：
h=fspecial(type,para);
参数type指定算子类型，参数para指定相应的参数，具体意义如下：
其一是type=’average’，为均值滤波器，对应参数为n，代表模板的尺寸，用向量表示，默认值为[3 3]。
其二是type=’laplacian’，为拉普拉斯算子，对应参数为alpha，用于控制算子的形状，取值范围为0到1，默认值为0.2。
其三是type=’prewit’，为Prewitt算子，用于边缘增强，无参数。
其四是type=’sobel’，为Sobel算子，用于边缘提取，无参数。
其五是type=’unsharp’，为对比度增强滤波器，对应参数为alpha，用于控制滤波器的形状，取值范围为0到1，默认值为0.2。

（2） 基于小波变换的图像去噪
小波去噪是小波变换较为成功的一类应用，其去噪的基本思路:带噪信号经过预处理，然后利用小波变换把信号分解到各尺度中，在每一尺度下把属于噪声的小波系数去掉，保留并增强属于信号的小波系数，最后再经过小波逆变换恢复检测信号。

*********************************************************************
clear;
X=imread('life.jpg');
X=rgb2gray(X);
subplot(221);
imshow(X);
title('原始图像');
% 生成含噪图像并图示
init=2055615866;
randn('seed',init);
X=double(X);
% 添加随机噪声
XX=X+8*randn(size(X));
subplot(222);
imshow(uint8(XX));
title(' 含噪图像 ');
%用小波函数coif2对图像XX进行2层
% 分解
[c,l]=wavedec2(XX,2,'coif2');
% 设置尺度向量
n=[1,2];
% 设置阈值向量 , 对高频小波系数进行阈值处理
p=[10.28,24.08];
nc=wthcoef2('h',c,l,n,p,'s');
% 图像的二维小波重构
X1=waverec2(nc,l,'coif2');
subplot(223);
imshow(uint8(X1));
%colormap(map);
title(' 第一次消噪后的图像 ');
%再次对高频小波系数进行阈值处理
mc=wthcoef2('v',nc,l,n,p,'s');
% 图像的二维小波重构
X2=waverec2(mc,l,'coif2');
subplot(224);
imshow(uint8(X2));
title(' 第二次消噪后的图像 ');
*********************************************************************

不变矩(Invariant moments)是一种高度浓缩的图像特征，具有平移、灰度、尺度、旋转等多畸变不变性，因此矩和矩函数被广泛用于图像的模式识别、图像分类、目标识别和场景分析中。M．K．Hu在 1961年首先提出不变矩的概念，并将几何矩(Geometric moments，GMg)用于图像描述。
求解Hu的七个不变矩的Matlab源程序
function inv_m7 = invariable_moment(in_image)
%%%========================================================
format long
image=rgb2gray(in_image); %将输入的RGB图像转换为灰度图像
image=double(image); %将图像矩阵的数据类型转换成双精度型
%%%=================计算 、 、 =========================
m00=sum(sum(image)); %计算灰度图像的零阶几何矩
m10=0;
m01=0;
[row,col]=size(image);
for i=1:row
for j=1:col
m10=m10+i*image(i,j);
m01=m01+j*image(i,j);
end
end
%%%=================计算 、 ================================
u10=m10/m00;
u01=m01/m00;
%%%=================计算图像的二阶几何矩、三阶几何矩============
m20 = 0;m02 = 0;m11 = 0;m30 = 0;m12 = 0;m21 = 0;m03 = 0;
for i=1:row
for j=1:col
m20=m20+i^2*image(i,j);
m02=m02+j^2*image(i,j);
m11=m11+i*j*image(i,j);
m30=m30+i^3*image(i,j);
m03=m03+j^3*image(i,j);
m12=m12+i*j^2*image(i,j);
m21=m21+i^2*j*image(i,j);
end
end
%%%=================计算图像的二阶中心矩、三阶中心矩============
y00=m00;
y10=0;
y01=0;
y11=m11-u01*m10;
y20=m20-u10*m10;
y02=m02-u01*m01;
y30=m30-3*u10*m20+2*u10^2*m10;
y12=m12-2*u01*m11-u10*m02+2*u01^2*m10;
y21=m21-2*u10*m11-u01*m20+2*u10^2*m01;
y03=m03-3*u01*m02+2*u01^2*m01;
%%%=================计算图像的归格化中心矩====================
n20=y20/m00^2;
n02=y02/m00^2;
n11=y11/m00^2;
n30=y30/m00^2.5;
n03=y03/m00^2.5;
n12=y12/m00^2.5;
n21=y21/m00^2.5;
%%%=================计算图像的七个不变矩======================
h1 = n20 + n02; h2 = (n20-n02)^2 + 4*(n11)^2;
h3 = (n30-3*n12)^2 + (3*n21-n03)^2; h4 = (n30+n12)^2 + (n21+n03)^2;
h5 = (n30-3*n12)*(n30+n12)*((n30+n12)^2-3*(n21+n03)^2)+(3*n21-n03)*(n21+n03)*(3*(n30+n12)^2-(n21+n03)^2);
h6 = (n20-n02)*((n30+n12)^2-(n21+n03)^2)+4*n11*(n30+n12)*(n21+n03);
h7 = (3*n21-n03)*(n30+n12)*((n30+n12)^2-3*(n21+n03)^2)+(3*n12-n30)*(n21+n03)*(3*(n30+n12)^2-(n21+n03)^2);

inv_m7= [h1 h2 h3 h4 h5 h6 h7];
*********************************************************************

图像压缩的目的就是改变图像的表达方式，以尽量少的比特数表征图像，同时保持复原图像的质量。
于小波变换进行图像压缩的基本原理是：根据二维小波分解算法，一幅图像做小波分解后，可得到一些列不同分辨率的图像，而表现一幅图像最主要的部分是低频部分，如果去掉图像的高频部分而只保留低频部分，则可以达到图像压缩的目的。基于小波分析的图像压缩方法有很多，比较成功的有小波包最优基方法，小波域纹理模型方法，小波变换零树压缩，小波变换向量量化压缩等。
X=imread('robot.jpg');
X=rgb2gray(X);
X1=X;
% 分解图像，提取分解结构中的第一层系数
[c,l]=wavedec2(X,2,'bior3.7');
cA1=appcoef2(c,l,'bior3.7',1);
cH1=detcoef2('h',c,l,1);
cD1=detcoef2('d',c,l,1);
cV1=detcoef2('v',c,l,1);

% 重构第一层系数

A1=wrcoef2('a',c,l,'bior3.7',1);
H1=wrcoef2('h',c,l,'bior3.7',1);
D1=wrcoef2('d',c,l,'bior3.7',1);
V1=wrcoef2('v',c,l,'bior3.7',1);
c1=[A1 H1;V1 D1];
subplot(221),imshow(X1),title('原始图像'); axis square;
subplot(222),image(c1);title('分解后的高频和低频信息');
axis square

% 对图像进行压缩，保留第一层低频信息并对其进行量化编码
ca1=wcodemat(cA1,440,'mat',0);
ca1=0.5*ca1;
subplot(223);image(ca1);
axis square;
title('第一次压缩图像的大小:');
% 压缩图像，保留第二层低频信息并对其进行量化编码
cA2=appcoef2(c,l,'bior3.7',2);
ca2=wcodemat(cA2,440,'mat',0);
ca2=0.5*ca2;
subplot(224);
image(ca2);
title('第二次压缩图像');
axis square;