Matlab DIP(瓦)ch11表示与描述练习

     这一章主要是练习一些图像的表示和描述，内容主要包括图像的一些常用表示方法，比如说链码，MPP，标记，边界片段，骨骼等。另一方面是图像的描述算子，比如说边界描述算子，区域描述算子等。有些课本练习的代码暂时还没有搞定。

     其练习过程和结果如下所示：

%%  数组单元概念
clc
clear
f=imread('..\images\dipum_images_ch11\Fig1102(a)(noisy_circular_stroke).tif');
imshow(f);
[a b c d]=image_stats(f);
a=image_stats1(f);
b=image_stats2(f);

%% Freeman链码及其某些变体
clc
clear
f=imread('..\images\dipum_images_ch11\Fig1102(a)(noisy_circular_stroke).tif');
subplot(231),imshow(f);
title('freeman原图');

h=fspecial('average',9);%9*9平滑算子
g=imfilter(f,h,'replicate');%滤波
subplot(232),imshow(g);
title('freeman中值滤波后图像');

g=im2bw(g,0.5);%将原图像转换为二值图像
subplot(233),imshow(g);
title('freeman图像阈值处理后');

%函数boundaries(BW,CONN,DIR)表示的是跟踪目标的边界，返回值为一个p*1的数组单元，p为目标的个数
%每一个单元又是一个Q*2的矩阵，即Q个点的x,y坐标
B=boundaries(g);

%函数cellfun(FUN,C)是指对一个单元数组中的每一个单元应用函数FUN
d=cellfun('length',B);%求B中每一个目标边界的长度,所以返回值d是一个向量

[max_d,k]=max(d);%返回向量d中最大的值，存在max_d中，k为其索引。
v=B{k(1)};%如果最大边界不止一条，则取出其中的一条即可。v是一个坐标数组
[M,N]=size(g);

%函数bound2im(b,M,N,x0,y0)是生成一副二值图像,大小为M*N，x0和y0是b中最小的x和y轴坐标
g=bound2im(v,M,N,min(v(:,1)),min(v(:,2)));

subplot(234),imshow(g);
title('取出最大边界后图像');

% [S, SU] = bsubsamp(B, GRIDSEP)指的是对边界B进行子采样，每个GRIDSEP个点采样一次
%输出S是采样后的值，输出SU是归一化的边界，什么意思？
[s,su]=bsubsamp(v,50);
g2=bound2im(s,M,N,min(s(:,1)),min(s(:,2)));
subplot(235),imshow(g2);
title('蒋采样后的边界图像');

%connectpoly(X,Y)函数是按照坐标(X,Y)顺时针或者逆时针连接成多边形
cn=connectpoly(s(:,1),s(:,2));
g2=bound2im(cn,M,N,min(cn(:,1)),min(cn(:,2)));
subplot(236),imshow(g2);
title('连接成多边形图像');
%freeman链码过程如下：

%函数fchcode为计算边界的freeman链码
c=fchcode(su);
%链码结果：

%% 获得一个区域的边界的细胞墙(木有搞定)
clc
clear
f=imread('..\images\dipum_images_ch11\Fig1107(a)(mapleleaf).tif');
subplot(231),imshow(f);
title('原图像');

B=bwperim(f,8);%返回一个包含二进制图像的边界图像，采用4近邻
subplot(232),imshow(B);
title('二进制边界图像');

Q=qtdecomp(B,0,2);%对B图进行四叉树分解，最小尺寸为2*2，用0做阈值
subplot(233),imshow(Q);
title('四叉树分解后图像')

%% 使用函数minperpoly
clc
clear
f=imread('..\images\dipum_images_ch11\Fig1107(a)(mapleleaf).tif');
subplot(231),imshow(f);
title('原图像');

b=boundaries(f,4,'cw');%找出B的边界，顺时针4邻接
b=b{1};%因为这幅图比较特殊，它只有一个边界物体
[M,N]=size(f);
xmin=min(b(:,1));
ymin=min(b(:,2));
bim=bound2im(b,M,N,xmin,ymin);%将边界物体转换成边界图像
subplot(232),imshow(bim);
title('找出的边界图像');

% [X, Y] = minperpoly(F, CELLSIZE)为计算二值图像的最小周长多边形
[x,y]=minperpoly(f,2)%2为用于形成边界方形的大小，x和y是MPP上顶点的坐标

b2=connectpoly(x,y);%将顶点依次连接成多边形
b2=bound2im(b2,M,N,xmin,ymin);
subplot(233),imshow(b2);
title('使用大小为2*2的方形边界单元获得的MPP');

[x,y]=minperpoly(f,3)%2为用于形成边界方形的大小，x和y是MPP上顶点的坐标
b3=connectpoly(x,y);%将顶点依次连接成多边形
b3=bound2im(b3,M,N,xmin,ymin);
subplot(234),imshow(b3);
title('使用大小为3*3的方形边界单元获得的MPP');

[x,y]=minperpoly(f,4)%2为用于形成边界方形的大小，x和y是MPP上顶点的坐标
b4=connectpoly(x,y);%将顶点依次连接成多边形
b4=bound2im(b4,M,N,xmin,ymin);
subplot(235),imshow(b4);
title('使用大小为4*4的方形边界单元获得的MPP');

[x,y]=minperpoly(f,8)%2为用于形成边界方形的大小，x和y是MPP上顶点的坐标
b8=connectpoly(x,y);%将顶点依次连接成多边形
b8=bound2im(b8,M,N,xmin,ymin);
subplot(236),imshow(b8);
title('使用大小为8*8的方形边界单元获得的MPP');
%使用minperpoly过程如下：

%% 标记(木有搞定)
clc
clear
f=imread('..\images\dipum_images_ch11\Fig1111(a)(boundary_sq.tif');
imshow(f);
title('标记原图');

[st,angle,x0,y0]=signature(f);
plot(angle,st)

%% 计算一个区域的骨骼
clc
clear
f=imread('..\images\dipum_images_ch11\Fig1113(a)(chromo_original).tif');
subplot(231),imshow(f);
whos f;%其实是个二值图像
title('人体染色体原始图');

f=im2double(f);%不再是二值图像了
h=fspecial('gaussian',25,25);
g=imfilter(f,h,'replicate');
subplot(232),imshow(g);
title('高斯平滑后的图像');%为什么平滑后那些杂点就没了呢，至少平滑到那点时像素值它自己的权值最大的啊？

g=im2bw(g,1.5*graythresh(g));%其中的graythresh()函数为自动阈值，处理后g为二值图像了
subplot(233),imshow(g);
title('经阈值处理后');

s=bwmorph(g,'skel',Inf);%提取骨骼,bwmorph为二进制形态学处理函数
subplot(234),imshow(s);
title('骨骼化后图像');

s1=bwmorph(s,'spur',8);%删除毛刺
subplot(235),imshow(s1);
title('8次去毛刺后');

s2=bwmorph(s1,'spur',7);%删除毛刺
subplot(236),imshow(s2);
title('再7次去毛刺后');
%提取染色体骨骼图过程如下：

%% 傅里叶描述子
clc
clear
f=imread('..\images\dipum_images_ch11\Fig1113(a)(chromo_original).tif');
subplot(331),imshow(f);
whos f;%其实是个二值图像
title('人体染色体原始图');

f=im2double(f);%不再是二值图像了
h=fspecial('gaussian',15,9);
g=imfilter(f,h,'replicate');
g=im2bw(g,0.7);%其中的graythresh()函数为自动阈值，处理后g为二值图像了
subplot(332),imshow(g);
title('经阈值处理后');

b=boundaries(g);
b=b{1};
bim=bound2im(b,344,270);%344,270为图片的尺寸
subplot(333),imshow(bim);
title('提取的边界图');

z=frdescp(b);%对边界坐标b进行傅里叶变换系数(有多少个点就有多少个系数)，将b的坐标点看成是复平面中的某个复数
z14=ifrdescp(z,546);%用50%的描述子进行逆变换
z546im=bound2im(z14,344,270);
subplot(334),imshow(z546im);
title('546个描述子恢复后');

z110=ifrdescp(z,110);%用50%的描述子进行逆变换
z110im=bound2im(z110,344,270);
subplot(335),imshow(z110im);
title('110个描述子恢复后');

z56=ifrdescp(z,56);%用50%的描述子进行逆变换
z56im=bound2im(z56,344,270);
subplot(336),imshow(z56im);
title('56个描述子恢复后');

z28=ifrdescp(z,28);%用50%的描述子进行逆变换
z28im=bound2im(z28,344,270);
subplot(337),imshow(z28im);
title('28个描述子恢复后');

z14=ifrdescp(z,14);%用50%的描述子进行逆变换
z14im=bound2im(z14,344,270);
subplot(338),imshow(z14im);
title('14个描述子恢复后');

z8=ifrdescp(z,8);%用50%的描述子进行逆变换
z8im=bound2im(z8,344,270);
subplot(339),imshow(z8im);
title('8个描述子恢复后');
%傅里叶描述子图过程：

%% 函数regionprops的使用
clc
clear
B =  [1     1     1     0     0     0     0     0
      1     1     1     0     1     1     0     0
      1     1     1     0     1     1     0     0
      1     1     1     0     0     0     0     0
      1     1     1     0     0     0     1     0
      1     1     1     0     0     0     1     0
      1     1     1     0     0     1     1     0
      1     1     1     0     0     0     0     0];
  B=bwlabel(B);%将矩阵B变成标记矩阵，默认为8连通，将连通区域用不同的整数给标记出来
  
  D=regionprops(B,'area','boundingbox')
  w=[D.Area]
  NR=length(w)
  
  V=cat(1,D.BoundingBox)%cat函数表示沿着第一个方向级联存放数组
%该函数的使用过程如下：

%% 统计纹理度量(木有搞定)
clc
clear
f=imread('..\images\dipum_images_ch11\Fig1119(a)(superconductor-with-window).tif');
sA=imhist(f);
statxture(sA)%计算图像的6个统计度量量

%% 计算频谱纹理
clc
clear

f1=imread('..\images\dipum_images_ch11\Fig1121(a)(random_matches).tif');
subplot(241),imshow(f1);
title('随机排列火柴');

f2=imread('..\images\dipum_images_ch11\Fig1121(b)(ordered_matches).tif');
subplot(245),imshow(f2);
title('规则排列火柴');

[srad1,sang1,s1]=specxture(f1);
subplot(242),imshow(s1);
title('随机排列火柴频谱图');

subplot(243),plot(srad1);
title('随机排列火柴幅度图');

subplot(244),plot(sang1);
title('随机排列火柴角度图');

[srad2,sang2,s2]=specxture(f2);
subplot(246),imshow(s2);
title('规则排列火柴频谱图');

subplot(247),plot(srad2);
title('规则排列火柴幅度图');

subplot(248),plot(sang2);
title('规则排列火柴角度图');
%计算频谱纹理图过程如下：

%% 不变矩
clc
clear
f=imread('..\images\dipum_images_ch11\Fig1123(a)(Original_Padded_to_568_by_568).tif');
imshow(f);
title('不变矩原图');
%不变矩原始图像如下：

%padarray函数为矩阵0填充，图像原始大小为400*400，但是需要旋转45度，所以最大变为568*568
fp=padarray(f,[84 84],'both');%因此四个方向每个方向都要扩充84个像素

fhs=f(1:2:end,1:2:end);%表示从第一个到最后一个，每2个采样一次
figure,subplot(221),imshow(fhs);
title('半大图像');
fhsp=padarray(fhs,[184 184],'both');

fm=fliplr(f);
subplot(222),imshow(fm);
title('镜面图像');
fmp=padarray(fm,[84 84],'both');

fr2=imrotate(f,2,'bilinear');%旋转2度，双线性插值
subplot(223),imshow(fr2);
title('旋转2度后图像');

fr45=imrotate(f,45,'bilinear');
subplot(224),imshow(fr45);
title('旋转45度后图像');
%不变矩变换图像如下：

phiorig=abs(log(invmoments(f)))
phihalf=abs(log(invmoments(fhs)))
phimirror=abs(log(invmoments(fm)))
phirot2=abs(log(invmoments(fr2)))
phirot45=abs(log(invmoments(fr45)))
%这5幅图像的7个不变矩图像如下：

%图像主分量
clc
clear
f1=imread('..\images\dipum_images_ch11\Fig1125(a)(WashingtonDC_Band1_512).tif');
subplot(231),imshow(f1);
title('可见蓝色波段图像');

f2=imread('..\images\dipum_images_ch11\Fig1125(b)(WashingtonDC_Band2_512).tif');
subplot(232),imshow(f2);
title('可见绿色波段图像');

f3=imread('..\images\dipum_images_ch11\Fig1125(c)(WashingtonDC_Band3_512).tif');
subplot(233),imshow(f3);
title('可见红色波段图像');

f4=imread('..\images\dipum_images_ch11\Fig1125(d)(WashingtonDC_Band4_512).tif');
subplot(234),imshow(f4);
title('近红外光图像');

f5=imread('..\images\dipum_images_ch11\Fig1125(e)(WashingtonDC_Band5_512).tif');
subplot(235),imshow(f5);
title('中红外光图像');

f6=imread('..\images\dipum_images_ch11\Fig1125(f)(WashingtonDC_Band6_512).tif');
subplot(236),imshow(f6);
title('热红外光图像');
%6幅多谱图像显示结果如下：

S=cat(3,f1,f2,f3,f4,f5,f6);%将6幅图像在第三维中连接起来
[X,R]=imstack2vectors(S);%将S组织到数组X中
p=princomp(X,6);%princomp用来获得6幅图像的主分量

g1=p.Y(:,1);%计算第1主维分量
g1=reshape(g1,512,512);
figure,subplot(231),imshow(g1,[]);

g2=p.Y(:,2);%计算第2主维分量
g2=reshape(g2,512,512);
subplot(232),imshow(g2,[]);

g3=p.Y(:,3);%计算第3主维分量
g3=reshape(g3,512,512);
subplot(233),imshow(g3,[]);

g4=p.Y(:,4);%计算第4主维分量
g4=reshape(g4,512,512);
subplot(234),imshow(g4,[]);

g5=p.Y(:,5);%计算第5主维分量
g5=reshape(g5,512,512);
subplot(235),imshow(g5,[]);

g6=p.Y(:,6);%计算第6主维分量
g6=reshape(g6,512,512);
subplot(236),imshow(g6,[]);
%主分量分解过程如下：

     图像的表示和描述对后面做计算机视觉处理起来很大的作用。当然针对不同的研究对象其表示和描述的特点也不同。