MatLab代码大全

图像获取

2.3.2 二维连续傅里叶变换

例2.2

figure(1); %建立图形窗口1

[u,v] = meshgrid(-1:0.01:1); %生成二维频域网格

F1 = abs(sinc(u.*pi));

F2 = abs(sinc(v.*pi));

F=F1.*F2; %计算幅度频谱F=|F(u,v)|

surf(u,v,F);

%显示幅度频谱，如图2.3(b)

shading interp;

%平滑三维曲面上的小格

axis off; %关闭坐标系

figure(2); %建立图形窗口2

F1=histeq(F); %扩展F的对比度以增强视觉效果

imshow(F1); %用图像来显示幅度频谱，如图2.3(c)

图像变换

3.4.4 二维FFT的MATLAB实现

例3.2

简单图像及其傅里叶变换

MATLAB程序：

%建立简单图像d并显示之

d = zeros(32,32); %图像大小3232

d(13:20,13:20) = 1; %中心白色方块大小为88

figure(1); %建立图形窗口1

imshow(d,'notruesize')； %显示图像d如图3.5(a)所示

%计算傅里叶变换并显示之

D = fft2(d); %计算图像d的傅里叶变换，fft2(d) = fft(fft(d).').'

figure(2); %建立图形窗口2

imshow(abs(D),[-1 5],'notruesize'); %显示图像d的傅里叶变换谱如3.5(b)所示

例3.3 MATLAB图像及其傅里叶变换谱

MATLAB程序：

figure(1);

load imdemos saturn2; %装入MATLAB图像saturn2

imshow(saturn2); %显示图像saturn2如图3.6(a)所示

figure(2);

S= fftshift(fft2(saturn2)); %计算傅里叶变换并移位

imshow(log(abs(S)),[ ]); %显示傅里叶变换谱如3.6(b)所示

例3.4

真彩图像及其傅里叶变换谱

MATLAB程序：

figure(1);

A=imread(''); %装入真彩图像，见图1.1(b)

B=rgb2gray(A); %将真彩图像转换为灰度图像

8

imshow(B); %显示灰度图像如图3.7(a)所示

C=fftshift(fft2(B)); %计算傅里叶变换并移位

figure(2);

imshow(log(abs(C)),[ ]); %显示傅里叶变换谱如3.7(b)所示

3.5.4 离散余弦变换的MATLAB实现

例3.5

计算并显示真彩图像余弦变换的MATLAB程序如下：

RGB=imread(''); %装入真彩图像

figure(1);

imshow(RGB); %显示彩色图像

GRAY=rgb2gray(RGB); %将真彩图像转换为灰度图像

figure(2);

imshow(GRAY); %显示灰度图像如图3.10(a)所示

DCT=dct2(GRAY); %进行余弦变换

figure(3);

imshow(log(abs(DCT)),[ ]); %显示余弦变换如图3.10(b)所示。

3.8.2 Radon变换的MATLAB实现

例3.8真彩图像的Radon变换

MATLAB程序如下：

RGB=imread(''); %装入真彩图像

GRAY=rgb2gray(RGB); %将真彩图像转换为灰度图像

figure(2);

imshow(GRAY); %显示灰度图像如图3.16(a)

[R,xp] = radon(GRAY,[0 45]); %计算变换角度为0°和45°的Radon变换

figure; plot(xp,R(:,1)); title('R_{0^o} (xprime)')

%显示0°方向上的Radon变换如图3.16(b)

figure; plot(xp,R(:,2)); title('R_{45^o} (xprime)')

%显示45°方向上的Radon变换如图3.16(c)

例3.9连续角度的Radon变换

对于一组连续角度的Radon变换通常用一幅图像来表示。本例先建立一幅简单图像，然后令变换角度从0° 以1°的增量变化到180°时的Radon变换情况。其MATLAB程序如下：

I = zeros(100,100);

%建立简单图像如图3.17(a)

I(25:75, 25:75) = 1;

figure(1);imshow(I);

theta = 0:180; %规定变换角度的范围

[R,xp] = radon(I,theta); %计算Radon变换

figure(2);

imagesc(theta,xp,R); %以图像方式显示变换结果R，

%其x轴和y轴分别为theta和xp

title(‘R_{theta} (Xprime)’); %显示图像标题R(x')

9

(degrees)xlabel(‘theta (degrees)’); %显示x坐标“”

ylabel(‘Xprime’); %显示y坐标“x'”

set(gca,’Xtick’,0:20:180); %设置x坐标刻度

colormap(hot); %设置调色板

colorbar； %显示当前图像的调色板

第4章 图像增强

4.2.1 直接灰度变换

Matlab程序实现图像求反：

I = imread(' ');

imshow(I)

I=double(I)

I=256-1-I

I=uint8(I)

figure

imshow(I)

例4.1

用Matlab程序实现线性灰度变换的图像增强：

%读入并显示原始图像

I = imread('');

imshow(I);

I=double(I);

[M,N]=size(I);

%进行线性灰度变换

for i=1:M

for j=1:N

if I(i,j)<=30

I(i,j)=I(i,j);

elseif I(i,j)<=150

I(i,j)=(200-30)/(150-30)*(I(i,j)-30)+30;

else

I(i,j)=(255-200)/(255-150)*(I(i,j)-150)+200;

end

end

end

例4.2

I=imread('');

figure;imshow(I);

I=double(I);

I2=41*log(1+I);

I2=uint8(I2);

figure;imshow(I2);

例4.3

灰度切割变换的Matlab的程序如下：

I=imread('');

figure;imshow(I);

I=double(I)

[M,N]=size(I);

10

for i=1:M

for j=1:N

if I(i,j)<=50

I(i,j)=40;

elseif I(i,j)<=180

I(i,j)=220;

else

I(i,j)=40;

end

end

end

I=uint8(I);

figure;imshow(I);

例4.4

具体Matlab程序如下：

I=imread('');

imshow(I);

I=double(I);

[M,N]=size(I);

for k=1:8

J=zeros(M,N);

for i=1:M

for j=1:N

temp=I(i,j);

s1=0;s2=0;

range=[k:-1:1];

for d=range

s1=2^(8-d)+s1;s2=2^(8-d+1);

if temp>=s1 & temp

J(i,j)=255; break;

end

end

end

end

J=uint8(J);

figure;imshow(J);

end

4.2.2 直方图修正

例4.6直方图均衡化效果实例

用Matlab中的histeq函数实现直方图均衡化的程序如下：I=imread('');

figure

subplot(221);imshow(I);

subplot(222);imhist(I)

I1=histeq(I);

figure;

subplot(221);imshow(I1)

subplot(222);imhist(I1)

例4.8：直方图规定效果实例

11

用matlab中的histeq函数实现直方图均衡化的程序如下：

I=imread('');

[M,N]=size(I);

for i=1:8:257

counts(i)= i;

end

Q=imread('');

N=histeq(Q,counts);

figure

subplot(221);imshow(N);

subplot(222);imhist(N);

axis([0 260 0 5000]);

4.2.3 图像间运算

例：用图像平均减少随机噪声

I=imread('');

[M,N]=size(I);

II1=zeros(M,N);

for i=1:16

II(:,:,i)=imnoise(I,'gaussian',0,0.01);

II1=II1+double(II(:,:,i));

if or(or(i==1,i==4),or(i==8,i==16));

figure;

imshow(uint8(II1/i));

end

end

4.3 空域滤波增强

Matlab实现的邻域平均法抑制噪声的程序：

I=imread('');

J=imnoise(I,'salt & pepper', 0.02);

subplot(231),imshow(I);title('原图像');

subplot(232),imshow(J);title('添加椒盐噪声图像')

k1=filter2(fspecial('average',3),J); %进行3×3模板平滑滤波

k2=filter2(fspecial('average',5),J); %进行5×5模板平滑滤波

k3=filter2(fspecial('average',7),J); %进行7×7模板平滑滤波

k4=filter2(fspecial('average',9),J); %进行9×9模板平滑滤波

subplot(233),imshow(uint8(k1));title('3×3模板平滑滤波');

subplot(234),imshow(uint8(k2));title('5×5模板平滑滤波');

subplot(235),imshow(uint8(k3));title('7×7模板平滑滤波');

subplot(236),imshow(uint8(k4));title('9×9模板平滑滤波')

例4.10：使用中值滤波降低图像噪声

I=imread('');

J=imnoise(I,'salt & pepper', 0.02);

subplot(231),imshow(I);title('原图像');

subplot(232),imshow(J);title('添加椒盐噪声图像')

k1=medfilt2(J); %进行3×3模板中值滤波

k2=medfilt2(J,[5 5]); %进行5×5模板中值滤波

12

k3=medfilt2(J,[7 7]); %进行7×7模板中值滤波

k4=medfilt2(J,[9 9]); %进行9×9模板中值滤波

subplot(233),imshow(k1);title('3×3模板中值滤波')

subplot(234),imshow(k2);title('5×5模板中值滤波')

subplot(235),imshow(k3);title('7×7模板中值滤波')

subplot(236),imshow(k4);title('9×9模板中值滤波')

例4.11：梯度锐化实例

I=imread('');

subplot(131),imshow(I)

H=fspecial('Sobel');

H=H'; %Sobel垂直模板

TH=filter2(H,I);

subplot(132),imshow(TH,[]);

H=H'; %Sobel水平模板

TH=filter2(H,I);

subplot(133),imshow(TH,[])

4.4 图像频域增强

例4.12：频域低通滤波所产生的模糊

%理想低通过滤波器所产生的模糊和振铃现象

J=imread('');

subplot(331);imshow(J);

J=double(J);

%采用傅里叶变换

f=fft2(J);

%数据矩阵平衡

g=fftshift(f);

subplot(332);imshow(log(abs(g)),[]),color(jet(64));

[M,N]=size(f);

n1=floor(M/2);

n2=floor(N/2);

% d0=5,15,45,65

d0=5;

for i=1:M

for j=1:N

d=sqrt((i-n1)^2+(j-n2)^2);

if d<=d0

h=1;

else

h=0;

end

g(i,j)=h*g(i,j);

end

end

g=ifftshift(g);

g=uint8(real(ifft2(g)));

subplot(333);

imshow(g);

例4.12：用巴特沃斯低通滤波器去除图像中的盐椒噪声

13

%实现Butterworth低通过滤波器

I=imread('');

J=imnoise(I,'salt & pepper',0.02); %给原图像加入椒盐噪声，如图4.33(a)所示

subplot(121);imshow(J);

tilte('含有盐椒噪声的图像')

J=double(J);

%采用傅里叶变换

f=fft2(J);

%数据矩阵平衡

g=fftshift(f)

[M,N]=size(f);

n=3;

d0=20

n1=floor(M/2)

n2=floor(N/2)

for i=1:M

for j=1:N

d=sqrt((i-n1)^2+(j-n2)^2)

h=1/(1+(d/d0)^(2*n));

g(i,j)=h*g(i,j);

end

end

g=ifftshift(g);

g=uint8(real(ifft2(g)));

subplot(122);

imshow(g); %结果如图4.33(b)所示

例4.13：频域高通滤波增强示例

J=imread('');

imshow(uint8(J));title('模糊图像')

J=double(J);

f=fft2(J); %采用傅里叶变换

g=fftshift(f);%数据矩阵平衡

[M,N]=size(f);

n1=floor(M/2);

n2=floor(N/2);

d0=20;

for i=1:M %进行理想高通滤波和理想高通加强滤波

for j=1:N

d=sqrt((i-n1)^2+(j-n2)^2);

if d>=d0

h1=1;

h2=1+0.5;

else

h1=0;

h2=0.5;

end

g1(i,j)=h1*g(i,j);

g2(i,j)=h2*g(i,j);

end

end

14

g1=ifftshift(g1);

g1=uint8(real(ifft2(g1)));

subplot(221);imshow(g1); %显示理想高通滤波结果

title('理想高通滤波结果');

g2=ifftshift(g2);

g2=uint8(real(ifft2(g2)));

subplot(222);imshow(g2); %显示理想高通加强滤波结果

title('理想高通加强滤波结果');

n=2;

d0=20;

for i=1:M %进行巴特沃斯高通滤波和巴特沃斯高通加强滤波

for j=1:N

d=sqrt((i-n1)^2+(j-n2)^2);

if d==0

h1=0;

h2=0.5;

else

h1=1/(1+(d0/d)^(2*n));

h2=1/(1+(d0/d)^(2*n))+0.5;

end

gg1(i,j)=h1*g(i,j);

gg2(i,j)=h2*g(i,j);

end

end

gg1=ifftshift(gg1);

gg1=uint8(real(ifft2(gg1)));

subplot(223);imshow(gg1); %显示巴特沃斯高通滤波结果

title('巴特沃斯高通滤波结果')

gg2=ifftshift(gg2);

gg2=uint8(real(ifft2(gg2)));

subplot(224);imshow(gg2); %显示巴特沃斯高通加强滤波结果

title('巴特沃斯高通加强滤波结果');

例4.14：同态滤波的增强效果

J=imread(''); %读入原图

subplot(121);imshow(J);

J=double(J);

f=fft2(J); %采用傅里叶变换

g=fftshift(f); %数据矩阵平衡

[M,N]=size(f);

d0=10;

rl=0.5;

rh=2

c=4；

n1=floor(M/2);

n2=floor(N/2);

for i=1:M

for j=1:N

d=sqrt((i-n1)^2+(j-n2)^2);

h=(rh-rl)*(1-exp(-c*(d.^2/d0.^2)))+rl;

15

g(i,j)=h*g(i,j);

end

end

g=ifftshift(g);

g=uint8(real(ifft2(g)));

subplot(122);imshow(g);

第5章 图像复原

例5.1

C=imread(''); %装入清晰图像

subplot(1,2,1)； %将图形窗口分成两个矩形平面

imshow(C); %在第一个矩形平面中显示装入的图像

LEN=30; %设置运动位移为30个象素

THETA=45; %设置运动角度为45o

PSF=fspecial('motion',LEN,THETA); %建立二维仿真线性运动滤波器PSF

MF=imfilter(C,PSF,'circular','conv'); %用PSF产生退化图像

subplot(1,2,2),imshow(MFUZZY); %在第二个矩形平面中显示模糊后的图像

imwrite(MF,''); %将运动模糊后的图像保存起来备用

例5.2

消除图5.4(b)的运动模糊，其MATLAB程序如下：

[MF,map]=imread(''); %装入运动模糊图像

figure(1);

imshow(MF); %显示模糊图像

LEN=30;

THETA=45;

INITPSF=fspecial('motion',LEN,THETA); %建立复原点扩散函数

[J P]= deconvblind(MF,INITPSF,30); %去卷积

figure(2); imshow(J); %显示结果图像如图5.6(a)

figure(3); imshow(P,[],'notruesize'); %显示复原点扩散函数如图5.6(b)

例5.3逆滤波与维纳滤波的比较

F = checkerboard(8); %生成原始图像F

figure(1);

imshow(F,[]);

PSF = fspecial('motion',7,45); %生成运动模糊图像MF

MF = imfilter(F,PSF,'circular');

noise = imnoise(zeros(size(F)),'gaussian',0,0.001); %生成高斯噪声

MFN = MF + noise; %生成运动模糊+高斯噪声图像MFN

figure(2);

imshow(MFN,[]);

NSR = sum(noise(:).^2)/sum(MFN(:).^2); %计算噪信比

figure(3);

imshow(deconvwnr(MFN,PSF),[]); %逆滤波复原

figure(4);

imshow(deconvwnr(MFN,PSF,NSR),[]); %维纳滤波复原

例5.7顺序统计滤波器比较

16

f=imread('');

figure(1); imshow(f);

title('原始图像');

g=imnoise(f,'salt & pepper',0.2);

figure(2); imshow(g);

title('椒盐噪声污染的图像');

g1=double(g)/255;

j1=medfilt2(g1,'symmetric');

figure(3); imshow(j1);

title('中值滤波图像');

j2=ordfilt2(g1,median(1:3*3),ones(3,3),'symmetric');

figure(4); imshow(j2);

title('中点滤波图像');

j3=ordfilt2(g1,1,ones(3,3));

figure(5); imshow(j3);

title('最小值滤波图像');

j4=ordfilt2(g1,9,ones(3,3));

figure(6); imshow(j4);

title('最大值滤波图像');

例5.8 简单图像的affine变换

MATLAB程序如下：

f=checkerboard(24); % 建立原始图像，如图5.16(a)

figure(1); imshow(f);

s=0.7;

theta=pi/6;

T=[s*cos(theta) s*sin(theta) 0 % 建立变换矩阵：旋转与尺度变换

-s*sin(theta) s*cos(theta) 0

0 0 1];

tform=maketform('affine',T);

g1=imtransform(f,tform,'nearest'); % 最近邻插值变换, 如图5.16(b)

figure(2); imshow(g1);

g2=imtransform(f,tform); % 双线性插值变换, 如图5.16(c)

figure(3); imshow(g2);

g3=imtransform(f,tform,'FillValue',0.5); % 修改双线性插值变换的背景色为灰色

figure(4); imshow(g3);

例5.9 利用“连接点”实施图像配准复原

MATLAB程序如下：

f=imread(''); % 读入256256原始图，如图5.17(a)

figure(1); imshow(f);

g=imread(''); % 读入几何失真图，，如图5.17(b)

figure(2); imshow(g);

17

%利用cpselect(g, f)函数交互选择“连接点”

base_points =[ 256.4000 256.1273; 1.5818 256.4182;

256.4182 1.0000; 200.5636 203.4727;

147.3273 183.9818; 96.4182 145.0000;

44.6364 35.0364; 157.5091 30.3818];

input_points =[ 280.0455 304.6182; 1.3455 255.2545;

255.2545 1.0000; 205.8545 225.8909;

145.7455 196.5273; 90.4727 146.7818;

38.6545 32.4364; 148.5091 31.0545];

tform = cp2tform(input_points, base_points, 'projective');

gp = imtransform(g, tform, 'XData', [1 256], 'YData', [1 256]);

figure(3); imshow(gp); % 显示复原图像，如图5.17(c)

第6章 彩色图像处理

例6.1 考虑生成一幅128*128的RGB图像，该图像左上角为红色，右上角为蓝色，左下角为绿色，右下角为黑色。其Matlab程序如下：

iR = zeros(128, 128); %生成一个128*128的零矩阵，作为R分量

iR(1:64, 1:64) = 1; %将左上角的64*64设置成1

iG = zeros(128, 128); %生成一个128*128的零矩阵，作为G分量

iG(65:128, 1:64) = 1; %将右下角的64*64设置成1

iB = zeros(128, 128); %生成一个128*128的零矩阵，作为B分量

iB(1:64, 65:128) = 1; %将右上角的64*64设置成1

I = cat(3, iR, iG, iB); %使用cat函数将三个分量组合

imshow(I) %显示生成的RGB图像，如图6.5所示

例6.2 将一幅RGB图像转换到CMY空间：

rgb_I = imread(''); %载入一幅彩色图像

cmy_I = imcomplement(rgb_I); %函数imcomplement转换到CMY空间

imshow(I); %显示原图，如图6.7（a）所示

figure, imshow(I2); %显示转换后图，如图6.7（b）所示

例6.3 将一幅三原色图像从RGB空间转换到HSI空间，其结果见图6.9。

rgb = imread('三原色.bmp'); %载入一幅图像

imshow(rgb); %见图6.9（a）

rgb = im2double(rgb); %将图像转换成double类型

r = rgb(:, :, 1); %提取图像的r分量

g = rgb(:, :, 2); %提取图像的g分量

b = rgb(:, :, 3); %提取图像的b分量

I = (r + g + b)/3; %计算I分量

%计算S分量

tmp1 = min(min(r, g), b);

tmp2 = r + g + b;

tmp2(tmp2 == 0) = eps;

18

S = 1 - 3.*tmp1./tmp2;

%计算H分量

tmp1 = 0.5*((r - g) + (r - b));

tmp2 = sqrt((r - g).^2 + (r - b).*(g-b));

theta = acos(tmp1./(tmp2 + eps));

H = theta;

H(b > g) = 2*pi - H(b > g);

H = H/(2*pi);

H(S == 0) = 0;

hsi = cat(3, H, S, I);

figure, imshow(H);

figure, imshow(S);

figure, imshow(I);

%见图6.9（b）

%见图6.9（c）

%见图6.9（d）

例6.4

rgb = imread(‘’); %载入一幅图像

imshow(rgb); %显示，见图6.10（a）

R = rgb(:, :, 1); %提取图像的R、G、B分量

G = rgb(:, :, 2);

B = rgb(:, :, 3);

figure, imshow(R); %分别显示图像的R、G、B分量。见图6.10(b)。

figure, imshow(G); %见图6.10(c)

figure, imshow(B); %见图6.10(d)

m = fspecial(‘average’); %生成一个空间均值滤波器

R_filtered = imfilter(R, m); %分别对图像的R、G、B分量进行滤波。

G_filtered = imfilter(G, m);

B_filtered = imfilter(B, m);

rgb_filtered = cat(3, R_filtered, G_filtered, B_filtered);

figure, imshow(rgb_filtered); %见图6.10（e）

例6.5 彩色图像锐化。

I = imread(‘’);

imshow(I);

lapMatrix = [1 1 1; 1 –8 1; 1 1 1];

I_tmp = imfilter(fb, lapMatrix, ‘replicate’);

I_sharped = imsubtract(I, I_tmp));

imshow(I_sharped);

%见图6.10（a）

% Laplacian模板

%滤波

%图像相减

%见图6.11

例6.6 空间滤波图像复原。

I ＝ imread(‘’);

imshow(I);

m = fspecial('motion', 20, 45);

%见图6.12（a）

%生成运动空间滤波器

I2 = imfilter(I, m, 'circular'); %用运动滤波器模糊原图像

noise = imnoise(zeros(size(I)), 'gaussian', 0, 0.05); %生成高斯噪声

19

figure, imshow(noise);

I3 = double(I2) + noise;

I3 = uint8(I3)

figure, imshow(I3)

I3_recovered = deconvwnr(I3, m);

figure, imshow(I3_recovered);

%见图6.12（b）

％将高斯噪声加入到图像中

%显示噪声图像如图6.12（c）

%维纳滤波复原

%显示复原图像，见图6.12（d）

例6.7 彩色图像的边缘检测。

I = imread(‘’);

sOpt = fspecial(‘sobel’); %生成sobel算子

% 计算R、G、B对x、y的偏导数

Rx = imfilter(double(I(:,:, 1)), sOpt, ‘replicate’);

Rx = imfilter(double(I(:,:, 1)), sOpt’, ‘replicate’);

Gx = imfilter(double(I(:,:, 2)), sOpt, ‘replicate’);

Gy = imfilter(double(I(:,:, 2)), sOpt’, ‘replicate’);

Bx = imfilter(double(I(:,:, 3)), sOpt, ‘replicate’);

By = imfilter(double(I(:,:, 3)), sOpt’, ‘replicate’);

% 计算点积gxx，gyy以及gxy

gxx = Rx.^2 + Gx.^2 + Bx.^2;

gyy = Ry.^2 + Gy.^2 + By.^2;

gxy = Rx.*Ry + Gx.*Gy + Bx.*By;

%计算变化率最大的角度

theta = 0.5*(atan(2*gxy./(gxx – gyy + eps)));

G1 = 0.5*((gxx + gyy) + (gxx - gyy).*cos(2*theta) + 2*gxy.*sin(2*theta));

% 由于tan函数的周期性，现旋转90度再次计算

theta = theta + pi/2;

G2 = 0.5*((gxx + gyy) + (gxx - gyy).*cos(2*theta) + 2*gxy.*sin(2*theta));

G1 = G1.^0.5;

G2 = G2.^0.5;

% 得到梯度向量

gradiant = mat2gray(max(G1, G2)); %通过上述代码得到的梯度见图6.14（b）。

例6.8 密度分层法MATLAB实现，其程序如下：

I=imread('');

imshow(I); %显示灰度图像，见图6.17(a)

G2C=grayslice(I,8); %密度分层

figure;

20

imshow (G2C, hot(8)); %显示伪彩色图像，见图6.17(b)

例6.9 灰度级-彩色变换法MATLAB实现，其程序如下：

I=imread(''); %读入灰度图像，见图6.17(a)

I=double(I);

[M,N]=size(I);

L=256;

for i=1:M %以下按图6.19的变换函数进行变换

for j=1:N

if I(i,j)

R(i,j)=0;

G(i,j)=4*I(i,j);

B(i,j)=L;

else if I(i,j)<=L/2

R(i,j)=0;

G(i,j)=L;;

B(i,j)=-4*I(i,j)+2*L;

else if I(i,j)<=3*L/4

R(i,j)=4*I(i,j)-2*L;

G(i,j)=L;

B(i,j)=0;

else

R(i,j)=L;

G(i,j)=-4*I(i,j)+4*L;

B(i,j)=0;

end

end

end

end

end

for i=1:M

for j=1:N

G2C(i,j,1)=R(i,j);

G2C(i,j,2)=G(i,j);

G2C(i,j,3)=B(i,j);

end

end

21

G2C=G2C/256;

figure;

imshow(G2C); %显示灰度级-彩色变换法结果，见图6.20

第7章 图像编码

7.2.6 行程编码(RLE）

下面给出对灰度图像和彩色图像的行程编码的MATLAB程序，并对girl图像进行行程编码。

clear all;

I=imread(''); %读入图像数据

[zipped,info]=RLEncode(I); %调用RLE进行编码

unzipped=RLEdecode(zipped,info) %调用解码程序进行解码

%显示原始图像和经编码解码后的图像,显示压缩比,并计算均方根误差得erms=0,表示%RLE是无失真编码

subplot(121);imshow(I);

subplot(122);imshow(unzipped);

erms=compare(I(:),unzipped(:))

cr=

whos I unzipped zipped

function [zipped,info]=RLEncode(vector)

[m,n]=size(vector);

vector=vector(:)';

vector=uint8(vector(:));

L=length(vector);

c=vector(1);e(1,1)=c;e(1,2)=0; %e(:,1)存放灰度,e(:,2)存放行程

t1=1

for j=1:L

if((vector(j)==c))

e(t1,2)=double(e(t1,2))+1;

else

c=vector(j);

t1=t1+1;

e(t1,1)=c;

e(t1,2)=1;

end

end

zipped=e;

=m;

=n;

[m,n]=size(e);

=m*n/(*);

function unzipped=RLEdecode(zip,info)

zip=uint8(zip);

[m,n]=size(zip);

unzipped=[];

for i=1:m

section=repmat(zip(i,1),1,double(zip(i,2)));

22

unzipped=[unzipped section];

end

unzipped=reshape(unzipped,,); %%程序结束

下面程序是对二值化后的图像用行程编码：

I=imread('');

I=im2bw(I,0.4); %图像二值化

[zipped,info]=RLEncode(I);

unzipped=RLEdecode(zipped,info);

I2=logical(unzipped)

erms=compare(I,I2)

cr=

例7.12考虑一幅大小512×512，灰度级256的标准图像Lena，用FFT实现图像数据的压缩。首先将图像分割成(512/8)2个8×8子图像，对每个子图像进行FFT，这样每个子图像有64个傅里叶变换系数。按照每个系数的方差来排序，由于图像是实值的，其64个复系数只有一半有差别的。舍去小的变换系数，就可以实现数据压缩。这里，我们保留32个系数，实现2:1的数据压缩，然后进行逆变换的程序。其Matlab程序如下：

% 设置压缩比cr

cr=0.5; % cr=0.5为2:1；cr=0.125为8:1

% 读入并显示原始图像

I1=imread(''); % 图像的大小为512×512

I1=double(I1)/255; % 图像为256级灰度图像，对图像进行归一化操作

figure(1); % 显示原始图像

imshow(I1);

% 对图像进行FFT

fftcoe=blkproc(I1,[8 8],'fft2(x)'); %将图像分割成8×8的子图像进行FFT

coevar=im2col(fftcoe,[8 8],'distinct'); %将变换系数矩阵重新排列

coe=coevar;

[y,ind]=sort(coevar);

[m,n]=size(coevar);

snum=64-64*cr; % 根据压缩比确定要变0的系数个数

% 舍去不重要的系数

for i=1:n

coe(ind(1:snum),i)=0; % 将最小的Snum个变换系数清0

end

B2=col2im(coe,[8 8],[512 512],'distinct'); % 重新排列系数矩阵

%对子图像块进行FFT逆变换获得各个子图像的恢复图像，并显示压缩图像

I2=blkproc(B2,[8 8],'ifft2(x)'); % 对截取后的变换系数进行FFT逆变换

figure(2); % 显示压缩后图像

imshow(I2);

%计算均方根误差erms

e=double(I1)-double(I2);

[m,n]=size(e);

erms=sqrt(sum(e(:).^2)/(m*n))

例7.13考虑例7.12相同的图像。首先将图像分割成(512/8)2个8×8子图像，对每个子图像进行DCT，这样每个子图像有64个变换系数，舍去50%小的变换系数，即保留32个系数，进行2:1的压缩。其Matlab程序如下：

% 设置压缩比cr

cr=0.5; % cr=0.5为2:1；cr=0.125为8:1

% 读入并显示原始图像

23

initialimage=imread('');

initialimage=double(initialimage)/255;

figure(1);

imshow(initialimage);

% 对图像进行DCT

t=dctmtx(8);

dctcoe=blkproc(initialimage,[8 8],'P1*x*P2',t,t');

coevar=im2col(dctcoe,[8 8],'distinct');

coe=coevar;

[y,ind]=sort(coevar);

[m,n]=size(coevar);

% 舍去不重要的系数

snum=64-64*cr;

for i=1:n

coe(ind(1:snum),i)=0;

end

b2=col2im(coe,[8 8],[512 512],'distinct');

% 对截取后的变换系数进行DCT逆变换

i2=blkproc(b2,[8 8],'P1*x*P2',t',t);

% 显示压缩图像

figure(2);

imshow(i2)

%计算均方根误差

e=double(initialimage)-double(i2);

[m,n]=size(e);

erms=sqrt(sum(e(:).^2)/(m*n))

例7.14仍使用例7.12中的图像。首先将图像分割成(512/8)2个8×8子图像，对每个子图像进行哈达玛变换，这样每个子图像有64个变换系数，舍去50%小的变换系数，即保留32个系数，进行2:1的压缩。其Matlab程序如下：

% 设置压缩比cr

cr=0.5; % cr=0.5为2:1；cr=0.125为8:1

% 读入并显示原始图像

I1=imread(''); % 图像的大小为512×512

I1=double(I1)/255; % 图像为256级灰度图像，对图像进行归一化操作

figure(1); % 显示原始图像

imshow(I1);

% 对图像进行哈达玛变换

T=hadamard(8); %用于产生8×8的哈达玛矩阵

htCoe=blkproc(I1,[8 8],'P1*x*P2',T,T);

CoeVar=im2col(htCoe,[8 8],'distinct');

Coe=CoeVar;

[Y,Ind]=sort(CoeVar);

[m,n]=size(CoeVar);

% 舍去不重要的系数

Snum=64-64*cr;

for i=1:n

Coe(Ind(1:Snum),i)=0;

end

B2=col2im(Coe,[8 8],[512 512],'distinct');

% 对截取后的变换系数进行哈达玛逆变换

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

% 显示压缩图像

I2=I2./(8*8);

24

figure(2);

imshow(I2);

%计算均方根误差

e=double(I1)-double(I2);

[m,n]=size(e);

erms=sqrt(sum(e(:).^2)/(m*n))

7.6.2 矢量量化过程

矢量量化分为三个主要步骤：一是训练码书；二是编码；三是解码。

1、训练码书

对图像进行矢量量化时，首先要选择码书的尺寸和码字的大小，这两个参数与图像压缩效果有直接的关系。另外，不同的码书训练方法生成的码书性能也有所不同。这里我们用最简单的LBG方法训练码书。码字为4×4的子图像块，码书尺寸为64，其Matlab程序如下：

function LBGdesign()

%读入标准图像，用于码书的训练

figure(1);

sig=imread('');

%用size函数得到图像的行数和列数

[m_sig,n_sign]=size(sig);

%设置码字的大小，4×4

siz_word=4;

%设置码书的大小

siz_book=64;

%将图像分割成4×4的子图像，作为码书训练的输入向量

num=m_sig/siz_word;

ss=siz_word*siz_word; %码字的大小

nn=num*num; %子图像个数，即输入矢量个数

re_sig=[];

for i=1:m_sig

for j=1:m_sig

f1=floor(i./siz_word);

m1=mod(i,siz_word);

if m1==0

m1=siz_word;

f1=f1-1;

end

f2=floor(j./siz_word);

m2=mod(j,siz_word);

if m2==0

m2=siz_word;

f2=f2-1;

end

re_sig(num*f1+f2+1,siz_word*(m1-1)+m2)=sig(i,j);

end

end

%码书初始化，从nn个输入矢量随机取siz_book个矢量作为初始矢量

codebook=[];

for i=1:siz_book

r=floor(rand*nn)+1;

25

codebook=[codebook;re_sig(r,:)];

end

%LBG训练算法

%d0,d1用于存放各训练矢量与其在码书中最相近的码字的距离平方之和

%sea用于存放迭代精度

d0=0;

for i=1:nn

d0=d0+vectordistance(ss,re_sig(i,:),codebook(1,:));

end

while 1

d1=0.0;

for i=1:siz_book

vectornumber(i)=0;

end

for i=1:nn

codenumber(i)=1;

min=vectordistance(ss,re_sig(i,:),codebook(1,:));

for j=2:siz_book

d=0.0;

for l=2:ss

d=d+(re_sig(i,l)-codebook(j,l))^2;

if d>=min

break;

end

end

if d

min=d;

codenumber(i)=j;

end

end

vectornumber(codenumber(i))=vectornumber(codenumber(i))+1;

d1=d1+min;

end

sea=(d0-d1)/d1;

if sea<=0.0001

break;

end

d0=d1;

for j=1:siz_book

if vectornumber[j]~=0

dd=zeros(1,ss);

for l=1:nn

if codenumber(l)==j

for k=1:ss

dd(k)=dd(k)+re_sig(l,k);

end

end

end

for k=1:ss

codebook(j,k)=dd(k)/vectornumber(j);

end

else

l=floor(rand*nn)+1;

codebook(j,:)=re_sig(l,:);

end

end

26

end

save codebook_kn codebook;

%函数vectordistance计算矢量间距离的函数

function z=vectordistance(siz_word,vector1,vector2)

z=0;

for i=1:siz_word

z=z+(vector1(i)-vector2(i))^2;

end

程序运行结果生成一个码书文件codebook_kn，是一个n行k列的矩阵。N为码书的大小，k为码书的维数，从上面程序知n=64，k=16。

2、编码

将任意一幅要压缩的图像分割成4×4的子图像块，而后按照码书中的码字索引进行编码，也就是说每一个子图像经过编码后仅用一个索引号表示，实现编码的Matlab程序如下：

function LBGencode()

%打开和显示要编码的图像

figure(1);

sig=imread('');

imagesc(sig);

colormap(gray);

axis square

axis off

[m_sig,n_sig]=size(sig);

%根据已有的码书设置分割子图像的大小和码书的大小

siz_word=4;

siz_book=64;

%调用码书

load codebook_kn

%根据码书的要求，分割要编码的图像

num=m_sig/siz_word;

ss=siz_word*siz_word;

nn=num*num;

re_sig=[];

for i=1:m_sig

for j=1:m_sig

f1=floor(i./siz_word);

m1=mod(i,siz_word);

if m1==0

m1=siz_word;

f1=f1-1;

end

f2=floor(j./siz_word);

m2=mod(j,siz_word);

if m2==0

m2=siz_word;

f2=f2-1;

end

re_sig(num*f1+f2+1,siz_word*(m1-1)+m2)=sig(i,j);

end

end

%用LBG算法编码

d1=0.0;

for i=1:nn

codenumber(i)=1;

min=vectordistance(ss,re_sig(i,:),codebook(1,:));

for j=2:siz_book

27

d=0.0;

for l=1:ss

d=d+(re_sig(i,l)-codebook(j,l))^2;

if d>=min

break;

end

end

if d

min=d;

codenumber(i)=j;

end

end

d1=d1+min;

end

3、解码

解码就是按照索引号将码书中的码字找出来，用找到的码字将图像重建出来。重建图像和原始图像之间存在一定的失真，只要其失真控制在一定的范围内，则认为该图像压缩是有效的。实现解码的Matlab程序如下：

function LBGdecode()

%解码算法

for i=1:nn

re_sig(i,:)=codebook(codenumber(i),:);

end

%重建图像

for ni=1:nn

for nj=1:ss

f1=floor(ni./num);

f2=mod(ni,num);

if f2==0

f2=num;

f1=f1-1;

end

m1=floor(nj./siz_word)+1;

m2=mod(nj,siz_word);

if m2==0

m2=siz_word;

m1=m1-1;

end

re_re_sig(siz_word*f1+m1,siz_word*(f2-1)+m2)=re_sig(ni,nj);

end

end

%显示解压后的图像，即压缩图像。

figure(2);

imagesc(re_re_sig);

colormap(gray);

axis square

axis off

附1：算术编码、解码

function y=arithmetic()

format long e;

clear all;

format long e;

symbol=['abcd']; %

ps=[0.1 0.4 0.2 0.3];

inseq=('cadacdbaaabbccccdddbbbabadba') ;

28

codeword=arithencode(symbol,ps,inseq)

outseq=arithdecode(symbol,ps,codeword,length(inseq))

%函数arithencode对symbol进行算术编码

function acode=arithencode(symbol,ps,inseq)

high_range=[];

for k=1:length(ps)

high_range=[high_range sum(ps(1:k))];

end

low_range=[0 high_range(1:length(ps)-1)];

sbidx=zeros(size(inseq));

for i=1:length(inseq)

sbidx(i)=find(symbol==inseq(i));

end

low=0;high=1;

for i=1:length(inseq)

range=high-low;

high=low+range*high_range(sbidx(i));

low=low+range*low_range(sbidx(i))

end

acode=low

%函数arithdecode对算术编码进行解码

function symbos=arithdecode(symbol,ps,codeword,symlen)

format long e;

high_range=[];

for k=1:length(ps)

high_range=[high_range sum(ps(1:k))];

end

low_range=[0 high_range(1:length(ps)-1)];

psmin=min(ps);

symbos=[];

for i=1:symlen

idx=max(find(low_range<=codeword));

codeword=codeword-low_range(idx);

if abs(codeword-ps(idx))<0.01*psmin

idx=idx+1;

codeword=0;

end

symbos=[symbos symbol(idx)];

codeword=codeword/ps(idx);

if abs(codeword)<0.01*psmin %用于退出循环

i=symlen+1;

end

end

附2：Huffman编码、解码

function HUFFMAN()

clear

load woman

%X=imread('','bmp')

data=uint8(X);

29

imshow(data)

[zipped,info]=huffencode(data); %调用Huffman编码程序进行压缩

unzipped=huffdecode(zipped,info); %调用Huffman解码程序进行解压

%显示原始图像和经编码解码后的图像,显示压缩比,并计算均方根误差得erms=0

%表示是Huffman是无编码

subplot(121);imshow(data);

subplot(122);imshow(unzipped);

erms=compare(data(:),unzipped(:))

cr=

whos data unzipped zipped

%huffencode函数对输入矩阵vector进行Huffman编码,返回编码后的向量(压缩后数据)及相关信息

function [zipped,info]=huffencode(vector)

%输入,输出都是uint8格式

%info返回解回需要的结构信息

%是添加的比特数

%des是Huffman码字

%是原始图像行数

%是原始图像列数

%是原始图像数据长度

%elen是最大码长

if ~isa(vector,'uint8')

error('input argument must be a uint8 vector');

end

[m,n]=size(vector);

vector=vector(:)';

f=frequency(vector); %计算各符号出现的概率

simbols=find(f~=0);

f=f(simbols);

[f,sortindex]=sort(f); %将符号按照出现的概率大小排列

simbols=simbols(sortindex);

len=length(simbols);

simbols_index=num2cell(1:len);

codeword_tmp=cell(len,1);

while length(f)>1 %生成Huffman树 ,得到码字得码表

index1=simbols_index{1};

index2=simbols_index{2};

codeword_tmp(index1)=addnode(codeword_tmp(index1),uint8(0));

codeword_tmp(index2)=addnode(codeword_tmp(index2),uint8(1));

f=[sum(f(1:2)) f(3:end)];

simbols_index=[{[index1,index2]} simbols_index(3:end)];

[f,sortindex]=sort(f);

simbols_index=simbols_index(sortindex);

end

codeword=cell(256,1);

codeword(simbols)=codeword_tmp;

len=0;

for index=1:length(vector) %得到整个图像所有比特数

len=len+length(codeword{double(vector(index))+1});

end

string=repmat(uint8(0),1,len);

pointer=1;

for index=1:length(vector) %对输入图像进行编码

code=codeword{double(vector(index))+1};

30

len=length(code);

string(pointer+(0:len-1))=code;

pointer=pointer+len;

end

len=length(string);

pad=8-mod(len,8); %非8整数倍 ,最后补pad

个0

if pad>0

string=[string uint8(zeros(1,pad))];

end

codeword=codeword(simbols);

codelen=zeros(size(codeword));

weights=2.^(0:23);

maxcodelen=0;

for index=1:length(codeword)

len=length(codeword{index});

if len>maxcodelen

maxcodelen=len;

end

if len>0

code=sum(weights(codeword{index}==1));

code=bitset(code,len+1);

codeword{index}=code;

codelen(index)=len;

end

end

codeword=[codeword{:}];

%计算压缩后的向量

cols=length(string)/8;

string=reshape(string,8,cols);

weights=2.^(0:7);

zipped=uint8(weights*double(string));

%码表存储到一个稀疏矩阵

huffcodes=sparse(1,1);

for index=1:nnz(codeword) %length(codeword) %nume1(codeword)

huffcodes(codeword(index),1)=simbols(index);

end

%填写解码时所需的结构信息

=pad;

des=huffcodes;

=cols./length(vector);

=length(vector);

elen=maxcodelen;

=m;

=n;

%huffdecode函数对输入矩阵vector进行Huffman解码,返回解压后的图像数据

function vector=huffdecode(zipped,info,image)

if ~isa(zipped,'uint8')

error('input argument must be a uint8 vector');

end

%产生01序列,每位占一个byte

len=length(zipped);

string=repmat(uint8(0),1,len.*8);

bitindex=1:8;

for index=1:len

31

string(bitindex+8.*(index-1))=uint8(bitget(zipped(index),bitindex));

end

string=logical(string(:)');

len=length(string);

string((+1):end)=[];

len=length(string);

%开始解码

weights=2.^(0:51);

vector=repmat(uint8(0),1,);

vectorindex=1;

codeindex=1;

code=0;

for index=1:len

code=bitset(code,codeindex,string(index));

codeindex=codeindex+1;

byte=decode(bitset(code,codeindex),info);

if byte>0

vector(vectorindex)=byte-1;

codeindex=1;

code=0;

vectorindex=vectorindex+1;

end

end

vector=reshape(vector,,);

%函数addnode添加节点

function codeword_new=addnode(codeword_old,item)

codeword_new=cell(size(codeword_old));

for index=1:length(codeword_old)

codeword_new{index}=[item codeword_old{index}];

end

%函数frequency计算各符号出现的概率

function f=frequency(vector)

if ~isa(vector,'uint8')

error('input argument must be a uint8 vector');

end

f=repmat(0,1,256);

len=length(vector);

for index=0:255

f(index+1)=sum(vector==uint8(index));

end

f=f./len;

%函数decode返回码字对应的符号

function byte=decode(code,info)

byte=des(code);

附3：JPEG 编码、解码

function d=d()

clear all

x=imread('');

figure(1);

subplot(121);

imshow(x);

y=jpegencode(x,5);

X=jpegdecode(y);

32

subplot(122);

imshow(X);

e=double(x)-double(X);

[m,n]=size(e);

erms=sqrt(sum(e(:).^2)/(m*n))

cr=imageratio(x,n)

%jpegencode函数用来压缩图像用一种近似的JPEG方法

function y=jpegencode(x,quality)

%x 为输入图像

%quality决定了截去的系数和压缩比

error(nargchk(1,2,nargin)); %检查输入参数

if nargin<2

quality=1; %缺省的quality=1

end

x=double(x)-128; %像素层次移动-128

[xm,xn]=size(x); %得到图像的尺寸

t=dctmtx(8); %得到8*8DCT变换矩阵

%将图像分割成8*8子图像,进行DCT变换,然后进行量化

y=blkproc(x,[8 8],'dct2(x)') % 'P1*x*P2',t,t');

m=[16 11 10 16 24 40 51 61 %JEPG标准化矩阵

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]*quality;

%用m标准化阵列对变换阵列进行量化,即根据式( )标准化

yy=blkproc(y,[8 8],'round(x./P1)',m);

y=im2col(yy,[8 8],'distinct'); %将图像块排列成向量

xb=size(y,2); %得到列数,也就子图像个数

order=[1 9 2 3 10 17 25 18 11 4 5 12 19 26 33 ... %变换系数排列次序

41 34 27 20 13 6 7 14 21 28 35 42 49 57 50 ...

43 36 29 22 15 8 16 23 30 37 44 51 58 59 52 ...

45 38 31 24 32 39 46 53 60 61 54 47 40 48 55 62 63 56 64];

%用Z形模板对变换系数重新排列

y=y(order,:);

eob=max(x(:))+1; %创建一个块结束符号

num=numel(y)+size(y,2);

r=zeros(num,1);

count=0;

%将非零元素重新排列放入r中.-26 -3 eob -25 1 eob

for j=1:xb %每次对一列(即一块)进行操作

i=max(find(y(:,j))); %找最后一个非零元素

if isempty(i) %没有非零元素

i=0;

end

p=count+1;

q=p+i;

r(p:q)=[y(1:i,j);eob]; %截去0 并加上结束符号

count=count+i+1;

end

33

r((count+1):end)=[]; %删除没有用的部分

r=r+128;

%保存编码信息

=uint16([xm,xn]);

cks=uint16(xb);

y=uint16(quality*100);

%对r进行huffman编码

[n ]=huffencode(uint8(r));

%jpegdecode函数jpegedcode的解码程序

function x=jpegdecode(y)

error(nargchk(1,1,nargin)); %检查输入参数

m=[16 11 10 16 24 40 51 61 %JEPG标准化矩阵

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];

order=[1 9 2 3 10 17 25 18 11 4 5 12 19 26 33 ... %变换系数排列次序

41 34 27 20 13 6 7 14 21 28 35 42 49 57 50 ...

43 36 29 22 15 8 16 23 30 37 44 51 58 59 52 ...

45 38 31 24 32 39 46 53 60 61 54 47 40 48 55 62 63 56 64];

rev=order; %计算逆序

for k=1:length(order)

rev(k)=find(order==k);

end

%ff=max(rev(:))+1;

m=double(y)/100*m;

xb=double(cks); %得到图像块数

sz=double();

xn=sz(1); %得到行数

xm=sz(2); %得到列数

x=huffdecode(n,); %huffman解码

x=double(x)-128;

eob=max(x(:)); %得到块结束符

z=zeros(64,xb);k=1;

for j=1:xb

for i=1:64

if x(k)==eob

k=k+1;break;

else

z(i,j)=x(k);

k=k+1;

end

end

end

z=z(rev,:); %恢复次序

x=col2im(z,[8 8],[xm xn],'distinct');%重新排列成图像块

x=blkproc(x,[8 8],'x.*P1',m); %反标准化

t=dctmtx(8);

x=blkproc(x,[8 8], 'P1*x*P2',t',t); %DCT反变换

x=uint8(x+128); %进行位移

34

第8章 小波图像编码

例8.2 用Matlab提供正交小波’db2’对大小为256×256的Cameraman图像作一级小波分解，其结果如图8.17所示，图8.17(a)为原图，(b)为经过一级分解后得到的图象，为了细节部分显示得清楚，图(b)中除了左上角近似图像外，其余数据均乘以4。

其Matlab程序如下：

I=imread('','tif'); %读入并显示原始图像

figure(1)；

subplot(121);imshow(I);

[ca1,ch1,cv1,cd1] = dwt2(I,'db2'); %用'db2'小波对图像进行一层小波解

I2=[ca1,ch1*4;cv1*4,cd1*4; %组成变换后的矩阵

%直接用小波系数矩阵作图像输出，imshow(I2),很多数据超范围，图像不能反映实际情况，要作一些处理。

min=min(I2(:));

max=max(I2(:));

subplot(122);imshow(I2,[min,max]); %显示变换后近似和细节图像

X=idwt2(ca1,ch1,cv1,cd1,'db2'); %用idwt2作反变换

erms=compare(I,X) %反变换结果与原始图像比较

例8.3 对大小为512×512，灰度为256的Lean用表8.4给出的小波基进行4层小波分解，然后置第一层小波系数为0，重构结果如图8.20(a)，其均方根误差erms =5.3182；接着再置第二层小波系数为0，重构结果如图8.20(b)所示，其均方根误差erms = 9.8357。

上例Matlab程序如下：

%使用表８.4给出的小波基

ld=[0 0.80976 -0.87495 ...

-0.98785 0.2668641184428723

0.2668641184428723 -0.98785 -

0.80976];

hd=[0 -0.24928

0.59470 -1.115

0.49957 -0.24928 0 0];

t=(0:9);

lr=cos(pi*(t+1)).*hd;

hr=cos(pi*t).*ld;

X=imread(''); %读入原始图像

[c,s]=wavedec2(X,4, ld,hd); %进行四层小波分解

znumb=prod(s(5,:))*3; %计算第一层小波系数个数

y=repmat(0,1,znumb);

c=[c(1:length(c)-znumb) y]; %将第一层小波系数置0

XX=waverec2(c,s,lr,hr); %进行反变换

figure(1);subplot(121);imshow(uint8(XX)) %显示没有第一层小波系数重构结果

title('无第一层小波系数重构图')

erms=compare(X,XX) %计算均方根误差

35

znumb=prod(s(5,:))*3+prod(s(4,:))*3;

y=repmat(0,1,znumb);

c=[c(1:length(c)-znumb) y];

XX=waverec2(c,s,lr,hr);

figure(1);subplot(122);imshow(uint8(XX))

title('无第一,第二层小波系数重构图')

erms=compare(X,XX)

例8.4 对大小256×256，灰度级为256的Woman图像应用全局阈值法，第一次阈值大小由函数ddencmp求得，由此压缩的图像如图8.21(b)所示；第二次阈值大小为100，由此压缩的图像如图8.21(c)所示；第三次为分层阈值法，各层各方向的阈值由函数wdcbm2计算得，由此压缩的图像如图8.21(d)所示。

%装入并显示原始图像

load woman;

nbc=size(map,1);

subplot(221);image(wcodemat(X,nbc));title('原始图像')

%用ddencmp函数求出图像的全局阈值

[thr,sorh,keepapp]=ddencmp('cmp','wv',X);

thr %显示求得的阈值

%对图像作用全局阈值

[xd,cxd,lxd,perf0,perfl2]=wdencmp('gbl',X,'bior3.5',3,thr,sorh,keepapp);

subplot(222);image(wcodemat(xd,nbc));title('全局化阈值压缩图像');

xlabel(['能量成分',num2str(perfl2),'%','零系数成分',num2str(perf0),'%']);

%用bior3.5 小波对图像进行3层分解

[c s]=wavedec2(X,3,'bior3.5');

%指定策略中的经验系数

alpha=1.5;

m=2.7*prod(s(1,:));

%用wdcbm2求出各层的各方向的阈值

[thr,nkeep1]=wdcbm2(c,s,alpha,m);

thr %显示求得的阈值求得的各层阈值

%对图像分层作阈值化

[xd1,cx1,sxd1,perf01,perfl21]=wdencmp('lvd',c,s,'bior3.5',3,thr,'s');

colormap(pink(nbc));

subplot(223);image(wcodemat(xd1,nbc));title('分层阈值化压缩图像');

xlabel(['能量成分',num2str(perfl21),'%','零系数成分',num2str(perf01),'%']);

第9章 图像分割

例9.1 Roberts，Sobel和Prewitt算子检测边缘

下面是采用Roberts、Sobel和Prewitt算子检测边缘的MATLAB程序，结果如图9.3所示。

f=imread('');

subplot(2,2,1);

imshow(f);

title('原图像');

[g,t]=edge(f,'roberts',[],'both');

subplot(2,2,2);

imshow(g);

title('Roberts算子分割结果');

36

[g,t]=edge(f,'sobel',[],'both');

subplot(2,2,3);

imshow(g);

title('Sobel算子分割结果');

[g,t]=edge(f,'prewitt',[],'both');

subplot(2,2,4);

imshow(g);

title('Prewitt算子分割结果');

例9.2 LOG算子检测边缘

下面是采用LOG算子检测边缘的MATLAB程序，结果如图9.5所示。

f=imread('');

subplot(1,2,1);

imshow(f);

title('原图像');

[g,t]=edge(f,'log');

subplot(1,2,2);

imshow(g);

title('LOG算子分割结果');

例9.3 Canny算子检测边缘

下面是采用Canny算子检测边缘的MATLAB程序，结果如图9.6所示。f=imread('');

subplot(1,2,1);

imshow(f);

title('原图像');

[g,t]=edge(f,'canny');

subplot(1,2,2);

imshow(g);

title('Canny算子分割结果');

例9.4 二值图像边界跟踪

下面是实现该算法的Matlab程序，结果如图9.8所示。

f=imread('');

subplot(2,2,1);

imshow(f);

title('原图像');

b=boundary_trace(f);

subplot(2,2,2);

imshow(b);

title('四连通边界跟踪结果');

%下面是实现边界跟踪的子函数，用于跟踪目标的外边界

function g=boundary_trace(f)

offsetr=[-1,0,1,0];

offsetc=[0,1,0,-1];

next_serach_dir_table=[4,1,2,3];

next_dir_table=[2,3,4,1]; %分别为搜索方向和顺序查找表

start=-1;

boundary=-2;

[rv,cv]=find((f(2:end-1,:)>0) & (f(1:end-2,:)==0)); %找出起始点

rv=rv+1;

startr=rv(1);

startc=cv(1);

f=im2double(f);

37

f(startr,startc)=start;

cur_p=[startr,startc];

init_departure_dir=-1;

done=0;

next_dir=2; %初始搜索方向

while ~done

dir=next_dir;

found_neighbour=0;

for i=1:length(offsetr) %四邻域方向寻找下一个边缘点

offset=[offsetr(dir),offsetc(dir)];

neighbour=cur_p+offset;

if(f(neighbour(1),neighbour(2)))~=0 %找到新的边缘点

if(f(cur_p(1),cur_p(2))==start) & (init_departure_dir==-1)

init_departure_dir=dir; %记下离开初始点时的方向

else if(f(cur_p(1),cur_p(2))==start) & (init_departure_dir==dir)

done=1;

found_neighbour=1;

break;

end

next_dir=next_serach_dir_table(dir);

found_neighbour=1; %下一个搜索方向

if f(neighbour(1),neighbour(2))~=start

f(neighbour(1),neighbour(2))=boundary;

end

cur_p=neighbour;

break;

end

dir=next_dir_table(dir);

end

end

bi=find(f==boundary);

f(:)=0;

f(bi)=1;

f(startr,startc)=1;

g=im2bw(f);

例9.5 灰度图像边界跟踪

下面是实现该算法的MATLAB程序，结果如图9.10所示。

f=imread('');

subplot(1,2,1);

imshow(f);

title('原图像');

g=im2double(f);

w=fspecial('laplacian',0); %利用laplacian微分算子计算图像梯度图

g=imfilter(g,w,'replicate');

h=boundary_trace2(g);

subplot(1,2,2); imshow(h); title('边界跟踪结果');

%下面是实现边界跟踪的子函数，用于跟踪目标的外边界

function g=boundary_trace2(f) %f为输入灰度图像，g为输出二值图像

cell_cp={[0,1],[-1,1],[-1,0],[-1,-1],[0,-1],[1,-1],[1,0],[1,1]};

cell_n={[-1,1],[1,1],[0,1];[-1,0],[0,1],[-1,1];[-1,-1],[-1,1],[-1,0];[0,-1],[-1,0],[-1,-1];[1,-1],…

[-1,-1],[0,-1];[1,0],[0,-1],[1,-1];[1,1],[1,-1], [1,0];[0,1],[1,0],[1,1]};

f=padarray(f,[1,1],0,'both');

boundaryval=-1000;

maxval=max(f(:));

38

[rv,cv]=find(f==maxval)=S=[rv(1),cv(1)];

P=S; %找出起始点

minval=min(f(:));

T=(maxval+minval)/2;

f(S(1),S(2))=boundaryval;

g=f(S(1)-1:S(1)+1,S(2)-1:S(2)+1);

maxval=max(g(:));

[rv,cv]=find(g==maxval); %在起始点的八邻域方向找出第二个边缘点

C=P+[rv(1)-2,cv(1)-2];

done=f(C(1),C(2))

f(S(1),S(2))=-boundaryval;

while ~done

f(C(1),C(2))=boundaryval;

c_p=C-P;

for dir=1:length(cell_cp)

if cell_cp{dir}==c_p

break;

end

end

maxval=boundaryval;

for i=1:size(cell_n,2) %在候选点中选出梯度值最大的像素点作为下一个边缘点

N2=C+cell_n{dir,i};

if(maxval

maxval=f(N2(1),N2(2));

N=N2;

end

end

if (f(N(1),N(2))

done=true;

break;

end

P=C;

C=N; %把前一个边缘点P设为当前点C，把新找出来的边缘点N设为当前点C

end

bi=find(f==boundaryval);

f(:)=0;

f(bi)=1;

f(S(1),S(2))=1;

f=f(2:end-1,2:end-1);

g=im2bw(f);

例9.6 Hough变换检测直线

下面是利用hough函数检测直线的MATLAB程序，结果如图9.13所示。

f=imread('');

subplot(2,2,1);

imshow(f);

title('原图像');

T=graythresh(f);

f=im2bw(f,T);

subplot(2,2,2);

imshow(f);

title('二值化图像');

f=bwmorph(f,'skel',Inf);

f=bwmorph(f,'spur',8);

subplot(2,2,3);

imshow(f);

39

title('细化图像');

[rodetect,tetadetect,Accumulator]=houghtrans(f,0.25,1,20);

subplot(2,2,4);

[m,n]=size(f);

for ln=1:length(rodetect)

if tetadetect(ln)~=0

x=0:n-1;

y=-cot(tetadetect(ln)*pi/180)*x+rodetect(ln)/sin(tetadetect(ln)*pi/180);

else

x=rodetect(ln)*ones(1,n); y=0:m-1;

end

xr=x+1;

yr=floor(y+1.0e-10)+1;

xidx=zeros(1,n);

xmin=0;

xmax=0;

for i=1:n

if(yr(i)>=1 & yr(i)<=m)

if f(yr(i),xr(i))==1

if xmin==0

xmin=i;

end

xmax=i;

end

end

end

if tetadetect(ln)~=0

x=xmin-1:xmax-1;

y=y(x+1);

else

y=xmin-1:xmax-1;

x=x(y+1);

end

y=m-1-y;

plot(x,y,'linewidth',1);

hold on;

end

axis([0,m-1,0,n-1]);

title('Hough变换检测出的直线');

例9.7 迭代式阈值选择

下面是实现该算法的MATLAB程序，结果如图9.16所示。

f=imread('');

subplot(1,2,1);

imshow(f);

title('原始图像');

f=double(f);

T=(min(f(:))+max(f(:)))/2;

done=false;

i=0;

while ~done

r1=find(f<=T);

r2=find(f>T);

Tnew=(mean(f(r1))+mean(f(r2)))/2;

done=abs(Tnew-T)<1;

T=Tnew;

40

i=i+1;

end

f(r1)=0;

f(r2)=1;

subplot(1,2,2);

imshow(f);

title('迭代阈值二值化图像');

例9.8 Otsu方法阈值选择

以下程序实现的功能是：先用graythresh函数求取阈值，然后用此阈值将灰度图像二值化，结果如图9.17所示。

f=imread('');

subplot(2,2,1);

imshow(f);

title('原始图像');

T=graythresh(f);

g=im2bw(f,T);

subplot(2,2,2);

imshow(g);

title('Otsu方法二值化图像');

例9.9分水岭算法分割图像

下面是利用分水岭算法分割图像的MATLAB程序，结果如图9.19所示。

f=imread('');

subplot(2,2,1);

imshow(f);

title('原始图像');

f=double(f);

hv=fspecial('prewitt');

hh=hv.'; %计算梯度图

gv=abs(imfilter(f,hv,'replicate'));

gh=abs(imfilter(f,hh,'replicate'));

g=sqrt(gv.^2+gh.^2); %计算距离

L=watershed(g);

wr=L==0;

subplot(2,2,2);

imshow(wr);

title('分水岭');

f(wr)=255;

subplot(2,2,3);

imshow(uint8(f));

title('分割结果');

rm=imregionalmin(g); %取出梯度图中局部极小值点

subplot(2,2,4);

imshow(rm);

title('局部极小值');

例9.10 改进的分水岭算法分割图像

下面是利用改进的分水岭算法分割图像的MATLAB程序，结果如图9.20所示。f=imread('');

subplot(2,3,1);

imshow(f);

title('原始图像');

f=double(f);

41

hv=fspecial('prewitt');

hh=hv.';

gv=abs(imfilter(f,hv,'replicate'));

gh=abs(imfilter(f,hh,'replicate'));

g=sqrt(gv.^2+gh.^2);

df=bwdist(f);

subplot(2,3,2);

imshow(uint8(df*8));

title('原图像的距离变换');

L=watershed(df); %计算外部约束

em=L==0;

subplot(2,3,3);

imshow(em);

title('标记外部约束');

im=imextendedmax(f,20); %计算内部约束

subplot(2,3,4);

imshow(im);

title('标记内部约束');

g2=imimposemin(g,im|em); %重构梯度图

subplot(2,3,5);

imshow(g2);

title('由标记内外约束重构的梯度图');

L2=watershed(g2); %调用watershed函数进行分割

wr2=L2==0;

f(wr2)=255;

subplot(2,3,6);

imshow(uint8(f));

title('分割结果');

例9.11 区域生长法分割图像

下面是利用区域生长法分割图像的MATLAB程序，结果如图9.22所示。f=imread('');

subplot(2,2,1);

imshow(f);

seedx=[30,76,86];

seedy=[110,81,110]; %选择种子点

hold on;

plot(seedx,seedy,'gs','linewidth',1);

title('原图像及种子点位置');

f=double(f);

markerim=f==f(seedy(1),seedx(1));

for i=2:length(seedx)

markerim=markerim | (f==f(seedy(i),seedx(i)));

end

thresh=[12,6,12];

maskim=zeros(size(f));

for i=1:length(seedx)

g=abs(f-f(seedy(i),seedx(i)))<=thresh(i);

maskim=maskim | g;

end

[g,nr]=bwlabel(imreconstruct(markerim,maskim),8);

g=mat2gray(g);

subplot(2,2,2);

imshow(g);

title('三个种子点区域生长结果');

42

例9.12 区域分裂合并法分割图像

下面是利用区域分裂合并法分割图像的MATLAB程序，结果如图9.24所示。

f=imread('');

[m,n]=size(f);

pow2size=2^nextpow2(max(m,n));

if m~=n | m~=pow2size

error('图像必须是方的且大小为2的整数次幂');

end

subplot(2,2,1);

imshow(f);

title('原图像');

std_thresh=10;

min_dim=2;

g=split_merge(f,min_dim,@predicate_fun,std_thresh);

g=mat2gray(g);

subplot(2,2,2);

imshow(g);

title('分裂最小子区域大小2×2');

%下面是实现区域分裂合并的子函数，利用基于四叉树的数据表示方式分裂图像function g=split_merge(f,min_dim,predicate_fun,std_thresh)

% f为输入图像，min_dim为分裂的最小子区域的大小

% predicate_fun为实现相同性质逻辑谓词P的函数，std_thresh为标准方差阈值

spare_qtim=qtdecomp(f,@split_test_fun,min_dim,@predicate_fun,std_thresh);

max_block_size=full(max(spare_qtim(:))); %取出最大块的大小

maskim=zeros(size(f));

markerim=zeros(size(f));

for i=1:max_block_size

[val,r,c]=qtgetblk(f,spare_qtim,i);

if numel(val)~=0

for j=1:length(r)

xlow=r(j);

ylow=c(j);

xhigh=xlow+i-1;

yhigh=ylow+i-1;

subblock=f(xlow:xhigh,ylow:yhigh);

lag=feval(predicate_fun,subblock,std_thresh);

if flag

maskim(xlow:xhigh,ylow:yhigh)=1;

markerim(xlow:xhigh,ylow:yhigh)=1;

end

end

end

end

g=bwlabel(imreconstruct(markerim,maskim),8);

%下面是判断是否分裂图像的子函数

function splitflag = split_test_fun(subblocks_im,min_dim,predicate_fun,std_thresh)

block_num=size(subblocks_im,3);

splitflag(1:block_num)=false; %取出子块数

for i=1:block_num

subblock=subblocks_im(:,:,i);

if(size(subblock,1))<=min_dim

splitflag(i)=false;

continue；

end

43

flag=feval(predicate_fun,subblock,std_thresh);

if flag

splitflag(i)=true;

end

end

%下面是判断图像子块一致性的子函数

function flag=predicate_fun(subblock_im,std_thresh)

stdval=std2(subblock_im);

flag=stdval>std_thresh;

例9.13 能量最小化形式分割的形变模型图像

下面是利用能量最小化形式的形变模型分割U形图像的MATLAB程序。

[I,map]=rawread(''); %读入U形图像

figure(1);

subplot(131);

imdisp(I);

title('original image');

f=1-I/255; %计算边缘图，使非边缘点灰度值小，边缘点灰度值大

subplot(132);

imdisp(f);

title('edge map');

f0=gaussianBlur(f,1); %与高斯函数做卷积运算，标准差为1

subplot(133);

imdisp(f0);

title('gaussianblurred edge map'); %去噪后的边缘图

[px,py]=gradient(f0); %计算高斯外力(梯度)

figure(2);

subplot(121);

imdisp(-f);

title('gaussian force');

subplot(122);

quiver(px,py); %用quiver函数显示箭头图

axis('square', 'equal', 'off', 'ij');

figure(2);

subplot(121);

cla;

colormap(gray(64));

image(((1-f)+1)*40);

axis('square', 'equal', 'off');

t=0:0.5:6.28;

x=32+22*cos(t);

y=32+22*sin(t);

snakedisp(x,y,'b'); %显示活动轮廓曲线

[x,y]=snakeinterp(x,y,2,0.5); %进行插值操作，2为最大间隔，0.5为最小间隔

snakedisp(x,y,'r'); %显示活动轮廓曲线

for i=1:ITER %活动轮廓曲线发生形变的迭代过程

[x,y]=snakedeform(x,y,1,0.2,10,px,py,5); %调用活动轮廓曲线发生形变子函数

[x,y]=snakeinterp(x,y,2,0.5);

snakedisp(x,y,'r');

title(['Deformation in progress,iter=' num2str(i*5)]);

pause(0.5);

end

figure(3);

colormap(gray(64));

image(((1-f)+1)*40);

44

axis equal

snakedisp(x,y,'r');

title(['Final result, iter=' num2str(40*5)]);

%下面是活动轮廓曲线发生形变的子函数

function [x,y]=snakedeform(x,y,alpha,beta,gamma,fx,fy,ITER)

% alpha为弹性系数，beta为刚性系数，gamma为权重因子，fx和fy为外力场

N=length(x);

alpha=alpha*ones(1,N);

beta=beta*ones(1,N);

% 以下产生5对角矢量

alpham1=[alpha(2:N) alpha(1)];

alphap1=[alpha(N) alpha(1:N-1)];

betam1=[beta(2:N) beta(1)];

betap1=[beta(N) beta(1:N-1)];

a=betam1;

b= -alpha-2*beta-2*betam1;

c=alpha+alphap1+betam1+4*beta+betap1;

d =-alphap1-2*beta-2*betap1;

e=betap1;

% 以下产生参数矩阵

for count=1:ITER

vfx=interp2(fx,x,y,'*linear');

vfy=interp2(fy,x,y,'*linear');

x=invAI*(gamma*x+vfx);

y=invAI*(gamma*y+vfy);

end

%下面是活动轮廓曲线显示的子函数

function snakedisp(x,y,style)

hold on

x=x(:);

y=y(:); %转为列数据

if nargin == 3

plot([x;x(1,1)],[y;y(1,1)],style,'MarkerSize',20);

hold off

else

disp('snakedisp.m: The input parameter is not correct!');

end

%下面是实现插值操作子函数

function [xi,yi]=snakeinterp(x,y,dmax,dmin)

% dmax为活动轮廓曲线上两个点之间的最大距离，dmin为最小距离

x=x(:);

y=y(:);

N=length(x);

d=abs(x([2:N 1])- x(:))+abs(y([2:N 1])- y(:)); %计算活动轮廓曲线上两个点之间的距离

IDX=(d

idx=find(IDX==0);

x=x(idx); y=y(idx);

N=length(x);

d=abs(x([2:N 1])- x(:))+abs(y([2:N 1])- y(:));

IDX=(d>dmax); %若距离大于dmax，则在两个点之间插入一个新的点

z=snakeindex(IDX); %用snakeindex函数为活动轮廓曲线上的点编号

p=1:N+1;

xi=interp1(p,[x;x(1)],z');

yi=interp1(p,[y;y(1)],z');

N=length(xi);

d=abs(xi([2:N 1])- xi(:))+abs(yi([2:N 1])- yi(:));

45

while (max(d)>dmax) %在活动轮廓曲线上增加一个新的点

IDX=(d>dmax);

z=snakeindex(IDX);

p=1:N+1;

xi=interp1(p,[xi;xi(1)],z');

yi=interp1(p,[yi;yi(1)],z');

N=length(xi);

d=abs(xi([2:N 1])- xi(:))+abs(yi([2:N 1])- yi(:));

end

例9.14 基于气球力的形变模型分割图像

下面是利用基于气球力的形变模型分割U形图像的MATLAB程序。

%除下述语句段外其余同例9.13

for i=1:ITER

[x,y]=snakedeform1(x,y,1,0.2,1,10,0.05,px,py,5); %气球力系数值设为0.05

[x,y]=snakeinterp(x,y,2,0.5);

snakedisp(x,y,'r');

title(['Deformation in progress, iter=' num2str(i*5)]);

pause(0.5);

end

%下面是活动轮廓曲线发生形变的子函数

function [x,y]=snakedeform1(x,y,alpha,beta,gamma,kappa,fx,fy,ITER)

% kappa为气球力权重因子

N=length(x);

alpha=alpha*ones(1,N);

beta=beta*ones(1,N);

%产生5对角矢量和参数矩阵的程序同例9.13

for count=1:ITER

vfx=interp2(fx,x,y,'*linear');

vfy=interp2(fy,x,y,'*linear');

xp=[x(2:N);x(1)];

yp=[y(2:N);y(1)]; %外力中增加了气球力

xm=[x(N);x(1:N-1)];

ym=[y(N);y(1:N-1)];

qx=xp-xm;

qy=yp-ym;

pmag=sqrt(qx.*qx+qy.*qy);

px=qy./pmag;

py=-qx./pmag;

x=invAI*(gamma*x+vfx+kappa.*px);

y=invAI*(gamma*y+vfy+kappa.*py);

end

例9.15 基于距离力的形变模型分割图像

下面是利用基于距离力的形变模型分割U形图像的MATLAB程序。

%除下述语句段外其余同例9.13

disp(' Compute distance force ...'); %计算由距离力产生的外力场

D=dt(f>0.5);

[px,py]=gradient(-D);

for i=1:ITER

[x,y]=snakedeform(x,y,0.05,0,0.5,px,py,5);

[x,y]=snakeinterp(x,y,2,0.5);

snakedisp(x,y,'r');

title(['Deformation in progress,iter='num2str(i*5)]);

pause(0.5);

end

46

%下面是欧氏距离变换的子函数

function D=dt(B) %B是一个二值映射图

[i,j]=find(B);

[n,m]=size(B);

for x=1:n

for y=1:m

dx=i-x;

dy=j-y;

dmag=sqrt(dx.*dx+dy.*dy);

D(x,y)=min(dmag);

end

end

例9.16 基于GVF的形变模型分割图像

下面是利用基于GVF的形变模型分割U形图像的MATLAB程序。

%除下述语句段外其余同例9.13

disp(' Compute GVF ...'); %计算梯度矢量流场

[u,v]=GVF(f,0.2,80);

mag=sqrt(u.*u+v.*v);

px=u./(mag+1e-10);

py=v./(mag+1e-10);

for i=1:ITER

[x,y]=snakedeform(x,y,0.05,0,0.5,px,py,5);

[x,y]=snakeinterp(x,y,2,0.5);

snakedisp(x,y,'r');

title(['Deformation in progress,iter=' num2str(i*5)]);

pause(0.5);

end

%下面是计算梯度矢量流的子函数

function [u,v]=GVF(f,mu,ITER)

% mu为规范化系数，ITER为迭代次数

[fx,fy]=gradient(f);

u=fx;v=fy;

SqrMagf = fx.*fx + fy.*fy; %外力场的幅度平方

for i=1:ITER

u = u + mu*4*del2(u) - SqrMagf.*(u-fx);

v = v + mu*4*del2(v) - SqrMagf.*(v-fy);

end

例9.17 基于GAC的形变模型分割图像，结果如图9.26所示。

Img=imread('');

Img=double(Img(:,:,1));

sigma=1.2;

G=fspecial('gaussian',15,sigma); % 高斯卷积做平滑处理

Img_smooth=conv2(Img,G,'same');

[Ix,Iy]=gradient(Img_smooth);

f=Ix.^2+Iy.^2;

g=1./(1+f);

alf=.4;

timestep=.1;

iterNum_evo=10;

iterNum_ri = 10;

[nrow, ncol]=size(Img); %将初始水平集函数定义为圆内带符号距离

initialLSF=sdf2circle(nrow,ncol, nrow/2,ncol/2,30);

u=initialLSF;

figure;

47

imagesc(Img);

colormap(gray);

hold on;

[c,h]=contour(u,[0 0],'g');

for n=1:ITER %水平集演化过程

u=e_GAC(u,g,alf,timestep,iterNum_evo);

u=reinit_SD(u,1,1, .5,iterNum_ri);

pause(0.05);

if mod(n,20)==0

imagesc(Img);

colormap(gray);

hold on;

[c,h]=contour(u,[0 0],'r');

iterNum=[num2str(n*iterNum_evo), 'iterations'];

title(iterNum);

hold off;

end

end

imagesc(Img);

colormap(gray);

hold on;

[c,h]=contour(u,[0 0],'r');

title('Final contour');

%下面是实现水平集演化过程的子函数

function u=e_GAC(initLSF, g, alf, delt, numIter)

%initLSF为初始化水平集函数，g为边缘检测算子，[gx, gy]为g的梯度

%alf为外力权重因子，delt为迭代时间步长，numIter为迭代次数

u=initLSF;

[gx,gy]=gradient(g);

for k=1:numIter

K=curvature(u); %求曲率

dx_f=Dx_forward(u);

dy_f=Dy_forward(u);

dx_b=Dx_backward(u);

dy_b=Dy_backward(u);

[dx_c,dy_c]=gradient(u);

norm_grad_p_u=gradientNorm_upwind(u,'+');

norm_grad_n_u=gradientNorm_upwind(u,'-');

u=u+delt*(g.*K.*sqrt(dx_c.^2+dy_c.^2)+...

+alf*(max(g,0).*norm_grad_p_u+min(g,0).*norm_grad_n_u)+...

+ (max(gx,0).*dx_b+min(gx,0).*dx_f+max(gy,0).*dy_b + min(gy,0).*dy_f));

end

例9.18 基于C-V水平集的形变模型分割图像，结果如图9.27所示。

Img=imread('');

U=Img(:,:,1);

[nrow,ncol]=size(U);

ic=nrow/2;

jc=ncol/2;

r=20;

phi_0=sdf2circle(nrow,ncol,ic,jc,r);

delta_t=0.1;

lambda_1=1;

lambda_2=1;

nu=0;

h=1;

48

epsilon=1;

mu=0.01*255*255;

I=double(U);

phi=phi_0;

figure;

imagesc(uint8(I));

colormap(gray);

hold on;

plotLevelSet(phi,0,'r'); %显示零水平集

numIter=1;

for k=1:ITER

phi=e_CV(I, phi, mu, nu, lambda_1, lambda_2, delta_t, epsilon, numIter); %水平集演化

if mod(k,10)==0

pause(.5);

imagesc(uint8(I));

colormap(gray);

hold on;

plotLevelSet(phi,0,'r');

end

end

%下面是计算圆内带符号的距离的子函数

function f=sdf2circle(nrow,ncol,ic,jc,r)

% nrow和ncol为行和列的数目，(ic,jc)为圆的中心位置，r为圆的半径

[X,Y]=meshgrid(1:ncol, 1:nrow);

f=sqrt((X-jc).^2+(Y-ic).^2)-r;

%下面是绘制零水平集的子函数

function [c,h]=plotLevelSet(u,zLevel, style)

[c,h]=contour(u,[zLevel zLevel],style); %用contour函数绘制零水平集

%下面是C-V水平集演化的子函数

function phi=e_CV(I, phi0, mu, nu, lambda_1, lambda_2, delta_t, epsilon, numIter);

%I为输入图像, phi0为水平集更新, mu为长度项权重因子, nu为面积项权重因子v

% lambda_1和lambda_2分别为目标区域能量项和背景区域能量项的权重因子o和b

% delta_t为时间步长, epsilon为计算Heaviside和Dirac函数的参数, numIter为迭代次数I=BoundMirrorExpand(I);

phi=BoundMirrorExpand(phi0);

for k=1:numIter

phi=BoundMirrorEnsure(phi);

delta_h=Delta(phi,epsilon);

Curv=curvature(phi);

[C1,C2]=binaryfit(phi,I,epsilon);

phi=phi+delta_t*delta_h.*(mu*Curv-nu-lambda_1*(I-C1).^2+lambda_2*(I-C2).^2);

end

phi=BoundMirrorShrink(phi);

例9.19

背景差值法分割图像

下面是利用背景差值法从静止背景中分割出运动目标的MATLAB程序，结果如图9.28所示。需要指出的是，将每一帧图像的灰度值减去背景图像的灰度值得到的差值图像并不完全等同于运动目标的图像，除非背景图像的像素点的灰度值全为零。

f=imread('');

subplot(2,2,1);

imshow(f);

title('原图像');

b=imread('');

subplot(2,2,2);

49

imshow(b);

title('背景图像');

df=im2double(f);

db=im2double(b);

c=df-db;

d=im2uint8(c);

subplot(2,2,3);

imshow(d);

title('差值图像');

T=50;

T=T/255;

i=find(abs(c)>=T);

c(i)=1;

i=find(abs(c)

c(i)=0;

subplot(2,2,4);

imshow(c);

title('二值化差值图像');

第10章 图像表示与描述

例10.1 计算图10.6（a）中的图像边界的链码。

图10.6（a）是一幅不规则闭合曲线的二值图像。首先求其边界。由于前面提到过的原因，为了减短码长，降低噪声干扰，在计算链码之前先对边界进行重采样。获取图像的边界并对其重采样的语句如下：

A = imread ('');

A = im2bw (A);

%载入图像

%将图像二值化。

imshow (A); %见图10.6（a）。

%获取图像的边界，见图10.6（b）

L = bwlabel(A, 8);

B = zeros(0, 2);

per = bwperim(L);

% 对图像进行标记

% 初始化B

％ 得到区域的周长的图像

％ 标记区域周长图像

％ 找到组成周长的像素的坐标

％ 初始化rr， cc

L2 = bwlabel(per,conn);

[r, c] = find(L2 == 1);

rr = zeros(length(r), 1);

cc = zeros(length(c), 1);

rr(1) = r(1); cc(1) = c(1);

r(1) = 0; c(1) = 0;

dir = 0;

％ 寻找边界

％ 方向

for j =1:1: length(r)

％ 查找当前点各个方向上的邻居点

[r1, c1] = find((r == rr(j)+1) & (c == cc(j)));

[r2, c2] = find((r == rr(j)+1) & (c == cc(j)-1));

[r3, c3] = find((r == rr(j)) & (c == cc(j)-1));

[r4, c4] = find((r == rr(j)-1) & (c == cc(j)-1));

[r5, c5] = find((r == rr(j)-1) & (c == cc(j)));

50

[r6, c6] = find((r == rr(j)-1) & (c == cc(j)+1));

[r7, c7] = find((r == rr(j)) & (c == cc(j)+1));

[r8, c8] = find((r == rr(j)+1) & (c == cc(j)+1));

％ 依次按1,8,7,6,5,4,3,2的方向把邻居点作为下一个当前点

x = 0; y = 0;

if ~isempty(r1)

x = r1; y = c1; dir = 1;

elseif ~isempty(r8)

x = r8; y = c8; dir = 8;

elseif ~isempty(r7)

x = r7; y = c7; dir = 7;

elseif ~isempty(r6)

x = r6; y = c6; dir = 6;

elseif ~isempty(r5)

x = r5; y = c5; dir = 5;

elseif ~isempty(r4)

x = r4; y = c4; dir = 4;

elseif ~isempty(r3)

x = r3; y = c3; dir = 3;

elseif ~isempty(r2)

x = r2; y = c2; dir = 2;

end ％ endif

if x == 0 & y == 0

break;

end

rr(j+1) = r(x); cc(j+1) = c(x);

r(x) = 0; c(x) = 0;

end; ％ endfor

rr(j + 1) = rr(1);

cc(j + 1) = cc(1);

B = [rr, cc];

[b1, b] = subsamp (B, 17);

％ 边界B， 见图10.6（b）

%重采样，结果见图10.6（c）

例10.3 下面是对两个不规则边界求其标记图的例子，其过程以及结果图如图10.10所示。首先读入两幅图像，然后将它们转化为二值图并且求其边界：

a1 = imread (''); ％ 载入图像

a1 = im2bw (a1);

％ 转化为二值图

％ 按例10.1的方法计算边界，见图10.10（a）

％ 载入第二幅图像

％ 转化为二值图

％ 按例10.1的方法计算边界，见图10.10（b）

b1 = a1的边界;

a2 = imread ('');

a2 = im2bw (a2);

b2 = a2的边界;

b = b1{1};

[nr, nc] = size(b);

51

接着分别计算两个边界的标记图，并把距离表示成角度的函数画出：

if isequal(b(1, :), b(nr, :)) % 若边界首尾是同一个点，则删除其中一个

b = b(1:np - 1, :);

nr = nr -1;

end

x0 = round(sum(b(:, 1))/nr);

y0 = round(sum(b(:, 2))/nr);

b(:, 1) = b(:, 1) - x0;

% 计算原点坐标

% 将边界起始点转化到上述坐标系中

b(:, 2) = b(:, 2) - y0;

% 将转换到极坐标系中

x = b(:, 2);

y = -b(:, 1);

[theta, r] = cart2pol(x, y);

theta = theta.*(180/pi);

% 将角度为负值的转化为正值

% 转换角度单位

j = theta == 0;

theta = theta.*(0.5*abs(1 + sign(theta))) - 0.5*(-1 + sign(theta)).*(360 + theta);

theta(j) = 0;

theta = round(theta);

signr = [theta, r];

[B, I, J] = unique(signr(:,1));

% 把角度和半径放在同一个矩阵中

% 删除重复的角度

% 如果最后的角度等于第一个的角度加360

% 那么删除

signr = signr(I, :);

if signr(end, 1) == signr(1, 1) + 360

signr = signr(1:end-1, :);

end

plot (signr(:, 1), signr(:, 2)) % 见图10.10（c）

例10.5 图10.15是一个对字母H计算傅立叶描述子以及通过反变换进行边界重构的过程。首先我们读入字母H的图像，并且提取其边界：

A = imread ('');

A = im2bw (A, 0.5);

imshow (A) % 见图10.15（a）

B = A的边界； % 按例10.1的方法获取边界，见图10.5（b）

% 对提取的边界计算其傅立叶描述子

[nr, nc] = size(B); ％ 边界B的大小

if nr/2 ~= round(nr/2) ％ 如果nr不是偶数

B(end + 1, :) = B(end, :);

nr = nr + 1;

end

x = 0:(nr - 1); ％ 计算系数

m = ((-1) .^ x)';

B(:, 1) = m .* B(:, 1);

B(:, 2) = m .* B(:, 2);

B = B(:, 1) + i*B(:, 2);

fd = fft(B); ％ 傅立叶变换

52

% 下面的代码是采用255项傅立叶描述子进行重构的过程：

nfd = length(fd); ％ 傅立叶描述子的长度

n = 255; ％ 用于恢复的项数

d = round((nfd - n)/2); % 计算需要清零的项数

fd(1:d) = 0; % 清零

fd(nfd - d + 1 : nfd) = 0;

iff = ifft(fd); ％ 傅立叶反变换

ifd(:, 1) = real(iff); ％ 分别提取实部和虚部

ifd(:, 2) = imag(iff);

ifd(:, 1) = m.*ifd(:, 1);

ifd(:, 2) = m.*ifd(:, 2);

将ifd转化为图像显示如图10.15（d）所示。

例10.6 图10.17（a）、（b）和（c）是三幅纹理图像，它们的归一化的直方图显示在图10.17（d）、（e）和（f）中。我们现在来计算它们的统计纹理度量，其结果列在表10.3中。

I1 = imread ('');

imshow(I1);

I2 = imread ('');

imshow(I2);

I3 = imread ('');

imshow(I3);

h = imhist(I1);

h = h/sum(h);

L = length(h);

% 见图10.17（a）

% 见图10.17（b）

% 见图10.17（c）

% I1的直方图

% 归一化，见图10.17（d）

% 直方图的长度

L = L-1;

%计算1到3阶统计矩

h = h(:);

rad = 0:L;

rad = rad./L;

m = rad * h;

rad = rad - m;

%计算统计矩

stm = zeros(1,3);

stm(1) = m;

for j = 2:n

stm(j) = (rad.^j) *h;

end

%获取非归一化的1～3阶统计矩

usm(1) = sm(1) * L;

usm(2) = sm(2) * L^2;

usm(3) = sm(3) * L^3;

%计算六个统计纹理测度

st(1) = usm(1);

st(2) = usm(2).^0.5;

st(3) = 1-1/(1+sm(2));

53

%转化为列向量

%生成随机数

%归一化

%均值

st(4) = usm(3)/(L^2);

st(5) = sum(h.^2);

st(6) = -sum(h .* log2(h + eps));

例10.7 图10.18（a）（b）给出了两幅纹理图像，一幅为鹅卵石纹理图像，另一幅为沙石纹理图像。下面的代码用于计算它们的频谱纹理特征。

I1 = imread ('');

I2 = imread ('');

％ 计算I1的频谱纹理度量

s = fftshift(fft2(I1)); ％ 傅立叶变换

s = abs(s); ％ 计算s的绝对值

[nc, nr] = size(s); ％ 计算s的大小

x0 = nc/2 + 1; ％ 计算原点坐标

y0 = nr/2 + 1;

rmax = min(nc, nr)/2 -1; %rmax是srad最大取值

srad = zeros(1, rmax);

srad(1) = s(x0, y0);

thetha = 91:270; %半圆

for r = 2:rmax ％ 对从2到rmax的半径

[x, y] = pol2cart(thetha, r); ％ 转换到极坐标

x = round(x)' + x0;

y = round(y)' + y0;

for j = 1:length(x)

srad(r) = sum(s(sub2ind(size(s), x, y)));

end

end

[x, y] = pol2cart(thetha, rmax);

x = round(x)' + x0;

y = round(y)' + y0;

sang = zeros(1, length(x));

for th = 1:length(x)

vx = abs(x(th) - x0);

vy = abs(y(th) - y0);

if ((vx == 0) & (vy == 0))

xr = x0;

yr = y0;

else

m = (y(th) - y0)/(x(th) - x0);

xr = (x0:x(th)).';

yr = round(y(th) + m * (xr - x0));

end

54

for j = 1:length(xr)

sang(th) = sum(s(sub2ind(size(s), xr, yr)));

end

end

s = mat2gray(log(1 + s));

imshow（s） ％ 见图10.18（c）

例10.8 图10.20（a）是一幅lena图，我们分别对其做旋转、镜像、尺度变换（缩小），变换后的图像见图10.20（b）、（c）以及（d）。然后计算原图及变换后图像的七个不变矩的值，对其进行比较，其结果见表10.4。从表中可见，这七个不变矩的值基本保持不变，也就是说，它们对旋转、镜像以及尺度变换不敏感。下面为相应的命令行：

I = imread ('');

I2 = imrotate (I, -4, 'bilinear'); % 逆时针旋转4度，图10.20（b）

I3 = fliplr (I);

I4 = imresize (I, 0.5, 'bilinear');

A = double(I);

% 垂直镜像，图10.20（c）

% 缩小为原图的二分之一，图10.20（d）

％ 转换为double类型

％计算七个不变矩

[nc, nr] = size(A);

[x, y] = meshgrid(1:nr, 1:nc);

x = x(:);

y = y(:);

A = A(:);

m.m00 = sum(A);

if m.m00 == 0

m.m00 = eps;

end

m.m10 = sum(x .* A);

m.m01 = sum(y .* A);

xmean = m.m10/m.m00;

ymean = m.m01/m.m00;

//计算中心矩

00 = m.m00;

02 = (sum((y-ymean).^2 .* A))/(m.m00^2);

03 = (sum((y-ymean).^3 .* A))/(m.m00^2.5);

11 = (sum((x-xmean) .* (y-ymean) .* A))/(m.m00^2);

12 = (sum((x-xmean) .* (y-ymean).^2 .* A))/(m.m00^2.5);

20 = (sum((x-xmean).^2 .* A))/(m.m00^2);

21 = (sum((x-xmean).^2 .* (y-ymean) .* A))/(m.m00^2.5);

30 = (sum((x-xmean).^3 .* A))/(m.m00^2.5);

im(1) = 20 + 02;

im(2) = (20 - 02)^2 + 4 * 11^2;

im(3) = (30 - 3 * 12)^2 + (3 * 21 - 03)^2;

55

% 得到网格

% 变为列向量

％ 计算均值

im(4) = (30 + 12)^2 + (21 + 03)^2;

im(5) = (30 - 3 * 12) * (30 + 12) …

* ((30 + 12)^2 - 3 * (21 + 03)^2)...

+ (3 * 21 - 03) * (21 + 03)...

* (3 * (30 + 12)^2 - (21 + 03)^2);

im(6) = (20 - 02)*((30 + 12)^2 - (21 + 03)^2)...

+ 4 * 11 * (30 + 12) * (21 + 03);

im(7) = (3 * 21 - 03) * (30 + 12) …

* ((30 + 12)^2 - 3 * (21 + 03)^2)...

+ (3 * 12 - 30) * (21 + 03) …

* (3 * (30 + 12)^2 - (21 + 03)^2);

例10.9 下面是一个膨胀操作的例子，其命令行如下：

I = imread (‘’); % 载入图像

imshow(I)

se = strel('ball',8, 8);

I2 = imdilate(I, se);

imshow(I2)

% 见图10.22（a）

％ 生成球形结构元素

％ 用生成的结构元素对图像进行膨胀

% 见图10.22（b）

例10.10下面的命令段给出了一个腐蚀操作的例子

I = imread(''); % 载入一幅图像

imshow(I) % 见图10.24（a）

se = strel('ball', 8, 8); ％ 生成球形结构元素

I2 = imerode(I, se); ％ 腐蚀操作

imshow(I2) % 见图10.24（b）

例10.11 图10.26（a）是一幅简单的二值图像，我们采用一个半径为5的圆作为结构元素，分别对该图像做开启和闭合操作，其开启和闭合结果分别如图10.26（b）和（c）所示：

I = imread (''); % 载入原图像

se = strel ('disk', 5, 4);

I2 = imopen (I, se);

I3 = imclose (I, se);

imshow (I);

figure, imshow (I2);

％ 生成圆形结构元素

％ 开启操作

％ 闭合操作

％ 见图10.26（a）

％ 见图10.26（b）

figure, imshow (I3); ％ 见图10.26（c）

例10.12 图10.28给出了“Hello World”图像的原图像及其边界的周长图像。其命令行在下面列出：

A = imread (‘’); ％ 载入“Hello World”图像

imshow(A) % 见图10.28（a）

A = bwperim (A); ％ 获取区域的周长

imshow (A) % 见图10.28（b）

例10.13 下面是一个细化操作的例子。图10.29（a）是‘Hello World’的原图像，图10.29（b）则是对（a）做一次细化的结果，图10.29（c）是对（a）做细化直至目标的宽度只有一个象素的结果。

56