十七、数字图像处理之表示与描述_图像处理之表述

（一）背景知识

（1）用于提取区域及其边界的函数

• 函数bwlabel计算二值图像中所有的连通分量（区域），调用语法为：

[L,num] = bwlabel(f,conn)
1

其中f是输入图像，conn指定了期望的连通数，num是找到的连通分量数，L是标记矩阵。

编写代码为：

f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0917(a).tif');
[L,num] = bwlabel(f,4);
1
2

（原始图像）
（找到的连通数）

• 函数bwperim可以用于检测连通性的值对区域数量的影响，调用函数为：

g = bwperim(f,conn)
1

其中，g是返回的二值图像，包含f中所有的边界像素。需要注意的是，这个函数具有特殊性，conn默认的连通数是4连通。为了获得4连通的区域边界，可以将conn指定为8；同理，为了获得8连通的区域边界，把conn设置为4即可。

• 函数bwboundaries用于提取二值图像f中所有区域的真实边界坐标，调用语法为：

B = bwboundaries(f,conn,options)
或 [B,L] = bwboundaries(f,conn,options)   %L是标记矩阵
或 [B,L,NR,A] = bwboundaries(f,conn,options)   %NR是孔洞数，A是逻辑稀疏矩阵
1
2
3

其中，f是二值图像，conn是8或4，options可以选择’holes’（提取区域和孔洞的边界、嵌套在区域内的区域边界）和’noholes’（智能提取区域或子区域的边界）。

当输出B是P×1的单位数组时，P是物体数。单位数组中的每个单元都包含一个np×2的矩阵，其中的行是边界像素的行和列的坐标，np是相应区域的边界像素数。每个边界的坐标都是以顺时针方向安排的，并且边界的最后一点与第一点相同，这样就提供了闭合的边界。

（提取区域或子区域的边界）

• 使用函数flipud改变边界B{k}的行进顺序，从顺时针改变成逆时针，调用代码：

Breversed(k) = flipud(B{k})
1

表示B{k}闭合的边界，描述父-子-孔洞的关系，调用语法：

boundaryEnclosed = find(A(:,k))
或 boundaryEnclosed = find(A(k,:))
1
2

• 函数bound2im用于构建和显示一幅包含感兴趣边界的二值图像，调用语法：

g = bound2im(b,M,N) 
1

g是一幅二值图像，产生在b中的坐标处为1、背景为0、大小尺寸为M×N。M=size(f,1) N=size(f,2)

• 函数cat的作用是，通过单位数组B的连通分量将函数bound2im用于坐标的单一数组b，调用语法为：

b = cat(1,B(:))  %1指出沿着第一维（垂直）串联
1

根据上面函数，编写实验代码：

f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0917(a).tif');
B = bwboundaries(f,'noholes');
numel(B)  %用于找出默认8连接的边界
b = cat(1,B{:});   
[M,N] = size(f);
image = bound2im(b,M,N)   %显示这些边界的二值图像
[B,L,NR,A] = bwboundaries(f) ;
numel(B)   %提取所有区域和孔洞的边界总数
numel(B)-NR  %计算空洞数
bR = cat(1,B{1:2},B{4});    %显示区域和孔洞的边界
imageBoundaries = bound2im(bR,M,N);
imageNumberedBoundaries = imageBoundaries.*L;  %显示有限的边界
1
2
3
4
5
6
7
8
9
10
11
12

（默认8连接的区域的边界）

（默认8连通提取所有区域和孔洞的边界）

（计算的空洞数）

bR = cat(1,B{:});    %显示所有的边界
imageBoundaries = bound2im(bR,M,N);
imageNumberedBoundaries = imageBoundaries.*L %显示所有的有限的边界
1
2
3

find(A(:,1))   %由B{1}得到的闭合边界数
find(A(1,:))   %由B{1}得到的边界数
A  %A中的元素
full(A)   %整个矩阵
1
2
3
4

（2）本章使用的MATLAB和IPT附加函数

• 函数imfill： 填充孔洞
  (二值函数）

gB = imfill(fB, locations, conn);
1

在输入二值图像fB的背景像素上从执行填充操作（即将背景像素值设为1），该操作从参数locations指定的点开始。当忽略conn和locations时，将会显示二值图像fB，并使用户用鼠标选择起始位置。单击鼠标左键来添加点，按backspace或Delete删除前面所选择的点，按住shift键单击或右击，或者双击可选择最后一个点，然后开始填充操作。按下Return键，可不添加任何点来结束选择。

编写实验代码：gB = imfill(fB, locations, conn); 填充二值图像中的孔洞

（灰度图像）

g = imfill(fI, conn, 'holes');  %填充灰度图像fI的孔洞
1

• 函数find：
  这个函数与bwlabel可以一起使用，返回组成某个指定物体的像素的坐标向量。

[gB,num] = bwlabel(fB);
[r, c] = find(g == 2);
1
2

• 函数sortrows：对数组进行排序

s = [1 1 1;2 2 1;8 1 2]
z = sortrows(s);  %参量s必须是矩阵或者列向量
1
2

• 函数unique：对数组排序后，去除重复行的操作

[z, m, n] = unique(S, 'rows');
1

• 函数circshift：对数组进行向上、向下和侧移这样的移位操作

z = circshift(S, [ud lr]);
1

其中，ud是S向上或向下移位的元素数目。当ud为正，移位操作为向下，否则向上；r1为正，移位操作为向右移动r1个元素，否则向左。

（3）一些基本的实用M-函数

• 函数bound2eight：从边界中移出4连通的像素，剩下8连通的像素。

b8 = bound2eight(b);
1

• 函数bound2four：从边界中移出8连通的像素，仅留下4连通的像素。

b4 = bound2four(b);
1

• 函数bsubsamp：

[s, su] = bsubsamp(b,gridsep);
1

把边界（单个的）b自取样为网络，那些网格线是由gridsep像素隔开的。输出s是具有相对b较少的像素的边界，这些点的数量由gridsep决定。输出su是已标定边界点的集合。

运用上面的函数编写实验代码：

image = bound2im(b);
b = bwboundaries(image,'noholes');  %把b变换成顺时针排序
z = connectpoly(s(:,1),s(:,2))  %使边界上的点重新连接，s(:,1)和s(:,2)分别是自取样边界的坐标
1
2
3

• 函数intline：计算连接两点间一条直线的整数坐标

[x,y] = intline(x1,x2,y1,y2)  %(x,y)是列向量，(x1,y1)和(x2,y2)是两个待连接点的整数坐标
1

（二）表示

（1）链码

概念：链码通过指定长度与方向的直线段的连接序列来表示边界。一条边界的链码取决于起点，代码可以将起点处理为方向数的循环序列和重新定义起点的方法进行归一化，因此产生最小量级的整数部分。

函数fchcode 用于计算链码，基本语法如下：

c = fchcode(b, conn, dir);  %dir指明输出链码的方向,有'same'（默认）和'reverse'两种
c = fchcode(b, conn);  %使用默认的方向
c = fchcode(b);  %使用默认的方向和连接
1
2
3

其中，输出c包含有以下字段结构：

• c.fcc = 佛雷曼链码（1×np）
• c.diff = c.fcc的一阶差分码（1×np）
• c.mm = 最小值的整数部分（1×np）
• c.diffmm = c.mm的一阶差分码（1×np）
• c.x0y0 = 链码的起点坐标（1×2）
  
  根据佛雷曼链码编写代码：

f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig1006(a).tif');
subplot(231),imshow(f),title('原始图像');
h = fspecial('average', 9);
g = imfilter(f, h ,'replicate');
subplot(232),imshow(g),title('9×9平均模板处理结果');
gB = im2bw(g,0.5);
subplot(233),imshow(gB),title('阈值处理后结果');
B = bwboundaries(gB,'noholes');
d = cellfun('length',B);
[max_d, k] = max(d);
b = B{1};
[M N] = size(g);
g = bound2im(b,M,N,min(b(:,1)),min(b(:,2)));
subplot(234),imshow(g),title('二值图像的边界图像');
[s,su] = bsubsamp(b,50);
g2 = bound2im(s,M,N);
subplot(235),imshow(g2),title('自取样的边界图像');
cn = connectpoly(s(:, 1),s(:, 2));
g3 = bound2im(cn,M,N,min(cn(:,1)),min(cn(:,2)));
subplot(236),imshow(g3),title('对上图的点进行连接');

c = fchcode(su);  %产生链码
c.x0y0
c.fcc
c.mm
c.diff
c.diffmm
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

（产生链码及其字段）

（2）使用最小周长多边形的多边形近似

数字边界能够用多边形以任意精度近似，对于闭合曲线，当多边形的顶点数目与边界点数目相同，并且每个顶点与边界点一致时，近似是精确的。
正如图中所看到的，形状的边是由4连接的直线段组成的，这种边界表示方法就是使用了最小周长多边形的多边形近似。

观察下图，内壁的凸点（白点）或外壁的凹点（黑点）对应MPP的顶点。由此引出MPP算法。
MPP算法：

1. 获取细胞联合体；Q = qtdecomp(B,threshold,[mindim maxdim])
2. 获取细胞联合体的内部区域；
3. 使用boundaries以4连接顺时针坐标序列的形式获得第2步骤的边界区域；
4. 使用函数fchcode获得该4连接序列的佛雷曼码；
5. 从链码中获得凸点（黑点）和凹点（白点）；
6. 使用黑点作为顶点构造一个初始多边形，在进一步的分析中删除位于该多边形之外的任何白点；
7. 用剩余的黑白点作为顶点构造一个多边形；
8. 删除所有的凹点的黑点；
9. 重复步骤7和8，直至变换停止，此时所有的180°的顶点都删除，留下的只有MPP的顶点。

实现MPP算法中使用到的一些M-函数：

1. 函数qtdecomp：将四叉树进行分解的函数，它是将一幅图像分解成等大小的平方块，再进行测试处理。

Q = qtdecomp(B,threshold,[mindim maxdim])
1

其中，Q是包含四叉树结构的稀疏矩阵。如果Q(k,m)非零，那么(k,m)为分解块的左上角，块的大小是Q(k,m)。

2. 函数qtgetblk：获取四叉树分解中块的值。

[vals,r,c] = qtgetblk(B,Q,mindim)
1

其中，vals是在B的四叉树分解中包含mindim×maxdim块值的数组，而Q是由函数qtdecomp返回的稀疏矩阵。参数r和c是包含块的左上角行和列坐标的向量。

编写实验代码：

f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0216(a).tif');
g = bwperim(f,8);
Q = qtdecomp(g,0,2);
gF = qtgetblk(g,Q,2);
R = imfill(gF,'holes') & g ;
B = bwboundaries(R,4,'noholes');
b = B(1);
1
2
3
4
5
6
7

3. 函数inpolygon用于决定点是否在多边形的边界外或边界内：

inpolygon IN = inpolygon(X,Y,xv,yv)
1

4. 函数im2minperpoly用于实现MPP算法：

[X,Y,R] = im2minperpoly(f,cellsize) 
1

其中，f是一幅包含单一区域或边界的二值输入图像，cellsize是指定用于围绕边界的细胞组合体中方形单元的大小。X和Y包含MPP顶点的x和y坐标。

编写代码：

f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig1107(a).tif');
subplot(321),imshow(f),title('原始图像');
B = bwboundaries(f,4,'noholes');
b = B{1};
[M,N] = size(f);
bim = bound2im(b,M,N);
subplot(322),imshow(bim),title('4连接边界的图像')
[X,Y] = im2minperpoly(f,2);
b2 = connectpoly(X,Y);
bCcellsize2 = bound2im(b2,M,N);
subplot(323),imshow(bCcellsize2),title('使用方形且大小为2的单元得到的MPP');
[X,Y] = im2minperpoly(f,3);
b3 = connectpoly(X,Y);
bCcellsize3 = bound2im(b3,M,N);
subplot(324),imshow(bCcellsize3),title('使用方形且大小为3的单元得到的MPP');
[X,Y] = im2minperpoly(f,16);
b16 = connectpoly(X,Y);
bCcellsize16 = bound2im(b16,M,N);
subplot(325),imshow(bCcellsize16),title('使用方形且大小为16的单元得到的MPP');
[X,Y] = im2minperpoly(f,32);
b32 = connectpoly(X,Y);
bCcellsize32 = bound2im(b32,M,N);
subplot(326),imshow(bCcellsize32),title('使用方形且大小为32的单元得到的MPP');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

实验分析： 实验结果显示了使用大小为2、3、16和32的方形单元得到的MPP，由于使用大于2*2的单元，导致了分辨率会降低，叶茎的细节丢失了，随着使用的方形的单元的大小增大，发现结果只能显示基本轮廓，叶子的细节基本都丢失了。

（3）标记

概念： 标记(signature)，图像处理术语。标记是边界的一维表达，可以用多种方法来产生。其基本思想是将原始的二维边界用元函数来表示，以降低表达难度。最简单的方法就是把从重心到边界的距离作为角度的函数来标记。标记与链码同样不受边界平移的影响，但是当边界旋转或者发生尺度变换时，标记格会发生改变。（来自百度百科）

函数signature用于寻找边界的标记，基本语法是：

[dist,angle] = signature(b,x0,y0)  %（x0,y0）是该点的坐标，度量该点到边界的距离，默认使用边界的质心坐标
1

函数cart2pol用于把笛卡尔坐标转换成极坐标：[THETA,RHO] = cart2pol(X,Y)

函数pol2cart用于反变换以返回到笛卡尔坐标：[X,Y] = pol2cart(THETA,RHO)

编写实验代码：

fsq= imread('C:\Users\Public\Pictures\Sample Pictures\Fig1111(a).tif');
subplot(221),imshow(fsq),title('非规则方形图像');
ftr= imread('C:\Users\Public\Pictures\Sample Pictures\Fig1111(b).tif');
subplot(222),imshow(ftr),title('非规则三角形图像');
bSq = bwboundaries(fsq,'noholes');
[distSq,angleSq] = signature(bSq{1});
subplot(223),plot(angleSq,distSq),title('非规则方形的标记');
btr = bwboundaries(ftr,'noholes');
[disttr,angletr] = signature(btr{1});
subplot(224),plot(angletr,disttr),title('非规则三角形的标记');
1
2
3
4
5
6
7
8
9
10

（4）边界片段

说明： 将边界分解为片段降低了边界的复杂度，还简化了描述过程。当边界包含一个或多个携带形状信息的重要凹面时，这种边界片段的分解方法是比较有用的。集合S的凸壳H是包含S的最小凸集，集合的差H-S被称为S的凸缺D。

（5）骨骼

说明：将平面区域结构形状简化为一幅图形，称为细化（或者骨骼化）。函数bwmorph用于产生二值图像B中所有区域的骨骼，调用语法为：

skeletonImage = bwmorph(B,'skel',Inf)
1

编写实验代码为：

f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig1113(a).tif');
subplot(231),imshow(f),title('原始的图像');
f=im2double(f);   %不是二值图像
h = fspecial('gaussian',25,15);     %加入高斯噪声
g = imfilter(f,h,'replicate');          %用高斯空域模板平滑后的图像
subplot(232),imshow(g),title('用高斯空域模板平滑后的图像');
g = im2bw(g,1.5*graythresh(g));    %对平滑后图像进行阈值处理
subplot(233),imshow(g),title('平滑后图像进行阈值处理');
s = bwmorph(g,'skel',Inf);   
subplot(234),imshow(s),title('骨骼处理后得到的图像');
s1 = bwmorph(s,'spur',8);
subplot(235),imshow(s1),title('刺状突起去除8次后骨骼化得到的图像');
s2 = bwmorph(s,'spur',7);
subplot(236),imshow(s2),title('刺状突起去除7次后骨骼化得到的图像');
1
2
3
4
5
6
7
8
9
10
11
12
13
14

（三）边界描述子

（1）一些简单的描述子

4连接边界长度被定义为边界像素的数目减1；
8连接则是以1计算垂直和水平过渡长度，以根号二计算对角过渡长度。

可能会用到的M-函数：

1. 函数bwperim用于提取包含在图像f中的物体的边界：g = bwperim(f,conn)
2. 函数diameter用于计算边界的直径、长轴、短轴和边界或区域的基本矩形：s = diameter(L)
   其中，L是标记矩阵，s是具有下面字段的结构：
   s.Diameter ：标量，边界或区域中任意两个像素之间的最大距离。
   s.MajorAxis：22矩阵，其中的行包含边界或区域的长轴端点的行列坐标。
   s.MimorAxis：22矩阵，其中的行包含边界或区域的短轴端点的行列坐标。
   s.BasicRectangle ：4*2矩阵，其中的每行包含基本矩阵的某个拐角的行列坐标。

（2）形状数

边界的形状数一般是以4方向的佛雷曼链码为基础的，形状数被定义为最小幅值的一阶差分，形状数的阶定义为表示形状数的数字的个数。 边界的形状数可以用函数fchcode中的c.diffmm给出：length(c.diffmm)。

（3）傅里叶描述子

傅里叶描述子的基本思想：假定物体的形状是一条封闭的曲线，沿边界曲线上的一个动点P(l)的坐标变化x(l)+iy(l)是一个以形状边界周长为周期的函数，这个周期函数可以用傅里叶级数展开表示，傅里叶级数中的一系列系数z(k)是直接与边界曲线的形状有关的，称为傅里叶描述子。使用函数frdescp用于计算边界S的傅里叶描述子，而函数ifrelescp用于做描述子的逆变换。

编写实验代码：

f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig1113(a).tif');
subplot(331),imshow(f),title('原始图像');
f=im2double(f);
h=fspecial('gaussian',15,9);
g=imfilter(f,h,'replicate');
g=im2bw(g,0.7);    %变成二值图像
subplot(332),imshow(g),title('阈值处理后');
b=boundaries(g);
b=b{1};
bim=bound2im(b,688,540);    %688×540是图像的尺寸
subplot(333),imshow(bim),title('提取的边界图');
z=frdescp(b);%对边界坐标b进行傅里叶变换系数(有多少个点就有多少个系数)，将b的坐标点看成是复平面中的某个复数
z546=ifrdescp(z,546);%用50%的描述子进行逆变换
z546im=bound2im(z546,688,540);
subplot(334),imshow(z546im),title('546个描述子恢复后');
z110=ifrdescp(z,110);%用10%的描述子进行逆变换
z110im=bound2im(z110,688,540);
subplot(335),imshow(z110im),title('110个描述子恢复后'); 
z56=ifrdescp(z,56);%用5%的描述子进行逆变换
z56im=bound2im(z56,688,540);
subplot(336),imshow(z56im),title('56个描述子恢复后');
z28=ifrdescp(z,28);%用2.5%的描述子进行逆变换
z28im=bound2im(z28,688,540);
subplot(337),imshow(z28im),title('28个描述子恢复后');
z14=ifrdescp(z,14);%用1.25%的描述子进行逆变换
z14im=bound2im(z14,688,540);
subplot(338),imshow(z14im),title('14个描述子恢复后');
z8=ifrdescp(z,8);%用0.70%的描述子进行逆变换
z8im=bound2im(z8,688,540);
subplot(339),imshow(z8im),title('8个描述子恢复后');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

实验分析： 实验将图像的边缘提取了之后，发现图像的一些细微细节已经丢失，只留下了图像的轮廓，分别使用不同的傅里叶描述子，随着描述子的减少会使得图像重构只留下了主要特征，细节特征基本丢失。当描述子减少到8个时，发现产生的图像不是原来的形状了，越来越少的描述子会使得图像变成椭圆，最后变成圆。

（4）统计矩

在数学和统计学中，矩（moment） 是对变量分布和形态特点的一组度量。n阶矩 被定义为一变量的n次方与其概率密度函数（Probability Density Function, PDF）之积的积分。在文献中n阶矩通常用符号μn表示，直接使用变量计算的矩被称为原始矩（raw moment），移除均值后计算的矩被称为中心矩（central moment）。变量的一阶原始矩等价于数学期望（expectation）、二至四阶中心矩被定义为方差（variance）、偏度（skewness）和峰度（kurtosis）。（来自百度百科）

一维边界形状的表示（如边界线段和信号波段）可以使用统计矩定量地进行描述，如均值、方差和高阶矩。可以使用函数statmoments来计算统计矩。

当g( r )是上面描绘的一维函数线段时，r是随机变量，则矩为：

其中，

（5）拐角

基本上之前所说的边界描述子都是全局性的，而拐角检测的方法是局部边角描述子，有两种方法常用在图像处理上，分别是：Harris-Stephens拐角检测器和最小特征值拐角检测器。

• Harris-Stephens拐角检测器

说明： Harris-Stephens拐角检测是一种简单的点特征提取 的算子，这种算子受信号处理中自相关函数的启发，给出与自相关函数相联系的矩阵M，M矩阵的特征值是自相关函数的一阶曲率，如果两个曲率值都高，那么就认为这点是特征点。 这个检测器是对Moarvec提出的基本技术的一种改进，Moarvec方法考虑一幅图像中的局部窗口，并且确定图像灰度的平均变化，这些变化是由在各个方向上小量的窗口移动产生的。

假设令w(x,y)表示空间平均（光滑）模板，其中所有的元素都是非负的，在图像f(x,y)的任何坐标(x,y)处，灰度上的平均变化E(x,y)可以定义为：

其中，(s,t)的值是，w和图像区域相当于方括号中表示的重叠部分。

对于小部分的移动，可以仅仅用线性项近似，表示为：

Harris-Stephens拐角检测器则可以近似为空间滤波模板[-1 0 1]^T和[-1 0 1]的偏微分，表示为：

也可以写成矩阵形式：

其中，
这个矩阵的元素是由平均模板w扫过子图像区域垂直和水平滤波的微分。

Harris-Stephens拐角检测器总结为：
1）对图像进行高斯滤波；
2）对每个像素，估计其垂直方向的梯度大小值，使用近似于导数的核做两次一维卷积即可；
3）对每一像素核给定的邻域窗口：计算局部结构矩阵A和计算响应函数R(A)；
4）选取响应函数R(A)的一个阈值，以选择最佳候选角点，并完成非最大化抑制。

• 最小特征值拐角检测器

说明： 该检测器就是利用**最小特征值算法(Minimum eigenvalue algorithm)**检测角点的一种方法。

调用语法：

points = detectMinEigenFeatures(I);
或者 points = detectMinEigenFeatures(I,Name,Value);
1
2

其中，I为灰度图像，points为返回的Harris角点检测算法检测到的特征点，Name必须为用单引号对包含的如下字符串名称，Vaule为对应Name的值。

上面两种方法在matlab图像处理当中使用函数cornermetric来实现角点检测，语法为：

C = cornermetric(f,method,param1,val1,param2,val2)
1

其中，f是输入图像，method是’Harris’或者’MinimumEigenvalue’，param1是’FilterCoefficients’，val1是包含一维空间滤波模板系数的向量来产生二维平方滤波器w，param2是’SenstiveityFactor’（仅适用于Harris），val2是敏感因子k的值（0<k<0.25，默认为0.04）。

把通过阈值测试的点作为拐角点，用函数cornerprocess来检测这些点，调用语法为：CP = cornerprocess(C,T,q)，其中C是cornermetric的输出，T是指定的阈值，q是用于减少拐角点的方形的形态学结构元的大小。

编写寻找拐角点的实验代码为：

f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig1119(a).tif');
subplot(331),imshow(f),title('原始图像');
CH = cornermetric(f,'Harris');
CH(CH < 0) = 0;
CH = mat2gray(CH);
subplot(332),imshow(imcomplement(CH),[]),title('Harris的原始输出结果');
CM = cornermetric(f,'MinimumEigenvalue');
CM = mat2gray(CM);
subplot(333),imshow(imcomplement(CH)),title('最小特征值检测器的输出结果');
hH = imhist(CH);
subplot(334),imhist(CH),title('Harris的原始输出直方图');
hM = imhist(CM);
subplot(335),imhist(CM),title('最小特征值检测器的输出直方图');
TH = percentile2i(hH,0.9945);
TM = percentile2i(hM,0.9970);

cpH = cornerprocess(CH,TH,1);
subplot(336),imshow(cpH),title('Harris的原始输出使用函数后的输出');
cpM = cornerprocess(CM,TM,1);
subplot(337),imshow(cpM),title('最小特征值原始输出使用函数后的输出');

[xH yH] = find(cpH);
subplot(338),imshow(f),title('Harris：包含在原始图像上的拐角点');
hold on
plot(yH(:)',xH(:)','wo')

[xM yM] = find(cpM);
subplot(339),imshow(f),title('最小特征值：包含在原始图像上的拐角点');
hold on
plot(yM(:)',xM(:)','wo')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

（四）区域描述子

（1）函数regionprops

这个函数是用于计算区域描绘子的主要工具，调用语法为：

D = regionprops(L,properties)  %L是标记矩阵，D是长度为max(L(:))的结构，properties是字符串
1

其中，变量properties包含字符串的单元数组，单个字符串’all’或者字符串’basic’。如果properties是字符串’all’，则计算所有描述子；若不指定或者是字符串’basic’，那么只计算’Area’（区域中的像素数）、‘Centroid’（区域的质心）、‘BoundingBox’（包含区域的最小矩阵）。

调用函数，编写实验代码：

f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig1107(a).tif');
B = bwlabel(f);
D = regionprops(B,'area','boundingbox')  %得到图像中每个区域的面积和边界限制盒
A = [D.Area]   %提取区域的面积
NR = numel(A)   %提取区域的个数
V = cat(1,D.BoundingBox)  %得到每个区域的边界限制盒
1
2
3
4
5
6

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];
D = regionprops(B,'area','boundingbox')  %得到图像中每个区域的面积和边界限制盒
A = [D.Area]   %提取区域的面积
NR = numel(A)   %提取区域的个数
V = cat(1,D.BoundingBox)  %得到每个区域的边界限制盒
1
2
3
4
5
6
7
8
9
10
11
12

（2）纹理

描述区域的一种重要方法就是量化区域的纹理内容，是基于统计方法和谱测度方法计算纹理。

• 统计法
  这种分析方法是基于灰度直方图的统计特性。自定义函数statxture用于计算纹理测量，语法为：

t = statxture(f,scale) 
1

其中，f是输入图像，t是含有6个元素的行向量，scale也是含有6个元素的行向量。

有以下六种基于灰度直方图的纹理描述子：
编写实验代码：

f1 = imread('C:\Users\Public\Pictures\Sample Pictures\Fig1121(a).tif');
f2 = imread('C:\Users\Public\Pictures\Sample Pictures\Fig1121(b).tif');
f3 = imread('C:\Users\Public\Pictures\Sample Pictures\Fig1121(c).tif');
subplot(131),imshow(f1),title('原始图像一');
subplot(132),imshow(f2),title('原始图像二');
subplot(133),imshow(f3),title('原始图像三');
t1 = statxture(f1,6);  %统计纹理的度量
t2 = statxture(f2,6);  %统计纹理的度量
t3= statxture(f3,6);  %统计纹理的度量
subplot(131),imhist(f1),title('原始图像一直方图');
subplot(132),imhist(f2),title('原始图像二直方图');
subplot(133),imhist(f3),title('原始图像三直方图');
1
2
3
4
5
6
7
8
9
10
11
12

仅仅用直方图计算的纹理度量并不携带关于像素彼此相对位置的信息，而这些信息在描述纹理时很重要，因此引入一种描述方法叫做灰度共生矩阵，这个矩阵不仅把它看成灰度分布，而且还有图像中像素的相对位置。产生共生矩阵示意图如下所示：
在图像处理工具箱中，使用函数graycomatrix计算共生矩阵，语法为：

[GS,FS] = graycomatrix(f,'NumLevels',n,'Offset',offsets)
1

其中，f是任何有效类的图像，FS是结果图像，使用函数产生GS。

编写实验代码：

f = [1 1 7 5 3 2;5 1 6 1 2 5;8 8 6 8 1 2;4 3 4 5 5 1;8 7 8 7 6 2;7 8 6 2 6 2];
f = mat2gray(f);
offsets = [0 1];
[GS,IS] = graycomatrix(f,'NumLevels',8,'Offset',offsets)    %计算共生矩阵
1
2
3
4

工具箱使用函数graycoprops产生描述子：

stats = graycoprops(GS,properties);
1

基于共生矩阵描述子的实验编写代码：

f1 = imread('C:\Users\Public\Pictures\Sample Pictures\Fig1124(a).tif');
f2 = imread('C:\Users\Public\Pictures\Sample Pictures\Fig1124(b).tif');
f3 = imread('C:\Users\Public\Pictures\Sample Pictures\Fig1124(c).tif');
subplot(311),imshow(f1),title('第一幅图像');
subplot(312),imshow(f2),title('第二幅图像');
subplot(313),imshow(f3),title('第三幅图像');
G2 = graycomatrix(f2,'NumLevels',256);   %对第二幅图进行计算其描述子
G2n = G2/sum(G2(:));
stats2 = graycoprops(G2,'all');
maxProbability2 = max(G2n(:))   %最大概率
contrast2 = stats2.Contrast  %相邻像素间灰度对比度的度量
corr2 = stats2.Correlation    %相关性
energy2 = stats2.Energy   % 能量
hom2 = stats2.Homogeneity   %同质
for I = 1:size(G2n,1);
sumcols(I) = sum(-G2n(I,1:end).*log2(G2n(I,1:end)+eps));
end
entropy2 = sum(sumcols)   %熵
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

（第一幅图的描述子参数）

（第二幅图的描述子参数）

（第三幅图的描述子参数）

实验分析： 以上是三幅图（随机的、周期的、混合纹理）的图像及其描述子的参数，看出了混合纹理图像的共生矩阵的最大概率是最大的，这是因为在混合纹理图像中，水平方向有以低灰度变换为特征的大的区域，所以其大概是最大的。相关性描述子 表明，第二幅图中的灰度是高度相关的。对比度描述子参数表明，对于非随机图像对比度比较低，随机性越高的图像的对比度反而最高。能量参数说明，周期性的图像能量是最高的，能量是以概率平方值为函数而增加的。同质性参数 用于度量关于主对角线G值的集中度，越接近对角线的最高概率值的矩阵就有高的同质值，这样的矩阵相当于具有丰富的灰度级内容和缓慢变化灰度值的区域。熵描述子在周期性图像中比较低，在随机和混合图像中，熵值比较高。

offsets = [zeros(50,1) (1:50)'];
G1 = graycomatrix(f1,'Offset',offsets);
stats1 = graycoprops(G1,'Correlation');
subplot(131),plot([stats1.Correlation]);
xlabel('Horizontal Offset'),ylabel('Correlation')
G2 = graycomatrix(f2,'Offset',offsets);
stats2 = graycoprops(G2,'Correlation');
subplot(132),plot([stats2.Correlation]);
xlabel('Horizontal Offset'),ylabel('Correlation')
G3 = graycomatrix(f3,'Offset',offsets);
stats3 = graycoprops(G3,'Correlation');
subplot(133),plot([stats3.Correlation]);
xlabel('Horizontal Offset'),ylabel('Correlation')
1
2
3
4
5
6
7
8
9
10
11
12
13

• 纹理的频谱度量
  纹理的频谱度量基于傅里叶频谱，频谱对于描述一幅图像中周期的或近似周期的二维模式的方向性是合适的。

用极坐标表示频谱得到的函数S(r,θ)可以简化对频谱特性的解释，S是频谱函数，r和θ是极坐标系统中自变量。对于每一个方向θ，S(r,θ)是一维函数，可以写为；类似地，对于每一个频率r，S(r,θ)可以表示为。将θ视为固定值分析，就可以得到沿着原点的径向频谱的状况（如峰值），同样的，将r视为固定值分析，就可以得到沿着以原点为中心的圆的频谱特性。

全局描述可以通过对这些函数的积分获得：

其中，R0是以原点为中心的圆的半径。

在图像处理工具箱中，可以使用函数specxture进行计算两个纹理度量，基本语法为：

[srad,sang,S] = specxture(f);   %srad是S(r),sang是S(θ),S是频谱图
1

编写实验代码：

f1 = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0807(a).tif');
subplot(221),imshow(f1),title('随机火柴的图像');
[srad1,sang1,S1] = specxture(f1);  
subplot(222),plot(srad1),title('S(r)图像');
subplot(223),plot(sang1),title('S(θ)图像');
subplot(224),imshow(S1),title('随机火柴的频谱图像');
axis([200 200 200 200])

f2 = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0807(c).tif');
subplot(221),imshow(f2),title('整齐排列的火柴图像');
[srad2,sang2,S2] = specxture(f2);  
subplot(222),plot(srad2),title('S(r)图像');
subplot(223),plot(sang2),title('S(θ)图像');
subplot(224),imshow(S2),title('整齐火柴的频谱图像');
axis([200 200 200 200])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

实验分析： 对于随机摆放和整齐摆放的物体图像，使用傅里叶频谱可以显示不同的图像。对于随机摆放的火柴图像来说，S( r )的曲线没有很强的周期分量，而整齐摆放的火柴图像来说，S( r )的曲线在0~50内的某个值有很强的周期分量。在随机摆放火柴图中，S(θ)的曲线很明显，但是能量分布比较混乱，但是在整齐摆放的物体图像中，S(θ)只有在接近数值100时有很强的能量分布，重点性比较强烈。

（3）不变矩

对于平移、缩放、镜像和旋转都不敏感的7个二维不变矩的集合可以用以下公式推导出来：
可以使用函数invmoments来实现这7个不变矩，语法为：

phi = invmoments(f);  %f是输入图像，phi是包含7个元素的行向量 
1

编写代码：

f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0807(a).tif');
subplot(241),imshow(f),title('原始图像');
fp = padarray(f,[84 84],'both');    %填充从过的原始图像
subplot(242),imshow(fp),title('填充原始图像');
ftrans = zeros(568,568,'uint8');    
ftrans(151:550,151:550) = f;     %平移过的图像
subplot(243),imshow(ftrans),title('平移原始图像');
fhs = f(1:2:end,1:2:end);
fhsp = padarray(fhs,[184 184],'both');     %拥有一半尺寸且对应的填充图像
subplot(244),imshow(fhsp),title('填充一半尺寸的图像');
fm = fliplr(f);
fmp = padarray(fm,[84 84],'both');     %获取镜像图像
subplot(245),imshow(fmp),title('镜像图像');
fr45 = imrotate(f,45,'bilinear');    %旋转45°线性内插后的图像
fr90 = imrotate(f,90,'bilinear');   %旋转90°线性内插后的图像
subplot(246),imshow(fr45),title('旋转45°线性内插后的图像');
subplot(247),imshow(fr90),title('旋转90°线性内插后的图像');
fr90p = padarray(fr90,[84 84],'both');   %对旋转90°图像进行填充
subplot(248),imshow(fr90p),title('旋转90°线性内插后并填充的图像');
phi = invmoments(f);  %计算原始图像的不变矩
format short e   
phi
format short
phinorm = -sign(phi).*(log10(abs(phi)))     %使用log10变换减少它们的动态范围
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

（原始图像的不变矩）
（使用log10变换减少它们的动态范围后的不变矩）

（五）主分量描述

概念： 主分量变换，是指由原始图像数据协方差矩阵的特征值和特征向量建立起来的变换核，将光谱特征空间原始数据向量投影到平行于地物集群椭球体各结构轴的主成分方向，突出和保留主要地物类别信息，用来进行图像增强、特征选择和图像压缩的处理方法。

进行主分量变换等特征变换的目的：将原始图像通过一定的数字变换生成一组新的特征图像，这一组新图像信息集中在少数几个特征图像上，这样，数据量有所减少，达到消除相关系数，进行有效的特征选择和减少波段特征空间维数，达到数据压缩的目的。

对于任意给定的坐标对(i,j)，都有n个像素，每一幅图像在那个位置上都有一个像素。这些像素以列向量的形式排列，如图所示：

若这些图像的尺寸是M×N，那么在这n幅图像中，包含所有像素的N维向量总共有MN个。

主分量变换（也称为霍特林变换），公式为：

矩阵A的行是归一化为单位长度的Cx的特征向量。

主分量变换计算步骤：
1、 计算多光谱图像的均值向量M和协方差矩阵∑。
2、计算矩阵∑的特征值λr和特征向量φr（r=1，2，…，m），m为多光谱图像的波段数。
3、将特征值λ，按由大到小的次序排列，即λ1>λ2>…>λm。
4、选择前n个特征值对应的n个特征向量构造变换矩阵Φn。
5、根据Y=ΦnX进行变换，得到的新特征影像就是变换的结果，X为多光谱图像的一个光谱特征矢量。

一组n个已配准图像的堆栈形式：S = cat(3,f1,f2,...,fn)

这个大小为M×N×n的图像堆栈数组可用函数转换成行是n维向量的数组，语法为：[X,R] = imstack2vectors(S,MASK)

编写使用主分量的实验代码：

f1 = imread('C:\Users\Public\Pictures\Sample Pictures\Fig1130(a).tif');
f2 = imread('C:\Users\Public\Pictures\Sample Pictures\Fig1130(b).tif');
f3 = imread('C:\Users\Public\Pictures\Sample Pictures\Fig1130(c).tif');
f4 = imread('C:\Users\Public\Pictures\Sample Pictures\Fig1130(d).tif');
f5 = imread('C:\Users\Public\Pictures\Sample Pictures\Fig1130(e).tif');
f6 = imread('C:\Users\Public\Pictures\Sample Pictures\Fig1130(f).tif');
subplot(231),imshow(f1),title('第一幅图像');
subplot(232),imshow(f2),title('第二幅图像');
subplot(233),imshow(f3),title('第三幅图像');
subplot(234),imshow(f4),title('第四幅图像');
subplot(235),imshow(f5),title('第五幅图像');
subplot(236),imshow(f6),title('第六幅图像');
S = cat(3,f1,f2,f3,f4,f5,f6);
[X,R] = imstack2vectors(S);
P = principalcomps(X,6);
g1 = P.Y(:,1);
g1 = reshape(g1,512,512);
subplot(231),imshow(g1,[]),title('第一幅图像分量产生的结果');   %重复六次获得每幅图分量产生的结果
d = diag(P.Cy);
P = principalcomps(X,2);    %使用小一点的值然后进行主分量重构
h1 = P.X(:,1);
h1 = mat2gray(reshape(h1,512,512));
D1 = tofloat(f1) - h1;
subplot(231),imshow(D1,[]),title('比较第一幅图分量的差异');   %也是重复六次可以得到六幅图的差异
P.ems

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)))    %7个不变矩图像
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

实验分析： 观察六幅图的主分量变换图像，发现图像中的对比度细节的主要部分包含在前两幅图中，而之后的主分量变换图像的对比度逐渐减少，原因是由于前两个特征值与其他特征值的数值要大得多，因此主特征值对应的图像显示很高的对比度。对比同一幅图的差异图发现，前两幅图的差异比较小，而后面的图像差异比较大，这是由于方差大，因此差异就比较大。

编写用主分量调整物体的实验代码：

f1 = im2bw(imread('C:\Users\Public\Pictures\Sample Pictures\Fig1134(a).tif'));
f2 = im2bw(imread('C:\Users\Public\Pictures\Sample Pictures\Fig1134(b).tif'));
f3 = im2bw(imread('C:\Users\Public\Pictures\Sample Pictures\Fig1134(c).tif'));
subplot(131),imshow(f1),title('第一幅图像');
subplot(132),imshow(f2),title('第二幅图像');
subplot(133),imshow(f3),title('第三幅图像');
[x1 x2] = find(f1);
X = [x1 x2];
P = principalcomps(X,2);
A = P.A;
Y = (A*(X'))';
miny1 = min(Y(:,1));
miny2 = min(Y(:,2));
y1 = round(Y(:,1) - miny1+min(x1));
y2 = round(Y(:,2) - miny2+min(x2));
idx = sub2ind(size(f1),y1,y2);
fout = false(size(f1));
fout(idx) = 1;
fout = imclose(fout,ones(3));
fout = rot90(fout,2);
subplot(131),imshow(fout),title('第一幅图主分量的变换');
subplot(132),imshow(fout),title('第二幅图主分量的变换');  %将上面的f1换成f2即可
subplot(133),imshow(fout),title('第三幅图主分量的变换');    %将上面的f1换成f3即可
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

实验分析： 观察图像发现，这个实验主要是将歪曲的图像变成正面的图像，比如说字符A是歪的，变换后的图像是正的。一般来说，主分量变换沿着主宽度方向排列数据，但不能保证排列不会向相反的方向旋转，我们使用旋转数据将保证字符的正确旋转方向。这是使用主分量变换的方法将随机排列的字符图像变成垂直排列的字符，使用该方法可以简化后续的图像处理问题，因此要把它运用于数字图像处理当中。