Matlab DIP瓦ch11表示与描述练习.docx

Matlab DIP瓦ch11表示与描述练习.docx

《Matlab DIP瓦ch11表示与描述练习.docx》由会员分享，可在线阅读，更多相关《Matlab DIP瓦ch11表示与描述练习.docx（17页珍藏版）》请在冰豆网上搜索。

MatlabDIP瓦ch11表示与描述练习

MatlabDIP（瓦）ch11表示与描述练习

    这一章主要是练习一些图像的表示和描述，内容主要包括图像的一些常用表示方法，比如说链码，MPP，标记，边界片段，骨骼等。

另一方面是图像的描述算子，比如说边界描述算子，区域描述算子等。

有些课本练习的代码暂时还没有搞定。

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

1%%数组单元概念

2clc

3clear

4f=imread（'..\images\dipum_images_ch11\Fig1102（a）（noisy_circular_stroke）.tif'）;

5imshow（f）;

6[abcd]=image_stats（f）;

7a=image_stats1（f）;

8b=image_stats2（f）;

9

10%%Freeman链码及其某些变体

11clc

12clear

13f=imread（'..\images\dipum_images_ch11\Fig1102（a）（noisy_circular_stroke）.tif'）;

14subplot（231）,imshow（f）;

15title（'freeman原图'）;

16

17h=fspecial（'average',9）;%9*9平滑算子

18g=imfilter（f,h,'replicate'）;%滤波

19subplot（232）,imshow（g）;

20title（'freeman中值滤波后图像'）;

21

22g=im2bw（g,0.5）;%将原图像转换为二值图像

23subplot（233）,imshow（g）;

24title（'freeman图像阈值处理后'）;

25

26%函数boundaries（BW,CONN,DIR）表示的是跟踪目标的边界，返回值为一个p*1的数组单元，p为目标的个数

27%每一个单元又是一个Q*2的矩阵，即Q个点的x,y坐标

28B=boundaries（g）;

29

30%函数cellfun（FUN,C）是指对一个单元数组中的每一个单元应用函数FUN

31d=cellfun（'length',B）;%求B中每一个目标边界的长度,所以返回值d是一个向量

32

33[max_d,k]=max（d）;%返回向量d中最大的值，存在max_d中，k为其索引。

34v=B{k

（1）};%如果最大边界不止一条，则取出其中的一条即可。

v是一个坐标数组

35[M,N]=size（g）;

36

37%函数bound2im（b,M,N,x0,y0）是生成一副二值图像,大小为M*N，x0和y0是b中最小的x和y轴坐标

38g=bound2im（v,M,N,min（v（:

1））,min（v（:

2）））;

39

40subplot（234）,imshow（g）;

41title（'取出最大边界后图像'）;

42

43%[S,SU]=bsubsamp（B,GRIDSEP）指的是对边界B进行子采样，每个GRIDSEP个点采样一次

44%输出S是采样后的值，输出SU是归一化的边界，什么意思？

45[s,su]=bsubsamp（v,50）;

46g2=bound2im（s,M,N,min（s（:

1））,min（s（:

2）））;

47subplot（235）,imshow（g2）;

48title（'蒋采样后的边界图像'）;

49

50%connectpoly（X,Y）函数是按照坐标（X,Y）顺时针或者逆时针连接成多边形

51cn=connectpoly（s（:

1）,s（:

2））;

52g2=bound2im（cn,M,N,min（cn（:

1））,min（cn（:

2）））;

53subplot（236）,imshow（g2）;

54title（'连接成多边形图像'）;

55%freeman链码过程如下：

56

57%函数fchcode为计算边界的freeman链码

58c=fchcode（su）;

59%链码结果：

60

61%%获得一个区域的边界的细胞墙（木有搞定）

62clc

63clear

64f=imread（'..\images\dipum_images_ch11\Fig1107（a）（mapleleaf）.tif'）;

65subplot（231）,imshow（f）;

66title（'原图像'）;

67

68B=bwperim（f,8）;%返回一个包含二进制图像的边界图像，采用4近邻

69subplot（232）,imshow（B）;

70title（'二进制边界图像'）;

71

72Q=qtdecomp（B,0,2）;%对B图进行四叉树分解，最小尺寸为2*2，用0做阈值

73subplot（233）,imshow（Q）;

74title（'四叉树分解后图像'）

75

76%%使用函数minperpoly

77clc

78clear

79f=imread（'..\images\dipum_images_ch11\Fig1107（a）（mapleleaf）.tif'）;

80subplot（231）,imshow（f）;

81title（'原图像'）;

82

83b=boundaries（f,4,'cw'）;%找出B的边界，顺时针4邻接

84b=b{1};%因为这幅图比较特殊，它只有一个边界物体

85[M,N]=size（f）;

86xmin=min（b（:

1））;

87ymin=min（b（:

2））;

88bim=bound2im（b,M,N,xmin,ymin）;%将边界物体转换成边界图像

89subplot（232）,imshow（bim）;

90title（'找出的边界图像'）;

91

92%[X,Y]=minperpoly（F,CELLSIZE）为计算二值图像的最小周长多边形

93[x,y]=minperpoly（f,2）%2为用于形成边界方形的大小，x和y是MPP上顶点的坐标

94

95b2=connectpoly（x,y）;%将顶点依次连接成多边形

96b2=bound2im（b2,M,N,xmin,ymin）;

97subplot（233）,imshow（b2）;

98title（'使用大小为2*2的方形边界单元获得的MPP'）;

99

100[x,y]=minperpoly（f,3）%2为用于形成边界方形的大小，x和y是MPP上顶点的坐标

101b3=connectpoly（x,y）;%将顶点依次连接成多边形

102b3=bound2im（b3,M,N,xmin,ymin）;

103subplot（234）,imshow（b3）;

104title（'使用大小为3*3的方形边界单元获得的MPP'）;

105

106[x,y]=minperpoly（f,4）%2为用于形成边界方形的大小，x和y是MPP上顶点的坐标

107b4=connectpoly（x,y）;%将顶点依次连接成多边形

108b4=bound2im（b4,M,N,xmin,ymin）;

109subplot（235）,imshow（b4）;

110title（'使用大小为4*4的方形边界单元获得的MPP'）;

111

112[x,y]=minperpoly（f,8）%2为用于形成边界方形的大小，x和y是MPP上顶点的坐标

113b8=connectpoly（x,y）;%将顶点依次连接成多边形

114b8=bound2im（b8,M,N,xmin,ymin）;

115subplot（236）,imshow（b8）;

116title（'使用大小为8*8的方形边界单元获得的MPP'）;

117%使用minperpoly过程如下：

118

119%%标记（木有搞定）

120clc

121clear

122f=imread（'..\images\dipum_images_ch11\Fig1111（a）（boundary_sq.tif'）;

123imshow（f）;

124title（'标记原图'）;

125

126[st,angle,x0,y0]=signature（f）;

127plot（angle,st）

128

129%%计算一个区域的骨骼

130clc

131clear

132f=imread（'..\images\dipum_images_ch11\Fig1113（a）（chromo_original）.tif'）;

133subplot（231）,imshow（f）;

134whosf;%其实是个二值图像

135title（'人体染色体原始图'）;

136

137f=im2double（f）;%不再是二值图像了

138h=fspecial（'gaussian',25,25）;

139g=imfilter（f,h,'replicate'）;

140subplot（232）,imshow（g）;

141title（'高斯平滑后的图像'）;%为什么平滑后那些杂点就没了呢，至少平滑到那点时像素值它自己的权值最大的啊？

142

143g=im2bw（g,1.5*graythresh（g））;%其中的graythresh（）函数为自动阈值，处理后g为二值图像了

144subplot（233）,imshow（g）;

145title（'经阈值处理后'）;

146

147s=bwmorph（g,'skel',Inf）;%提取骨骼,bwmorph为二进制形态学处理函数

148subplot（234）,imshow（s）;

149title（'骨骼化后图像'）;

150

151s1=bwmorph（s,'spur',8）;%删除毛刺

152subplot（235）,imshow（s1）;

153title（'8次去毛刺后'）;

154

155s2=bwmorph（s1,'spur',7）;%删除毛刺

156subplot（236）,imshow（s2）;

157title（'再7次去毛刺后'）;

158%提取染色体骨骼图过程如下：

159

160%%傅里叶描述子

161clc

162clear

163f=imread（'..\images\dipum_images_ch11\Fig1113（a）（chromo_original）.tif'）;

164subplot（331）,imshow（f）;

165whosf;%其实是个二值图像

166title（'人体染色体原始图'）;

167

168f=im2double（f）;%不再是二值图像了

169h=fspecial（'gaussian',15,9）;

170g=imfilter（f,h,'replicate'）;

171g=im2bw（g,0.7）;%其中的graythresh（）函数为自动阈值，处理后g为二值图像了

172subplot（332）,imshow（g）;

173title（'经阈值处理后'）;

174

175b=boundaries（g）;

176b=b{1};

177bim=bound2im（b,344,270）;%344,270为图片的尺寸

178subplot（333）,imshow（bim）;

179title（'提取的边界图'）;

180

181z=frdescp（b）;%对边界坐标b进行傅里叶变换系数（有多少个点就有多少个系数），将b的坐标点看成是复平面中的某个复数

182z14=ifrdescp（z,546）;%用50%的描述子进行逆变换

183z546im=bound2im（z14,344,270）;

184subplot（334）,imshow（z546im）;

185title（'546个描述子恢复后'）;

186

187z110=ifrdescp（z,110）;%用50%的描述子进行逆变换

188z110im=bound2im（z110,344,270）;

189subplot（335）,imshow（z110im）;

190title（'110个描述子恢复后'）;

191

192z56=ifrdescp（z,56）;%用50%的描述子进行逆变换

193z56im=bound2im（z56,344,270）;

194subplot（336）,imshow（z56im）;

195title（'56个描述子恢复后'）;

196

197z28=ifrdescp（z,28）;%用50%的描述子进行逆变换

198z28im=bound2im（z28,344,270）;

199subplot（337）,imshow（z28im）;

200title（'28个描述子恢复后'）;

201

202z14=ifrdescp（z,14）;%用50%的描述子进行逆变换

203z14im=bound2im（z14,344,270）;

204subplot（338）,imshow（z14im）;

205title（'14个描述子恢复后'）;

206

207z8=ifrdescp（z,8）;%用50%的描述子进行逆变换

208z8im=bound2im（z8,344,270）;

209subplot（339）,imshow（z8im）;

210title（'8个描述子恢复后'）;

211%傅里叶描述子图过程：

212

213%%函数regionprops的使用

214clc

215clear

216B=[11100000

21711101100

21811101100

21911100000

22011100010

22111100010

22211100110

22311100000];

224B=bwlabel（B）;%将矩阵B变成标记矩阵，默认为8连通，将连通区域用不同的整数给标记出来

225

226D=regionprops（B,'area','boundingbox'）

227w=[D.Area]

228NR=length（w）

229

230V=cat（1,D.BoundingBox）%cat函数表示沿着第一个方向级联存放数组

231%该函数的使用过程如下：

232

233%%统计纹理度量（木有搞定）

234clc

235clear

236f=imread（'..\images\dipum_images_ch11\Fig1119（a）（superconductor-with-window）.tif'）;

237sA=imhist（f）;

238statxture（sA）%计算图像的6个统计度量量

239

240%%计算频谱纹理

241clc

242clear

243

244f1=imread（'..\images\dipum_images_ch11\Fig1121（a）（random_matches）.tif'）;

245subplot（241）,imshow（f1）;

246title（'随机排列火柴'）;

247

248f2=imread（'..\images\dipum_images_ch11\Fig1121（b）（ordered_matches）.tif'）;

249subplot（245）,imshow（f2）;

250title（'规则排列火柴'）;

251

252[srad1,sang1,s1]=specxture（f1）;

253subplot（242）,imshow（s1）;

254title（'随机排列火柴频谱图'）;

255

256subplot（243）,plot（srad1）;

257title（'随机排列火柴幅度图'）;

258

259subplot（244）,plot（sang1）;

260title（'随机排列火柴角度图'）;

261

262[srad2,sang2,s2]=specxture（f2）;

263subplot（246）,imshow（s2）;

264title（'规则排列火柴频谱图'）;

265

266subplot（247）,plot（srad2）;

267title（'规则排列火柴幅度图'）;

268

269subplot（248）,plot（sang2）;

270title（'规则排列火柴角度图'）;

271%计算频谱纹理图过程如下：

272

273%%不变矩

274clc

275clear

276f=imread（'..\images\dipum_images_ch11\Fig1123（a）（Original_Padded_to_568_by_568）.tif'）;

277imshow（f）;

278title（'不变矩原图'）;

279%不变矩原始图像如下：

280

281%padarray函数为矩阵0填充，图像原始大小为400*400，但是需要旋转45度，所以最大变为568*568

282fp=padarray（f,[8484],'both'）;%因此四个方向每个方向都要扩充84个像素

283

284fhs=f（1:

2:

end,1:

2:

end）;%表示从第一个到最后一个，每2个采样一次

285figure,subplot（221）,imshow（fhs）;

286title（'半大图像'）;

287fhsp=padarray（fhs,[184184],'both'）;

288

289fm=fliplr（f）;

290subplot（222）,imshow（fm）;

291title（'镜面图像'）;

292fmp=padarray（fm,[8484],'both'）;

293

294fr2=imrotate（f,2,'bilinear'）;%旋转2度，双线性插值

295subplot（223）,imshow（fr2）;

296title（'旋转2度后图像'）;

297

298fr45=imrotate（f,45,'bilinear'）;

299subplot（224）,imshow（fr45）;

300title（'旋转45度后图像'）;

301%不变矩变换图像如下：

302

303phiorig=abs（log（invmoments（f）））

304phihalf=abs（log（invmoments（fhs）））

305phimirror=abs（log（invmoments（fm）））

306phirot2=abs（log（invmoments（fr2）））

307phirot45=abs（log（invmoments（fr45）））

308%这5幅图像的7个不变矩图像如下：

309

310%图像主分量

311clc

312clear

313f1=imread（'..\images\dipum_images_ch11\Fig1125（a）（WashingtonDC_Band1_512）.tif'）;

314subplot（231）,imshow（f1）;

315title（'可见蓝色波段图像'）;

316

317f2=imread（'..\images\dipum_images_ch11\Fig1125（b）（WashingtonDC_Band2_512）.tif'）;

318subplot（232）,imshow（f2）;

319title（'可见绿色波段图像'）;

320

321f3=imread（'..\images\dipum_images_ch11\Fig1125（c）（WashingtonDC_Band3_512）.tif'）;

322subplot（233）,imshow（f3）;

323title（'可见红色波段图像'）;

324

325f4=imread（'..\images\dipum_images_ch11\Fig1125（d）（WashingtonDC_Band4_512）.tif'）;

326subplot（234）,imshow（f4）;

327title（'近红外光图像'）;

328

329f5=imread（'..\images\dipum_images_ch11\Fig1125（e）（WashingtonDC_Band5_512）.tif'）;

330subplot（235）,imshow（f5）;

331title（'中红外光图像'）;

332

333f6=imread（'..\images\dipum_images_ch11\Fig1125（f）（WashingtonDC_Band6_512）.tif'）;

334subplot（236）,imshow（f6）;

335title（'热红外光图像'）;

336%6幅多谱图像显示结果如下：

337

338S=cat（3,f1,f2,f3,f4,f5,f6）;%将6幅图像在第三维中连接起来

339[X,R]=imstack2vectors（S）;%将S组织到数组X中

340p=princomp（X,6）;%princomp用来获得6幅图像的主分量

341

342g1=p.Y（:

1）;%计算第1主维分量

343g1=reshape（g1,512,512）;

344figure,subplot（231）,imshow（g1,[]）;

345

346g2=p.Y（:

2）;%计算第2主维分量

347g2=reshape（g2,512,512）;

348subplot（232）,imshow（g2,[]）;

349

350g3=p.Y（:

3）;%计算第3主维分量

351g3=reshape（g3,512,512）;

352subplot（233）,imshow（g3,[]）;

353

354g4=p.Y（:

4）;%计算第4主维分量

355g4=reshape（g4,512,512）;

356subplot（234）,imshow（g4,[]）;

357

358g5=p.Y（:

5）