[0030] 若ei的圆心不在ej之内且ej的圆心也不在ei之内，则ei与ej彼此互不包含，ei与ej均为中心椭圆；
[0031] 其中aj为第j个椭圆ej的长半轴；ai为第i个椭圆ei的长半轴；
[0032] 定义如下嵌套矩阵
[0033]
[0034] 其中，mij表示第i个椭圆ei与第j个椭圆ej的关系，即
[0035]
[0036] 根据上述嵌套矩阵滤除嵌套椭圆，得到中心椭圆，即如果mij值为1，则删除第i个椭圆ei；如果mij值为-1，则删除第j个椭圆ej，剩余全部为中心椭圆。
[0037] 如果上述各判断中任何一个不满足，则删除该椭圆。
[0038] 所述对中心椭圆进行参数修正包括以下步骤；
[0039] 按照编码点设计尺寸，根据当前椭圆D计算其编码带的内边界椭圆C和外边界椭圆A；并计算出位于内边界椭圆C和外边界椭圆A正中间的椭圆B；
[0040] 计算当前椭圆D以内的像素点灰度平均值，作为前景灰度阈值Tf；
[0041] 计算当前椭圆D以外和其内边界椭圆C以内的像素点灰度平均值，作为背景灰度阈值Tb；
[0042] 计算图像分割的阈值T＝(Tf+Tb)/2；
[0043] 把当前椭圆D以内像素点的灰度置为Tf；
[0044] 利用分割的阈值T对外边界椭圆A内的区域进行二值化，从而得到编码带的白色像素点，形成像素点集合S；
[0045] 利用椭圆的光栅化方法计算出正中间的椭圆B对应的轮廓EB，将像素点集合S和正中间的椭圆B对应的轮廓EB组成新的数据点集ESB；
[0046] 对数据集ESB进行椭圆拟合，得到新的椭圆B’；
[0047] 对前两个步骤进行迭代计算得到接近实际的编码带正中间的椭圆B’；
[0048] 由正中间的椭圆B’根据设计尺寸反算出当前椭圆D，作为修正后的中心椭圆D’。
[0049] 所述对编码点进行解码包括以下步骤：
[0050] 由修正后的中心椭圆D’计算出编码带的内边界椭圆C和外边界椭圆A，并计算出位于内边界椭圆C和外边界椭圆A正中间的椭圆B，得到正中间的椭圆B对应的轮廓EB；
[0051] 对轮廓EB中各点的灰度值进行线性滤波：
[0052] 设PB为轮廓EB上的一个点， 是过椭圆圆心O与PB的直线， 与内边界椭圆C和外边界椭圆A的交点分别为PC和PA，将线段 上的像素按灰度值进行排序，取灰度值居中的像素灰度作为PB的灰度值；
[0053] 对正中间的椭圆B上的点通过以下公式作逆变换：
[0054]
[0055] 式中，xo为正中间的椭圆B的圆心坐标，x为正中间的椭圆B上点PB的坐标，x′为与点PB对应的单位圆上的点的坐标，α为正中间的椭圆B的倾角；a、b分别为正中间的椭圆B的长半轴及短半轴；
[0056] 利用图像分割的阈值T，对单位圆上的点作二值化处理；
[0057] 利用二值化后的单位圆上的点集求解最优的初始搜索点Popt；
[0058] 从初始搜索点Popt开始，按顺时针方向在单位圆上每隔单位角度取一个二进制制位，计算制位内所有像素灰度的平均值，若所某位灰度平均值大于图像分割的阈值T，则取该位的二进制值为“1”，反之取“0”；
[0059] 依次处理完整个单位圆的圆周后，得到当前编码点的所有二进制编码；
[0060] 取上述所有二进制编码中最小的一个所对应的十进制数作为编码的ID号。
[0061] 所述求解最优的初始搜索点Popt包括以下步骤：
[0062] 对于单位圆上的任意一点Pi，计算出与它的前向k邻域{Pi-k，…，Pi-2，Pi-1}和后向k邻域{Pi+1，…，Pi-k-1，Pi+k}，其中，k的取值范围为
[0063]
[0064] 式中N为单位圆上总的点数，符号 表示向下取整；
[0065] 计算前向k邻域的平均灰度值Tpre和后向k邻域的平均灰度值Tpost；
[0066] 求Tpre与Tpost差的绝对值Tabs；
[0067] 取单位圆上Tpre与Tpost差的绝对值Tabs最大的点作为最优的初始搜索点。
[0068] 所述编码点由比色差明显的图案构成，每个图案中心为圆点，圆点周围为不连续的同心环形形成的编码带，每组编码点中，各编码点的编码均不相同；
[0069] 所述非编码点为周围无编码带的单纯圆点。
[0070] 本发明的有益效果是：
[0071] 1.鲁棒性强，精确度高。本发明采用鲁棒的Canny算子提取图像的边缘轮廓，采用尺寸准则、形状准则和嵌套性准则对各轮廓进行逐渐过滤，然后利用中心椭圆的编码带信息，对中心椭圆的参数进行重新估计和修正，提高了定位精度，在编码点解码时，采用径向线性中值滤波方法，消除了孤立噪声对编码带中间环形线上像素灰度的影响，并设计了一种新的最优初始搜索点估计策略，从而提高了整个编码点识别方法的鲁棒性，多次现场拍摄图像数据评测结果表明了本发明的优异性能。
[0072] 2.本发明可以直接嵌入到三维数字测量软件的数据预处理系统中，为后续的三维重建提供可靠的输入数据。
附图说明
[0073] 图1A～1D为本发明编码点和非编码点的基本规格形状示意图；
[0074] 图2为本发明编码点和非编码点在被测物体上的布置方法示意图；
[0075] 图3为本发明方法流程图；
[0076] 图4为本发明方法中提取中心椭圆的示意图；
[0077] 图5为本发明方法中中心椭圆修正示意图；
[0078] 图6为本发明方法中椭圆光栅化的处理区域示意图；
[0079] 图7为本发明方法中对轮廓EB中各点的灰度值进行线性滤波示意图；
[0080] 图8为本发明方法中最优初始搜索点选择示意图；
具体实施方式
[0081] 下面结合附图对本发明方法作进一步说明。
[0082] 如图3所示，本发明复杂背景下图像中鲁棒的编码点的精确识别与定位方法具有以下步骤：
[0083] 在被测物体表面均匀地布置多个不同的编码点和非编码点；
[0084] 利用摄像头或数码相机采集图像；
[0085] 利用上述采集的图像提取图像中编码点的中心椭圆；
[0086] 根据中心椭圆得到的该中心椭圆的编码带信息，并利用编码带信息对中心椭圆进行参数修正；
[0087] 根据修正后的中心椭圆参数对编码点进行解码，并输出各编码点的ID号。
[0088] 本发明在识别编码点时，只需用一台普通手持数码相机，一套编码点、非编码点在一台普通个人计算机上运行本发明方法程序即可完成。
[0089] 以下就图3中的几个关键步骤作详细说明：
[0090] 1.编码点的制作与布置
[0091] 本发明方法采用的编码点和非编码点如图1A～1D所示。每个编码点由于具有与其他编码点不同的图案，使之有唯一的ID号。本实施例中，如图1B所示，编码点外形尺寸为70×70mm，黑色无反光底色，编码点的中心为白色圆形图案，中心圆周围为与其同心的分段环状区域，用来确定编码点的身份信息，称为编码带；如图1A所示，本实施例将该圆环按角度平均分为12份，每份30°，相当于一个二进制位，每位可取白色(二进制码为“1”)或黑色(二进制码为“0”)；编码点各部分尺寸如图1C所示；如图1D所示，非编码点的尺寸为18×18mm，黑底色，中间的白色为直径10mm的圆形。编码点与非编码点用普通打印机打印出来后，均匀地布置在被测物体上，尽量避免编码点的局部被别的物体遮挡，否则容易导致误识别。
[0092] 2.拍摄图像
[0093] 本发明方法对数码相机的分辨率不作具体要求，目前，市场上数码相机的分辨率一般都在500万像素以上，完全可以满足拍摄要求。拍摄图像时角度尽量不要过偏或离物体太远，否则图像中编码点太扁或太小，影响识别效果。拍摄完成后，将照片输入到计算机中，进入到关键的软件识别阶段，主要包括提取中心椭圆、修正中心椭圆以及编码点解码等主要步骤。
[0094] 3.利用上述采集的图像提取图像中编码点的中心椭圆
[0095] 编码点和非编码点的中心圆经成像后变成了椭圆，在此，本发明方法的主要目的是在图像中提取出编码点的中心椭圆信息，如图4所示。
[0096] 首先利用鲁棒性较强的Canny算子提取出所有的边缘轮廓；
[0097] 以上述边缘轮廓的尺寸、形状以及嵌套性作为约束条件，对椭圆的边缘轮廓进行逐渐过滤，得到中心椭圆。
[0098] 本实施例采用下列准则逐渐对各轮廓进行过滤：
[0099] 3.1尺寸准则
[0100] 合格边缘轮廓的长度L应当满足如下条件
[0101] Ldown≤L≤Lup (1)
[0102] Ldown和Lup分别为事先设定的长度下限和长度上限阈值，阈值的设定一般取决于经验和识别任务。长度不满足式(1)的轮廓不是目标椭圆的轮廓，应当删除。
[0103] 3.2形状准则
[0104] 删除所有非封闭的轮廓。采用最小二乘方法分别对每条轮廓进行椭圆拟合。对于非常扁的椭圆，应当抛弃，因此，合格椭圆的长半轴a和短半轴b应满足下列关系[0105]
[0106] Raxis为长半轴和短半轴比值的阈值，一般取大于2的实数。
[0107] 合格椭圆的椭圆率ε还应该满足下列关系式
[0108] ε≤εellipticity (3)
[0109] εellipticity为椭圆率阈值。
[0110] 在本发明中，椭圆率是轮廓上各像素点到拟合椭圆的距离的均值。理想情况下，轮廓上的像素点完全落在拟合椭圆上，此时ε＝0。在实际应用中，由于离散误差的存在，ε不可能为0，因此可用ε的值来度量轮廓与椭圆的相似度。
[0111] 3.3嵌套性准则
[0112] 若某个椭圆的内部包含别的椭圆，则该椭圆肯定不是编码点的中心椭圆，称为嵌套椭圆，应当排除。经过尺寸准则和形状准则过滤后，对剩余的部分椭圆中任意两个椭圆ei和ej之间只存在三种位置关系：
[0113] 若aj
[0114] 若ai
[0115] 若ei的圆心不在ej之内且ej的圆心也不在ei之内，则ei与ej彼此互不包含，ei与ej均为中心椭圆；
[0116] 其中aj为第j个椭圆ej的长半轴；ai为第i个椭圆ei的长半轴；
[0117] 为此，定义如下嵌套矩阵
[0118]
[0119] 其中，mij表示第i个椭圆ei与第j个椭圆ej的关系，即
[0120]
[0121] 通过计算矩阵M，就可找出嵌套椭圆，即如果mij值为1，则删除第i个椭圆ei；如果mij值为-1，则删除第j个椭圆ej，剩余全部为中心椭圆。
[0122] M是一个反对称矩阵，因此只需计算上三角部分，下三角部分可通过反对称关系得到。
[0123] 4.对中心椭圆进行参数修正
[0124] 当图像中椭圆较小时，即用于拟合椭圆的轮廓像素数目较少，由于存在离散误差，拟合椭圆常常难以反映真实椭圆的形状和位置，若此时直接按设计尺寸关系用拟合椭圆估算编码带的位置，误差将被进一步放大，造成编码带与其实际位置严重偏离，最终导致识别错误。因此，必须对候选编码椭圆进行修正。对于当前中心椭圆D，如图5所示，其修正步骤如下：
[0125] 1)按照编码点设计尺寸，根据中心椭圆D计算其编码带的内边界椭圆C和外边界椭圆A；并计算出位于内边界椭圆C和外边界椭圆A正中间的椭圆B；
[0126] 2)计算当前中心椭圆D以内的像素点灰度平均值，作为前景灰度阈值Tf；
[0127] 3)计算当前中心椭圆D以外和其内边界椭圆C以内的像素点灰度平均值，作为背景灰度阈值Tb；
[0128] 4)计算图像分割的阈值T＝(Tf+Tb)/2；
[0129] 5)把当前中心椭圆D以内像素点的灰度置为Tf；
[0130] 6)利用分割的阈值T对外边界椭圆A内的区域进行二值化，从而得到编码带的白色像素点，形成像素点集合S；
[0131] 7)利用椭圆的光栅化方法计算出正中间的椭圆B对应的轮廓EB，将像素点集合S和正中间的椭圆B对应的轮廓EB组成新的数据点集ESB；
[0132] 8)对数据集ESB进行椭圆拟合，得到新的椭圆B’；
[0133] 9)对前两个步骤进行迭代计算得到接近实际的编码带正中间的椭圆B’；本实施例只需对步骤7)和8)5-7次迭代，就可得到满意的效果，此时的B’已经非常接近实际的编码带外边界椭圆；
[0134] 10)由正中间的椭圆B’根据设计尺寸反算出当前椭圆D，作为修正后的中心椭圆D’。
[0135] 在步骤7)中，由椭圆求对应的轮廓时需要用到椭圆的光栅化技术。为了便于描述，这里只考虑中心在原点，沿x方向的半轴为rx、沿y方向的半轴为ry的标准椭圆。以第一象限为例，光栅化算法分两部分应用于第一象限。图6中示例了依据椭圆斜率(rx
[0136] 椭圆光栅化方法可以概括为以下步骤：
[0137] a)取椭圆上的第一个点：
[0138] (x0，y0)＝(0，ry) (6)
[0139] b)计算区域1中决策参数的初始值p10：
[0140]
[0141] c)在区域1中的每个xk位置，从k＝0开始，完成下列测试：如果p1k<0，沿椭圆的下一个点为(xk+1，yk)，并且
[0142]
[0143] 否则，沿椭圆的下一个点为(xk+1，yk-1)，并且
[0144]
[0145] 其中：
[0146]
[0147] 并且直到
[0148] d)使用区域1中计算的最后点(x0，y0)来计算区域2中参数的初始值
[0149]
[0150] e)在区域2的每个yk位置处，从k＝0开始，完成下列测试：假如第k个参数值p2k>0，则沿椭圆的下一个点为(xk，yk-1)，并且
[0151]
[0152] 否则，沿椭圆的下一个点为(xk+1，yk-1)，并且
[0153]
[0154] 使用与区域1中相同的x和y增量进行计算。
[0155] f)确定其他三个象限中的对称点。
[0156] g)重复区域1的步骤，直到
[0157] 5.对编码点进行解码
[0158] 设计出编码点之后，每个编码点的图案都不一样，即每个编码点具有唯一的ID号。编码点解码是指根据编码椭圆识别出图像中每个编码点的ID号，为此，本发明采用如下策略：
[0159] 1)由修正后的中心椭圆D’计算出编码带的内边界椭圆C和外边界椭圆A，并计算出位于内边界椭圆C和外边界椭圆A正中间的椭圆B，得到正中间的椭圆B对应的轮廓EB；
[0160] 2)对轮廓EB中各点的灰度值进行线性滤波：
[0161] 如图7所示，设PB为轮廓EB上的一个点， 是过椭圆圆心O与PB的直线， 与内边界椭圆C和外边界椭圆A的交点分别为PC和PA，将线段 上的像素按灰度值进行排序，取灰度值居中的像素灰度作为PB的灰度值。
[0162] 这种方法相当于对点PB在一个线性窗口内进行一次中值滤波。经中值滤波后，可以消除孤立噪声的影响，提高识别算法的鲁棒性。
[0163] 3)对正中间的椭圆B上的点通过以下公式作逆变换：
[0164]
[0165] 式中，xo为正中间的椭圆B的圆心坐标，x为正中间的椭圆B上点PB的坐标，x′为与点PB对应的单位圆上的点的坐标，α为正中间的椭圆B的倾角；a、b分别为正中间的椭圆B的长半轴及短半轴；
[0166] 4)利用图像分割的阈值T，对单位圆上的点作二值化处理；
[0167] 5)利用二值化后的单位圆上的点集求解最优的初始搜索点Popt；
[0168] 6)从初始搜索点Popt开始，按顺时针方向在单位圆上每隔单位角度(30°)取一个二进制制位，计算制位内所有像素灰度的平均值，若某位灰度平均值大于图像分割的阈值T，则取该位的二进制值为“1”，反之取“0”；
[0169] 依次处理完整个单位圆的圆周后，得到当前编码点的所有二进制编码；选取的初始搜索点不同，得出的二进制编码也不一样。由于每个编码只对应一个十进制数，因此它不可能对应12个不同的二进制码。为此，取12个二进制码中最小的一个所对应的十进制数作为编码的ID。本实施例事先将所有的编码点的ID都保存在一个查表中，若某个椭圆的ID与查表的某个ID对应，说明该椭圆是编码点，否则就是非编码点。
[0170] 7)输出识别出的编码点ID。
[0171] 对于步骤5)，如图8所示，设X为像素点序列，Y为像素点灰度，T为灰度阈值；
[0172] 一般的做法是任选一个点作为初始搜索点，这是很危险的。比如，当起始点恰好选在单位圆上某个二进制位的中央附近时(如图8所示的三角形所处的位置)，该位本为“1”或“0”，由于噪声的存在，很可能被错判为“0”或“1”，从而导致识别错误。因此，初始搜索点的选择不能是随意的。直观地讲，如果把单位圆上每个点的灰度看作是一个离散信号，那么初始搜索点的最优位置应该是信号的突变位置。但是，当噪声比较大的时候，突变点并不明显，可能存在一定的漂移。为此，本发明方法提出了如下的初始搜索点优化策略：
[0173] a)对于单位圆上的任意一点Pi，计算出与它的前向k邻域{Pi-k，…，Pi-2，Pi-1}和后向k邻域{Pi+1，…，Pi-k-1，Pi+k}，其中，k的取值范围为
[0174]
[0175] 式中，N是单位圆上总的点数，符号 表示向下取整。在取值范围内，k值越大，起始点越逼近最优值，但计算时间越长，一般地，在满足识别要求的前提下，尽量取较小的k值；
[0176] b)计算前向k邻域的平均灰度值Tpre和后向k邻域的平均灰度值Tpost；
[0177] c)求Tpre与Tpost差的绝对值Tabs；
[0178] d)取单位圆上Tpre与Tpost差的绝对值Tabs最大的点作为最优的初始搜索点。
[0179] 综上所述，本发明方法能在具有复杂背景的图像中鲁棒地识别出编码点，在实验中获得了优异的识别和定位精度，具有广泛的应用前景。
