原文：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倍。
利用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。
利用imwrite完成图像的输出和保存操作，完全支持和支持上述各种操作
图像文件的格式。语法格式如下：
imwrite(I,‘filename’,‘fmt’)或imwrite(I,‘filename.fmt’)；
其中的I、filename和fmt上述意义。
注：当使用时imwrite当函数保存图像时，MATLAB默认保存方法简化为uint8数据类型。类似于读取文件的类型，MATLAB文件保存时还支持16位PNG和TIFF因此，当用户保存此类文件时，MATLAB存储它uint16中。
图像的实际过程是将数字图像从一组离散数据还原为可见图像的过程。
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个区域作为当前绘图区。
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。
的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）空域滤波
空域滤波可用来对图像的某些特征进行增强，而去除其他特征。空域滤波是一种模板操作，即输出图像的某一像素值由输入图像相应的像素及其邻域的像素值决定。模板操作可看成图像和模板的卷积。而不同的滤波操作，其模板就不同。
计算和显示直方图的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);