§10.1 图象特征参数
图象的参数测量基于图象的边界提取和跟踪。在对二值化图象进行边界提取后,目标区域的边界点坐标信息以结点的形式保存于链表中,然后以某个结点为边界的起始点,对链表中的结点进行遍历,寻找与该结点所对应象素点属于同一区域边界的结点,并重新排列结点顺序,直到完成对所有链表结点的遍历。通过边界提取和跟踪,将图象中各个有效区域的边界点信息以链表的形式表达,不仅可以实现参数测量的数据上的准备,而且在数据结构上也是一种有效的表达方式,有利于采用快速算法进行数据处理。
图象的参数反映了图象中的各种特征。例如,在图象中有一些水果,如果希望从该图象中提取出香蕉,对于计算机而言,它无法分辩出哪一个是香蕉,我们只能通过提供所要提取物体的一些特征来指示计算机,如香蕉的大小、形状等,使它依据这些特征来寻找满足条件的物体。
图象的特征参数包括:区域面积、周长、圆形度、重心、孔径等。
⑴ 图象的区域面积
区域面积的定义为:目标区域中所包含的象素数。
⑵ 图象的周长
图象的周长定义为:目标区域的边界线上象素间距离之和。其中,象素间的距离有以下两种情形:
① 边界上的两个象素点A与B以如下两种方式排列:
<!--[if !vml]--><!--[endif]-->
若象素点A的坐标为 <!--[if !vml]--><!--[endif]-->,象素点B的坐标为 <!--[if !vml]--><!--[endif]-->,那么象素A,
B间的距离为:
<!--[if !vml]--><!--[endif]-->
通常,可以认为以这两种方式排列的象素间距离(d)为1个象素。
② 边界上的两个象素点A与B以如下两种方式排列:
<!--[if !vml]--><!--[endif]-->
同①,象素A,B间的距离为:
<!--[if !vml]--><!--[endif]--> (10-1)
根据上图中两个象素的排列方式,也可以认为象素间距离(d)为 <!--[if !vml]--><!--[endif]-->个
象素。
在进行图象的周长测量时,顺序遍历属于同一区域的边界点坐标,根据上述两种情形下的象素点排列方式,分别计算象素间的距离,并将结果进行累加以得到目标区域的周长。
⑶ 图象的圆形度
图象的圆形度定义为:在目标区域面积和周长的基础上,计算该区域的形状复杂程度的特征量。圆形度的数学表达式为:
<!--[if !vml]--><!--[endif]--> (10-2)
其中, <!--[if !vml]--><!--[endif]-->为目标区域的圆形度, <!--[if !vml]--><!--[endif]-->为区域面积, <!--[if !vml]--><!--[endif]-->为区域周长。
若目标区域为圆形,其半径为 <!--[if !vml]--><!--[endif]-->,那么 <!--[if !vml]--><!--[endif]-->, <!--[if !vml]--><!--[endif]-->,圆形度为 <!--[if !vml]--><!--[endif]-->;若目标区域为正方形,其边长为 <!--[if !vml]--><!--[endif]-->,那么 <!--[if !vml]--><!--[endif]-->, <!--[if !vml]--><!--[endif]-->,圆形度为 <!--[if !vml]--><!--[endif]-->;若目标区域为等边三角形,其边长为 <!--[if !vml]--><!--[endif]-->,那么 <!--[if !vml]--><!--[endif]-->, <!--[if !vml]--><!--[endif]-->,圆形度为 <!--[if !vml]--><!--[endif]-->。
因此,圆形度反映了区域的形状复杂程度:区域形状越接近于圆, <!--[if !vml]--><!--[endif]-->越大,其最大值为1;区域形状越复杂, <!--[if !vml]--><!--[endif]-->越小,其最小值趋向0。
⑷ 图象的重心
图象的重心定义为目标区域中所有象素坐标的平均值。若某区域由 <!--[if !vml]--><!--[endif]-->个象素点组成,象素点坐标为 <!--[if !vml]--><!--[endif]-->,那么其重心坐标 <!--[if !vml]--><!--[endif]-->由下式求得:
<!--[if !vml]--><!--[endif]--> (10-3)
在本文中,测量的目标区域为近似圆形,除上述四个参数以外,还需进行孔径的测量。若目标区域为理想的圆形( <!--[if !vml]--><!--[endif]-->),那么孔径的计算方法为:若遍历链表结点,对属于同一圆形区域的边界点进行周长计算,求得周长为 <!--[if !vml]--><!--[endif]-->,那么目标区域的孔径( <!--[if !vml]--><!--[endif]-->)的计算公式为:
<!--[if !vml]--><!--[endif]--> (10-4)
在一般情况下,目标区域多为近似圆形,由于目标区域本身精度或噪声干扰影响,其圆形度 <!--[if !vml]--><!--[endif]-->,无法直接使用上述孔径计算公式。因此,我们需先对目标区域的边界进行圆拟合,根据拟合结果进行目标区域的孔径计算。圆拟合的数学公式如下所示:
<!--[if !vml]--><!--[endif]--> (10-5)
其中,目标区域的边界点坐标为 <!--[if !vml]--><!--[endif]-->,拟合后的圆心坐标为 <!--[if !vml]--><!--[endif]-->,拟合半径为 <!--[if !vml]--><!--[endif]-->。
§10.2 区域面积测量算法实现
图象的区域面积参数测量由函数getArea实现。函数getArea中入口参数定义如下:
⑴ head linked型形参,此为输入多区域边界点链表
⑵ nBegin 整型形参,此为输入链表结点位置头索引
⑶ nEnd 整型形参,此为输入链表结点位置尾索引
函数返回值为:
nResult 整型类型,此为返回由索引nBegin和nEnd指定边界所在区域的象素点统计数(区域面积)
函数getArea的流程图如图10.1所示:
<!--[if !vml]--><!--[endif]-->
图10.1 getArea流程图
其中,链表结点排序模块通过调用函数getSelectSort,实现对给定区域边界点的排序。目标区域的面积参数计算是基于如下方法:对某一区域,在跟踪提取边界后,统计边界所包含的象素点数。由于目标区域以链表结点的形式表达,因此,区域内部象素点数目的统计通过结点(边界点)坐标计算而实现。具体实现方法为:
对属于同一区域的边界点A和B,其坐标分别为 <!--[if !vml]--><!--[endif]-->和 <!--[if !vml]--><!--[endif]-->,若有如下关系式:
<!--[if !vml]--><!--[endif]-->,且 <!--[if !vml]--><!--[endif]--> (10-6)
那么,由这两个边界点所确定区域的象素点数为:
<!--[if !vml]--><!--[endif]-->
所以,需对链表进行遍历,以找到满足上述关系式的边界点。若链表为非排序表,其算法时间复杂度为 <!--[if !vml]--><!--[endif]-->,搜索时间会因图象数据量增加而呈平方级增长。因此,依据上述关系式对链表结点进行排序,可使算法的时间复杂度降为 <!--[if !vml]--><!--[endif]-->,提高程序运行效率。
函数getSelectSort实现链表结点的排序,采用的排序算法为选择排序。该算法适用于数据项比较大,键又比较小的随机文件。选择排序算法的执行过程为:首先,选出数组中最小的项,将它与数组的第一个成员交换位置。然后选出次小的项,将它与数组的第二个成员交换位置。按这种算法一直做下去,直到整个数组排序完毕。
对数据项较大键较小的文件,选择排序的运行时间是一个线性数。若 <!--[if !vml]--><!--[endif]-->表示数据项大小相对于键大小的比例,那么假设比较操作的花费为1单位时间,则交换操作的花费为 <!--[if !vml]--><!--[endif]-->单位时间。选择排序在比较操作上的花费约为 <!--[if !vml]--><!--[endif]-->单位时间,而交换操作的花费则是 <!--[if !vml]--><!--[endif]-->单位时间。如果 <!--[if !vml]--><!--[endif]-->比 <!--[if !vml]--><!--[endif]-->的常数倍要大,则 <!--[if !vml]--><!--[endif]-->超过 <!--[if !vml]--><!--[endif]-->,这种情况下运行时间与 <!--[if !vml]--><!--[endif]-->成比例,即和移动所有数据所需要的时间总量成比例。
函数getSelectSort中入口参数定义如下:
⑴ head linked型形参,此为输入待排序的链表
⑵ nBegin 整型形参,此为输入待排序的链表的起始结点位置
⑶ nEnd 整型形参,此为输入待排序的链表的终止结点位置
函数返回值:void
函数getSelectSort的部分实现代码如下所示:
for(i = nBegin, k = nBegin+1; i < nEnd; i++, k++)
for(j = k; j <= nEnd; j++)
if(head[i].item_x < head[j].item_x){ //①
exch(&(head[i].item_x), &(head[j].item_x)); //②
exch(&(head[i].item_y), &(head[j].item_y));
}else{
if(head[i].item_x == head[j].item_x){
if(head[i].item_y < head[j].item_y){
exch(&(head[i].item_x), &(head[j].item_x));
exch(&(head[i].item_y), &(head[j].item_y));
}
}
}
该代码段实现选择排序算法。其中注释①所指的if-else条件语句实现边界点坐标判断,即规定了链表结点排序准则:优先考虑结点横坐标大小,若某个结点的横坐标小于其后面结点的横坐标,那么交换这两个结点中的横纵坐标值;若两个结点的横坐标相等,那么比较它们的纵坐标大小,若该结点的纵坐标小于其后面结点的纵坐标,那么交换这两个结点的横纵坐标值。遍历整个链表结点完成选择排序。注释②处的函数exch实现不同结点间横纵坐标值的交换。
区域面积计算模块实现对目标区域中象素数的统计。目标区域由若干象素点组成,在目标区域范围内,每一对横坐标相等,纵坐标不等的边界点之间包含了一定数目的象素点,通过统计若干组具有这种特征的边界点间的象素数,并将结果累加即可求得区域面积。区域面积计算示例如图10.2所示:
<!--[if !vml]--><!--[endif]-->
图10.2区域面积计算示例
在上图所示区域中,边界点坐标为 <!--[if !vml]--><!--[endif]-->,区域面积计算方法如下:
① 边界点 <!--[if !vml]--><!--[endif]-->与 <!--[if !vml]--><!--[endif]-->之间,存在 <!--[if !vml]--><!--[endif]-->,象素数为: <!--[if !vml]--><!--[endif]-->
② 边界点 <!--[if !vml]--><!--[endif]-->与 <!--[if !vml]--><!--[endif]-->之间,存在 <!--[if !vml]--><!--[endif]-->,象素数为: <!--[if !vml]--><!--[endif]-->
③ 边界点 <!--[if !vml]--><!--[endif]-->与 <!--[if !vml]--><!--[endif]-->之间,存在 <!--[if !vml]--><!--[endif]-->,象素数为: <!--[if !vml]--><!--[endif]-->
④ 边界点 <!--[if !vml]--><!--[endif]-->与 <!--[if !vml]--><!--[endif]-->之间,存在 <!--[if !vml]--><!--[endif]-->,象素数为: <!--[if !vml]--><!--[endif]-->
目标区域面积 <!--[if !vml]--><!--[endif]-->为:
<!--[if !vml]--><!--[endif]--> (10-7)
区域面积计算模块的部分实现代码如下所示:
for(i = nBegin, nIndex = nBegin; i <= nEnd; i++)
if((i == nEnd) || ((i < nEnd) && (head[i].item_x != head[i+1].item_x))){//①
nResult += head[nIndex].item_y-head[i].item_y+1;
nIndex=i+1;
}
注释①处的if判断语句实现目标区域的象素点统计。其中,nBegin和nEnd为链表中某一区域边界的起始点和终止点位置,nResult为目标区域面积。
<!--[if !supportEmptyParas]--> <!--[endif]-->
§10.3 区域周长测量算法实现
图象的周长参数测量由函数getPerim实现。函数getPerim中入口参数定义如下:
⑴ head linked型形参,此为输入多区域边界点链表
⑵ nBegin 整型形参,此为输入链表结点位置头索引
⑶ nEnd 整型形参,此为输入链表结点位置尾索引
函数返回值:
fResult float类型,此为返回由索引nBegin和nEnd指定区域的周长
函数getPerim的流程图如图10.3所示:
<!--[if !vml]--><!--[endif]-->
图10.3 getPerim流程图
其中,链表结点排序模块通过调用函数getLinkedSort,实现对给定区域边界点的排序。其排序方法为:针对某一结点head+i(区域边界点),以该结点为起始点遍历整个链表,若找到这样一个结点head+j,它位于起始点head+i的8-邻域范围内,那么将该结点与起始点之后的结点head+i+1进行互换,再以该结点为起始点进行链表遍历,寻找它的8-邻域点,依次进行下去,直到到达链表末尾。如图10.4所示:
<!--[if !vml]--><!--[endif]-->
图10.4链表结点排序
函数getLinkedSort中入口参数定义如下:
⑴ head linked型形参,此为输入待排序的链表
⑵ nBegin 整型形参,此为输入待排序的链表的起始结点位置
⑶ nEnd 整型形参,此为输入待排序的链表的终止结点位置
函数返回值:void
函数getLinkedSort的部分实现代码如下所示:
for(i = nBegin, k = nBegin+1; i < nEnd; i++, k++)
for(j = k; j <= nEnd; j++)
if(isLinked(head[i], head[j])){ //①
exch(&(head[i+1].item_x), &(head[j].item_x)); //②
exch(&(head[i+1].item_y), &(head[j].item_y));
break;
}
其中,注释①处的if判断语句通过调用函数isLinked进行边界点的8-邻域点判断,若边界点head[j]为head[i]的8-邻域边界点,那么通过调用注释②处所指函数exch,交换链表结点head[i+1]与head[j]的横纵坐标值,否则继续遍历链表或以下一个结点为起始点进行8-邻域点判断。
函数isLinked实现边界点的8-邻域点判断。它的入口参数定义如下:
⑴ x node型形参,此为输入8-邻域的中心点
⑵ t node型形参,此为输入待判断的邻域点
函数返回值:整数类型
函数isLinked通过一系列if判断语句,对边界点t是否位于中心点x的8-邻域内进行判断,并返回相应的判断值。如图10.5所示:
<!--[if !vml]--><!--[endif]-->
图10.5边界点的8-邻域点判断
将边界点t的横纵坐标值以顺时针方向与中心点 <!--[if !vml]--><!--[endif]-->的8-邻域点坐标进行比较,若坐标值相等,那么边界点t即为中心点的邻域点,并再以边界点t为中心点进行8-邻域点比较。
区域周长计算模块在完成链表结点排序后,实现目标区域的周长计算。链表结点排序完成对目标区域边界点的跟踪,使得在给定区域边界起始点后,可以连续追踪下一个区域边界点,直到遍历所有边界点。因此,区域周长计算由整个目标区域缩小为对两个相邻边界点间距离的计算,并对所有两点间距离进行累加,最终得到目标区域周长。如图10.6所示:
<!--[if !supportEmptyParas]--> <!--[endif]-->
<!--[if !supportEmptyParas]--> <!--[endif]-->
<!--[if !supportEmptyParas]--> <!--[endif]-->
<!--[if !supportEmptyParas]--> <!--[endif]-->
<!--[if !vml]--><!--[endif]-->
图10.6区域周长计算
其中,目标区域边界点为 <!--[if !vml]--><!--[endif]-->,共9个。若以边界点 <!--[if !vml]--><!--[endif]-->为起始点 <!--[if !vml]--><!--[endif]-->进行边界跟踪,那么跟踪方向为逆时针方向,所得链表结点排列顺序为:
<!--[if !vml]--><!--[endif]-->
因此,可以通过区域周长计算模块求得其周长C为:
<!--[if !vml]--><!--[endif]--> (10-8)
其中, <!--[if !vml]--><!--[endif]-->, <!--[if !vml]--><!--[endif]-->
区域周长计算模块的部分实现代码如下所示:
for(i = nBegin; i <= nEnd; i++)
if((head[i].item_x != head[i+1].item_x) //①
&& (head[i].item_y != head[i+1].item_y))
fResult += 1.414;
else
fResult += 1;
其中,注释①处的if-else判断语句实现区域边界周长计算。nBegin与nEnd为链表中某一区域边界点的起始位置和终止位置,在完成从nBegin到nEnd的边界点遍历后,fResult即为最终的目标区域边界周长值。
<!--[if !supportEmptyParas]--> <!--[endif]-->
§10.4 圆形度测量算法实现
图象的圆形度参数测量由函数getCompact实现。函数getCompact中入口参数定义如下:
⑴ nArea 整型形参,此为输入目标区域面积
⑵ fPerim double型形参,此为输入目标区域周长
函数返回值:
fResult float类型,此为返回目标区域的圆形度
函数getCompact的实现代码如下所示:
double getCompact(int nArea, double fPerim)
{
double fResult=0.0;
return fResult=4*3.1415926*nArea/(fPerim*fPerim); //①
}
注释①处的代码实现圆形度计算,其所依据的圆形度计算公式如下所示:
<!--[if !vml]--><!--[endif]--> (10-9)
其中, <!--[if !vml]--><!--[endif]-->为目标区域的圆形度, <!--[if !vml]--><!--[endif]-->为区域面积, <!--[if !vml]--><!--[endif]-->为区域周长。
通过目标区域的圆形度计算,我们可以了解目标区域的边界复杂程度,并通过比较不同区域的圆形度,选取合适的目标区域。如图10.7是三种不同图象的圆形度对比:
<!--[if !supportEmptyParas]--> <!--[endif]-->
<!--[if !vml]--><!--[endif]-->
图10.7不同图象圆形度对比
由上表可得,圆形度可作为目标区域提取的依据。若所提取区域为圆形,则选取圆形度最接近于1的目标区域,并结合其他参数(如面积,周长等)最终确定目标区域。
<!--[if !supportEmptyParas]--> <!--[endif]-->
§10.5 重心测量算法实现
图象的重心参数测量由函数getCentroid实现。函数getCentroid中入口参数定义如下:
⑴ head linked型形参,此为输入目标区域边界点链表
⑵ nBegin 整型形参,此为输入目标区域边界起始点位置
⑶ nEnd 整型形参,此为输入目标区域边界终止点位置
⑷ nArea 整型形参,此为输入目标区域的面积
函数返回值:
fResult Centr类型,此为返回目标区域的重心坐标
其中,Centr类型为自定义类型,其定义如下所示:
struct centroid{
double x;
double y;
};
typedef struct centroid Centr;
结构centroid中两个元素x和y分别代表重心的横坐标和纵坐标,并将该结构通过typedef命名为Centr类型。
在进行重心参数计算时,均匀物体的重心和形心重合。因此,通过形心计算公式可以直接求得均匀物体的重心坐标。对于形心公式可通过如下例子进行解释。
一个质量忽略不计的刚性杆下有一个支点,在支点两边分别在杆上放有6个质量均匀,尺寸相同,边长为L的正方体。如图10.8所示:
<!--[if !vml]--><!--[endif]-->
图10.8重心位置计算
设正方体的重量为W,当杆平衡时,力矩平衡方程为:
<!--[if !vml]--><!--[endif]-->
所以,可求得 <!--[if !vml]--><!--[endif]-->, <!--[if !vml]--><!--[endif]-->即为右边物体在水平方向上的重心位置。
对于平面图象的形心计算公式如下所示:
<!--[if !vml]--><!--[endif]-->, <!--[if !vml]--><!--[endif]--> (10-10)
其中,形心坐标为 <!--[if !vml]--><!--[endif]-->, <!--[if !vml]--><!--[endif]-->表示图象内部的象素点坐标, <!--[if !vml]--><!--[endif]-->为图象面积。若将一个象素作为一个 <!--[if !vml]--><!--[endif]-->,且为1个单位,则图象的面积即为象素点的个数,分子上的积分则为图象内部象素坐标的和。因此,在给定目标区域的面积和区域内部象素点坐标后,目标区域的重心坐标可根据上式求得。
函数getCentroid的部分实现代码如下所示:
for(i = nBegin, nIndex = nBegin; i <= nEnd; i++){
if((i == nEnd)||((i < nEnd)&&(head[i].item_x != head[i+1].item_x))){ //①
k=head[nIndex].item_y-head[i].item_y+1;
fResult.x+=head[i].item_x*k;
fResult.y+=head[i].item_y*k+((k-1)*k)/2;
nIndex=i+1;
}
}
fResult.x /= nArea; //②
fResult.y /= nArea; //③
其中,注释①处的if判断语句实现对目标区域内部象素点的横坐标和纵坐标累加,最终累加结果存放于fResult.x和fResult.y中。注释②和③所指语句依据重心计算公式求得目标区域重心坐标(fResult.x,fResult.y)。
<!--[if !supportEmptyParas]--> <!--[endif]-->
§10.6 孔径测量算法实现
图象的孔径参数测量由函数getDiam实现。函数getDiam中入口参数定义如下:
⑴ head linked型形参,此为输入多区域边界点链表
⑵ nBegin 整型形参,此为输入链表中目标区域边界点起始位置
⑶ nEnd 整型形参,此为输入链表中目标区域边界点终止位置
函数返回值:
dResult double类型,此为返回目标区域的孔径大小
函数getDiam的流程图如图10.9所示:
<!--[if !vml]--><!--[endif]-->
图10.9 getDiam流程图
其中,提取区域边界点坐标模块根据入口参数nBegin和nEnd,确定边界点链表head中目标区域边界点的起始和终止位置,并由此提取相应的坐标值,为后继边界圆拟合做准备。
区域边界圆拟合及孔径计算模块中,通过调用函数getApprox实现边界圆拟合。由于存在噪声干扰以及目标区域本身加工误差等原因,无法直接使用孔径计算公式:
<!--[if !vml]--><!--[endif]--> (10-11)
其中, <!--[if !vml]--><!--[endif]-->为目标区域孔径, <!--[if !vml]--><!--[endif]-->为目标区域周长。因此,需要对目标区域边界进行曲线拟合,使其依据关系式:
<!--[if !vml]--><!--[endif]--> (10-12)
拟合为圆,进而求得目标区域的孔径。
曲线拟合根据一组二维数据,即平面上的若干点,确定一个一元函数 <!--[if !vml]--><!--[endif]-->,即拟合曲线,使这些数据点与曲线总体上尽量接近。最常用的曲线拟合方法是最小二乘法。其原理如下:
求一元函数 <!--[if !vml]--><!--[endif]-->,使 <!--[if !vml]--><!--[endif]-->为最小。如图10.10所示:
<!--[if !vml]--><!--[endif]-->
图10.10曲线拟合示例
其中, <!--[if !vml]--><!--[endif]-->为点 <!--[if !vml]--><!--[endif]-->与曲线 <!--[if !vml]--><!--[endif]-->的距离。曲线拟合即寻求一个函数 <!--[if !vml]--><!--[endif]-->,使 <!--[if !vml]--><!--[endif]-->在某种准则下与所有数据点最为接近。最小二乘法就是使所有散点到曲线的距离平方和最小。拟合时选用一定的拟合函数 <!--[if !vml]--><!--[endif]-->形式,设拟合函数可由一些简单的“基函数”(如幂函数,三角函数等) <!--[if !vml]--><!--[endif]-->来线性表示,如下所示:
<!--[if !vml]--><!--[endif]--> (10-13)
将上式代入 <!--[if !vml]--><!--[endif]-->中,可得 <!--[if !vml]--><!--[endif]-->,即 <!--[if !vml]--><!--[endif]-->成为
<!--[if !vml]--><!--[endif]-->的函数。若要使 <!--[if !vml]--><!--[endif]-->为最小,则令 <!--[if !vml]--><!--[endif]-->对 <!--[if !vml]--><!--[endif]--> <!--[if !vml]--><!--[endif]-->的偏导数等于零,于是得到 <!--[if !vml]--><!--[endif]-->个方程组,从中求解出 <!--[if !vml]--><!--[endif]--> <!--[if !vml]--><!--[endif]-->。若取基函数为 <!--[if !vml]--><!--[endif]-->,那么拟合函数 <!--[if !vml]--><!--[endif]-->为多项式函数。
图10.11为微孔边界象素点集合,以及进行边界圆拟合后圆心坐标位置,如下所示:
<!--[if !vml]--><!--[endif]-->
图10.11边界圆拟合
函数getApprox实现目标区域边界圆拟合。它的入口参数定义如下:
⑴ iX int*型形参,此为输入目标区域边界点横坐标数组首地址
⑵ iY int*型形参,此为输入目标区域边界点纵坐标数组首地址
⑶ iSize_Coord 整型形参,此为输入数组iX、iY的大小
⑷ iCenter int*型形参,此为输入存放整型元素的数组首地址,该数组用于存放目标区域边界拟合圆的圆心坐标
⑸ iSize_Center 整型形参,此为输入数组iCenter的大小
⑹ dR double*型形参,此为输入double型变量的地址,该变量用于存放目标区域边界拟合圆的半径
函数返回值:void
函数getApprox的流程图如图10.12所示:
<!--[if !vml]--><!--[endif]-->
图10.12 getApprox流程图
其中,函数getMean实现边界点坐标的平均值计算。函数getMean中入口参数定义如下:
⑴ piItem int*型形参,此为输入待求平均值的数组元素首地址
⑵ iSize 整型形参,此为输入数组大小
函数返回值:
double类型,此为返回数组元素的平均值
函数getMean的部分实现代码如下所示:
for(i = 0; i < iSize; i++){
dTemp += piItem[i]; //①
}
return dTemp/iSize; //②
其中,注释①处代码实现对输入的数组元素进行累加,并将累加和存入dTemp,注释②处代码实现数组元素平均值计算。计算公式如下所示:
<!--[if !vml]--><!--[endif]--> (10-14)
其中, <!--[if !vml]--><!--[endif]-->为数组元素平均值, <!--[if !vml]--><!--[endif]--> <!--[if !vml]--><!--[endif]-->为数组元素, <!--[if !vml]--><!--[endif]-->为数组元素个数。
在边界点坐标平均值计算中,dMean_x, dMean_y, dMean_z分别存放x,y,z方向边界点坐标平均值。由于目标区域为二维图象,所以z方向坐标平均值为0。
边界点坐标与均值之差模块的实现代码如下所示:
for(i = 0; i < iSize_Coord; i++){
pdX2[i]=iX[i]-dMean_x;
pdY2[i]=iY[i]-dMean_y;
pdZ2[i]=iZ[i]-dMean_z;
}
其中,iX[i],iY[i],iZ[i] (i=0,1, …, iSize_Coord-1)为边界点x,y,z方向上的坐标,且iZ[i]=0 (i=0,1, …, iSize_Coord-1)。
中间变量dY6,dZ6_1,dZ6_2计算模块的实现代码如下所示:
for(i = 0; i < iSize_Coord; i++){
dY6+=iY[i]*pdY2[i];
dZ6_1+=iZ[i]*pdY2[i];
dZ6_2+=iZ[i]*pdZ2[i];
}
其中,pdY2[i],pdZ2[i] (i=0,1, …, iSize_Coord-1)分别为边界点坐标与均值之差模块中的计算结果。
边界点坐标分量平方和计算模块的实现代码如下所示:
for(i = 0; i < iSize_Coord; i++){
pdA[i]=iX[i]*iX[i]+iY[i]*iY[i]+iZ[i]*iZ[i];
}
其中,iX[i],iY[i],iZ[i] (i=0,1, …, iSize_Coord-1)为边界点x,y,z方向上的坐标,且iZ[i]=0 (i=0,1, …, iSize_Coord-1)。
中间变量dX9,dY9,dZ9计算模块的实现代码如下所示:
for(i = 0; i < iSize_Coord; i++){
dX9+=pdX2[i]*pdA[i];
dY9+=pdY2[i]*pdA[i];
dZ9+=pdZ2[i]*pdA[i];
}
其中,pdX2[i],pdY2[i],pdZ2[i] (i=0,1, …, iSize_Coord-1)为边界点坐标与均值之差模块计算结果,pdA[i] (i=0,1, …, iSize_Coord-1)为边界点坐标分量平方和计算模块的结果。
建立3×1矩阵pdB模块,中间变量dX4,dY4,dZ4计算模块和建立3×3矩阵ppdArray10模块将上述模块的计算结果以矩阵形式进行表述,为矩阵计算及拟合圆半径与圆心坐标计算做准备。
矩阵ppdArray10求逆计算模块通过调用函数getInvMatrix实现矩阵求逆运算。函数getInvMatrix中入口参数定义如下:
⑴ ppdMatrix double*型形参,此为输入原始矩阵 <!--[if !vml]--><!--[endif]-->的元素;返回时存放其
逆矩阵 <!--[if !vml]--><!--[endif]-->的元素
⑵ iLevel 整型形参,此为输入矩阵 <!--[if !vml]--><!--[endif]-->的阶数
函数返回值:
整数类型,该返回值为一个整型标志值,若返回的标志值为0,则表示矩阵 <!--[if !vml]--><!--[endif]-->奇异,并输出相应信息;若返回的标志值不为0,则表示正常返回。
函数getInvMatrix采用高斯-约当(Gauss-Jordan)消去法求 <!--[if !vml]--><!--[endif]-->阶实矩阵 <!--[if !vml]--><!--[endif]-->的逆矩阵 <!--[if !vml]--><!--[endif]-->。高斯-约当法求逆矩阵的步骤如下。
首先,对于 <!--[if !vml]--><!--[endif]-->从 <!--[if !vml]--><!--[endif]-->到 <!--[if !vml]--><!--[endif]-->作如下几步:
⑴ 从第 <!--[if !vml]--><!--[endif]-->行、第 <!--[if !vml]--><!--[endif]-->列开始的右下角子阵中选取绝对值最大的元素,并记住此元素所在的行号和列号,再通过行交换与列交换将它交换到主元素位置上。这一步为全选主元。
⑵ <!--[if !vml]--><!--[endif]-->
⑶ <!--[if !vml]--><!--[endif]-->, <!--[if !vml]--><!--[endif]-->; <!--[if !vml]--><!--[endif]-->
⑷ <!--[if !vml]--><!--[endif]-->, <!--[if !vml]--><!--[endif]-->; <!--[if !vml]--><!--[endif]-->
⑸ <!--[if !vml]--><!--[endif]-->, <!--[if !vml]--><!--[endif]-->; <!--[if !vml]--><!--[endif]-->
最后,根据在全选主元过程中所记录的行、列交换的信息进行恢复,恢复的原则为:在全选主元过程中,先交换的行、列后进行恢复;原来的行(列)交换用列(行)交换来恢复。
矩阵pdB与ppdArray10的逆矩阵相乘模块通过调用函数getMulpMatrix实现矩阵相乘运算。函数getMulpMatrix中入口参数定义如下:
⑴ ppdMatr_A double*型形参,此为输入矩阵 <!--[if !vml]--><!--[endif]-->的元素
⑵ ppdMatr_B double*型形参,此为输入矩阵 <!--[if !vml]--><!--[endif]-->的元素
⑶ iRow_A 整型形参,此为输入矩阵 <!--[if !vml]--><!--[endif]-->的行数
⑷ iCol_A 整型形参,此为输入矩阵 <!--[if !vml]--><!--[endif]-->的列数
⑸ iCol_B 整型形参,此为输入矩阵 <!--[if !vml]--><!--[endif]-->的列数
⑹ ppdMatr_C double*型形参,返回乘积矩阵 <!--[if !vml]--><!--[endif]-->的元素
函数返回值:void
函数getMulpMatrix对 <!--[if !vml]--><!--[endif]-->阶矩阵 <!--[if !vml]--><!--[endif]-->与 <!--[if !vml]--><!--[endif]-->阶矩阵 <!--[if !vml]--><!--[endif]-->进行相乘运算,所得乘积矩阵为 <!--[if !vml]--><!--[endif]-->。矩阵相乘运算公式如下:
<!--[if !vml]--><!--[endif]-->, <!--[if !vml]--><!--[endif]-->; <!--[if !vml]--><!--[endif]--> (10-15)
其中, <!--[if !vml]--><!--[endif]-->为矩阵 <!--[if !vml]--><!--[endif]-->中的元素, <!--[if !vml]--><!--[endif]-->为矩阵 <!--[if !vml]--><!--[endif]-->中的元素, <!--[if !vml]--><!--[endif]-->为矩阵 <!--[if !vml]--><!--[endif]-->中的元素。
拟合圆半径与圆心坐标计算模块根据上述模块结果,实现半径和圆心坐标计算。该模块的部分实现代码如下所示:
for(i = 0; i < iSize_Coord; i++){ //①
dTemp1=iX[i]-ppdArray12[0][0];
dTemp2=iY[i]-ppdArray12[1][0];
dTemp3=iZ[i]-ppdArray12[2][0];
dX16+=dTemp1*dTemp1;
dY16+=dTemp2*dTemp2;
dZ16+=dTemp3*dTemp3;
}
*dR=sqrt((dX16+dY16+dZ16)/iSize_Coord); //②
iCenter[0]=(int)ppdArray12[0][0];
iCenter[1]=(int)ppdArray12[1][0];
iCenter[2]=(int)ppdArray12[2][0];
其中,注释①处的for循环语句内,iX[i],iY[i],iZ[i] (i=0,1, …, iSize_Coord-1) 为边界点x,y,z方向上的坐标,且iZ[i]=0 (i=0,1, …, iSize_Coord-1);矩阵ppdArray12保存了矩阵pdB与ppdArray10的逆矩阵相乘的运算结果,即为拟合圆的圆心坐标 <!--[if !vml]--><!--[endif]--> <!--[if !vml]--><!--[endif]-->。
注释②处的代码实现拟合圆半径计算。其中,dX16为x方向边界点坐标与圆心坐标之差的平方的累加和;dY16为y方向边界点坐标与圆心坐标之差的平方的累加和;dZ16为z方向边界点坐标与圆心坐标之差的平方的累加和,即为0。通过上述运算,求得拟合圆半径 <!--[if !vml]--><!--[endif]-->在x和y两个方向上的分量 <!--[if !vml]--><!--[endif]-->和 <!--[if !vml]--><!--[endif]-->,如图10.13所示:
<!--[if !supportEmptyParas]--> <!--[endif]-->
<!--[if !vml]--><!--[endif]-->
图10.13拟合圆半径计算示例
所以,拟合圆半径为: <!--[if !vml]--><!--[endif]-->。
<!--[if !supportEmptyParas]--> <!--[endif]-->
§10.7 图象标定
在进行区域测量时所得到的面积、周长、孔径等参数值基于实际目标区域在CCD光敏单元上的面积,即象素数。因此,为了得到所测参数的实际值,必须在图象测量之前,确定图象(象素)与真实尺寸之间的比例,即进行图象标定。
在进行图象标定之前,通过CCD摄像头在测量显微镜上捕获标准图象,如图10.14所示:
<!--[if !vml]--><!--[endif]-->
图10.14图象标定
其中,D为标准图象中十字两个端线间的实际距离。在捕获到静态标准图象后,软件中通过鼠标操作可以获取距离为D的端线间象素数P,因此,在测量显微镜的某一放大倍率下,标定常数计算公式为:
<!--[if !vml]--><!--[endif]--> (10-16)
根据上述计算公式,可分别求得水平方向和竖直方向上的标定常数 <!--[if !vml]--><!--[endif]-->与 <!--[if !vml]--><!--[endif]-->,进而可求得面积标定常数为:
<!--[if !vml]--><!--[endif]--> (10-17)
由此,我们可以计算得到目标区域的参数实际值。