基于劣弧演化的点到隐式曲线距离计算方法

摘要

本发明提供一种平面内的点到隐式曲线距离的数值计算方法，具有精确、鲁棒、算法效率高以及对隐式函数的光滑性要求不高等优点。具体步骤是：从一个以所述点为圆心、与隐式曲线相交的圆开始，先不断缩小半径直至同心圆与隐式曲线不相交，再不断二分同心圆的半径，让不相交的内圆与相交的外圆无限接近，最终求得距离。其中圆与隐式曲线的相交性判断采用劣弧演化算法，即对圆周的稠密等分点按其二进制编码逆序码的顺序，判断函数值是否与圆心异号，若异号，则判为相交，不然，根据函数值对劣弧集内的圆弧进行删除、收缩或分裂等演化处理，缩小劣弧集的覆盖范围，直至劣弧集变空，判为不相交。

说明书

技术领域

本发明涉及数值计算领域，具体地说是一种为电子计算机提供点到隐式曲线距离的数值计算方法。

背景技术

点到隐式曲线的距离计算方法在模式识别、几何建模、计算机视觉、医学影像处理等领域具有重要应用价值，也是数值分析、计算机科学、逆向工程等学科的重要研究基础。由于点到隐式曲线的距离的解析式一般很难得到甚至不可能得到，所以距离的计算通常是数值计算。

设所求距离是点p₀到隐式曲线f(x，y)＝0的距离，点到隐式曲线距离的数值计算通常有局部法和全局法之分。局部法含有一个迭代过程，如牛顿迭代法，它先以迭代方式求出隐式曲线上与p₀达到最小距离的点(称为足点)，也就是先求解下列关于p的方程的根p^*：

$>\left(\begin{array}{c}(p-p_0)\times ▿f\left(p\right)=0\\ f\left(p\right)=0\end{array}\right)$>

其中，表示梯度算子，×为矢量叉乘。于是所求距离就是||p*-p₀||。迭代法收敛速度较快，但经常受奇点和初值选择的影响，存在不收敛、陷入局部极值点等缺点，且对f(x，y)的光滑性要求较高。全局法则根据隐式函数j(x，y)的某种全局特性对足点进行二维搜索，如文献“Geometric constraint solver using multivariate rational spline functions”(Elber，G.and Kim，M.-S.，In：Proceedings of the Sixth ACM Symposium on Solid Modeling and Applications，ACM，2001：pp.1-10)和文献“Computation of the solutions of nonlinear polynomial systems”(Sherbrooke，E.C.and Patrikalakis，N.M.，Computer Aided Geometric Design，v10，1993：pp.379-405)所述方法，根据基函数的凸包性质，采用Bézier裁剪方法不断排除非足点区域，得到候选足点，最后在这些候选足点中选取与p₀距离最小的点。全局法能保证鲁棒计算，但计算量较大，对隐式函数有特殊要求(如凸性)。

此外，还有作为全局方法和局部方法之折中的方法，基本框架仍为迭代法，但每步迭代根据径向导数搜索下一个足点的估计位置，参见文献“Robust Computation of Foot Points onImplicitly Defined Curves”(Aigner，M.and Jittler，B.，In：Mathematical Methods for Curves andSurfaces：2004，Nashboro Press，2005：pp.1-10)。

当然，由于收敛速度和鲁棒计算很难两全，在许多需要计算点到隐式曲线距离的应用场合常常采用估算法，如基于曲率的二阶估计、单纯形估计。具体可参见文献“Fitting B-splinecurves to point clouds by curvature-based squared distance minimization”(Wang，Wenping；Pottmann，H.；et a1.，ACM Transactions on Graphics，25(2)，2006：214-238)和文献“Implicitpolynomial representation through a fast fitting error estimation”(Rouhani，M.；Sappa，A.D.，IEEE Transactions on Image Processing，21(4)，2012：2089-2098)。距离的估计不是本发明的设计目标，本发明的目标是精确计算与鲁棒性。

发明内容

为了精确计算平面内的点到隐式曲线的距离，达到事先给定的任意小的容许误差，并且不受曲线几何形状的影响，避免落入局部极值、依赖初值选择等问题，本发明采用如下技术方案。

给定二维平面内的一个点p₀及一条隐式曲线，其中隐式曲线为某有界区域Ω上有定义且方向导数有界的二元函数f(x，y)依其零值集确定的单分支平面曲线，本发明的点p₀到所述隐式曲线距离的数值计算方法，包含如下步骤：

寻找一个函数值f(p)与f(p₀)异号或等于零的初始点p，此时以p₀为圆心、以||p₀-p||为半径的圆必与隐式曲线相交；置内半径为0，外半径为||p₀-p||；再取中值半径为内半径与外半径之平均值，然后判断以p₀为圆心、以中值半径为半径的圆与隐式曲线是否相交，若相交，外半径下调为中值半径，若不相交，内半径上调为中值半径；再取中值半径、判断相交性、调整内半径或外半径，如此反复，直至外半径与内半径之差小于容许误差，返回内半径与外半径之平均值作为点p₀到所述隐式曲线的距离。

在上述数值计算过程中，会产生一系列的同心圆。一开始外圆与隐式曲线相交，内圆为一个退化的点，中值半径为外半径的1/2，以此半径作新的同心圆，若该同心圆与隐式曲线相交，则再取外半径的1/4，若还相交，再取1/8……，所以一开始有一个圆的倍缩过程。

上述数值计算的整个过程始终保持外圆与隐式曲线相交、内圆与隐式曲线不相交的状态。整个循环过程就是通过在外圆和内圆之间不断插入二分的同心圆，在一直保持外圆与隐式曲线相交、内圆与隐式曲线不相交的前提下让内圆与外圆无限接近。

更形式化的步骤描述如下：

(a1)取f的函数值f(p)与f(p₀)异号或零值的点p；

(a2)取初始内半径为r₁＝0，初始外半径为r₂＝||p₀-p||；

(a3)r＝(r₁+r₂)/2；

(a4)判断以点p₀为圆心、r为半径的圆与隐式曲线是否相交；

(a5)若相交，置r₂＝r，否则，置r₁＝r；

(a6)若r₂-r₁＜ε(ε为容许误差)，返回(r₁+r₂)/2作为点p₀到所述隐式曲线的距离，算法结束，否则，转至步骤(a3)。

在本发明的点到隐式曲线距离的数值计算方法的实施方案中，判断圆与隐式曲线是否相交的方法包含以下步骤：

(b1)取最小的正整数n，使圆周2ⁿ等分后每段弧长小于容许误差ε，圆周等分后依次分布的等分点为$>q_0,q_1,......,q_{2^n-1},q_{2^n}=q_0.$>

(b2)置初始值i＝0，初始劣弧集A＝{[0，2ⁿ]}，劣弧集是圆上函数值可能与f(p₀)异号的圆弧的集合，其元素形如[s，t]，其中s、t分别是圆弧两端的等分点标号，[s，t]表示由等分点q_s，q_s+1，……，q_t连成的圆弧，初始劣弧集只有一个从q₀连到的单个圆弧，其实为整个圆周。

(b3)将i的n位二进制数表示b₁b₂……b_n逆序，得到二进制数b_n……b₂b₁，其表示的整数为j。

(b4)若q_j不落在劣弧集A的任一圆弧内，跳至步骤(b7)。

(b5)计算f(q_j)，若该函数值与f(p₀)异号，结束判断，返回结果为“圆与隐式曲线相交”，否则，进行劣弧演化，即根据f(q_j)值与方向导数上界M排除f(q)肯定与f(p₀)同号的等分点q，对A中的圆弧进行删除、收缩和分裂处理，缩小劣弧集的覆盖范围。

(b6)若劣弧集结束判断算法，返回结果为“圆与隐式曲线不相交”。

(b7)置i为i+1，若i＝2ⁿ，结束判断算法，返回结果为“圆与隐式曲线不相交”，否则转至步骤(b3)。

其中步骤(b5)中的劣弧演化可以更加具体化，具体包含以下步骤：

(c1)若|f(q_j)|＞2Mr，置结束；

其中：

r为圆的半径；

M为f(x，y)的方向导数的上界，即

$>M=\underset{\Omega }{\text{max}}\left|\right|▿f(x,y)\left|\right|=\underset{\Omega }{\max }\left|\right|(\frac{\partial f(x,y)}{\partial x},\frac{\partial f(x,y)}{\partial y})\left|\right|<+\infty .$>

(c2)计算$>a=2\arcsin \frac{\left|f\left(q_j\right)\right|}{2\mathrm{Mr}}.$>

(c3)计算其中表示取下整数，此时f(q_i)必与f(p₀)同号(j-k≤i≤j+k)，即得2k+1个同号等分点。

(c4)若j＝0，直接将圆弧[0，2ⁿ]收缩为[k+1，2ⁿ-k-1]，否则，对劣弧集A的每个圆弧[s，t]，分如下五种情形分别作演化处理：

(c4.1)j-k＞t或j+k＜s：保持[s，t]不变；

(c4.2)j+k≥t≥j-k＞s：[s，t]收缩为[s，j-k-1]；

(c4.3)t＞j+k≥s≥j-k：[s，t]收缩为[j+k+1，t]；

(c4.4)j+k≥t≥s≥j-k：将[s，t]从A中删除；

(c4.5)t＞j+k≥j-k＞s：[s，t]分裂为[s，j-k-1]和[j+k+1，t]两个圆弧。

(c5)结束。

因为初始劣弧集的唯一元素为一个圆周，首次演化(i＝0)后成为真正的圆弧，但有可能是优弧(圆心角大于π)，并不一定是劣弧。但是在第二次演化(i＝1)时，由于本发明特殊的j值计算方法，由i＝1＝(00...001)₂，得j＝(100...00)₂＝2ⁿ/2，第2个处理的等分点q_j必为第1个处理等分点q₀的对径点。若该点为同号点，则算法结束，否则进行第2次演化，演化后的两个圆弧(若有的话)必为劣弧。这说明第2次演化之后劣弧集内都一直只有真正的劣弧了，而且这些劣弧不断地收缩、分裂或消失，总覆盖范围不断缩小，故称本算法为劣弧演化算法。

本发明的技术效果：

①精确性和鲁棒性。本发明的计算方案可达到事先给定的任意小的容许误差，并且不受曲线几何形状的影响，不存在落入局部极值、依赖初值选择、不收敛等问题。

②算法效率高。表面上看本发明的算法是逐圈搜索、逐点扫描的“蛮力法”，但实际算法速度并不慢。圆与隐式曲线的相交判断是算法效率的主要决定因素，无论圆与隐式曲线相交还是不相交，均能快速判断。具体地讲，在圆与隐式曲线相交的情况下，由于采用二进制逆序码(如步骤(b3))技术，算法会按等分点逐渐均匀加密的顺序处理，可以证明能在O(1)时间内遇到异号点而从步骤(b5)结束；在圆与隐式曲线不相交的情况下，算法会随着劣弧演化的进行，劣弧集快速变空而从步骤(b6)结束。所以算法会在很短时间内结束相交性判断，从而快速求得点到隐式曲线的距离。

③对隐式曲线的函数f(x，y)的光滑性要求不高。牛顿迭代法要求f(x，y)的一阶微分存在且连续，其它全局法要求f(x，y)具备某种全局性质，而本算法只要求可导，且方向导数有界，上界为M，实际上本算法适用的条件可以更松，只需满足Lipschitz连续条件：|f(p)-f(q)|≤M||p-q||。

④由于本发明的算法总是从异号点开始，先执行一个“不断缩小，直至遇到不相交同心圆”的倍缩过程，所以不可能存在“不断扩大，直至遇到相交同心圆”的方法所可能存在的圆半径扩大一倍就越过整个曲线分支的“扑空”现象。

附图说明

图1是本发明较佳实施例的算法总流程图。

图2是本发明较佳实施例中的圆与隐式曲线相交判断算法流程图。

图3是本发明较佳实施例中的劣弧演化流程图。

图4是本发明较佳实施例中的劣弧演化示意图。

图5是本发明较佳实施例中不相交情形的劣弧演化过程实例图。

图6是本发明较佳实施例中相交情形的劣弧演化过程实例图。

具体实施方式

给定二维平面内的一个点p₀＝(x₀，y₀)及一条有界区域Ω上有定义的隐式曲线f(x，y)＝0，点p₀到隐式曲线距离的数值计算的总体步骤如下：

(1)取f的函数值f(p)与f(p₀)异号或零值的点p；

(2)取初始内半径为r₁＝0，初始外半径为r₂＝||p₀-p||；

(3)r＝(r₁+r₂)/2；

(4)判断以点p₀为圆心、r为半径的圆与隐式曲线是否相交；

(5)若相交，置r₂＝r，否则，置r₁＝r；

(6)若r₂-r₁＜ε(ε为容许误差)，返回(r₁+r₂)/2作为点p₀到所述隐式曲线的距离，算法结束，否则，转至步骤(3)。

参见图1，其为本发明较佳实施例的算法总流程图。为求点p₀与隐式曲线f(x，y)＝0的距离，首先，执行步骤101，寻找一个f(p)与f(p₀)异号或等于零的点p，即若f(p₀)＞0，取f(p)≤0，若f(p₀)＜0，取f(p)≥0(f(p₀)＝0的情况不必考虑，此时点在隐式曲线上，所求距离为0)。寻找异号点方法有随机投点、直接观察取特定点等方法，后面另有说明。接着执行步骤102，给内半径r₁和外半径r₂分别赋初值为0和||p₀-p||，以p₀为圆心、r₁为半径的圆(此时退化为点)与隐式曲线不相交，而以r₂为半径的同心圆与隐式曲线相交。然后进入一个循环过程，先执行步骤103，取中值半径r，再通过步骤104判断以r为半径的圆与隐式曲线是否相交，若相交，执行步骤105，重置外半径r₂为r，内半径r₁保持不变，否则，执行步骤106，重置内半径r₁为r，外半径r₂保持不变。内半径或外半径重置之后，执行步骤107，判断r₂-r₁是否小于容许误差ε，若r₂-r₁＜ε，执行步骤108，返回r₁与r₂的平均值作为所求的距离，算法结束，否则，转至步骤103继续循环执行。

在上述实施步骤101中，寻找异号点很容易实现。至少有两个办法，第一个办法是随机选点，不妨设在区域Ω内f(p)的正值区域和负值区域的面积比为u：v(0＜u，v＜1，u+v＝1)，则n次随机选点得到f(p)非负点的概率为1-vⁿ，得到非正点的概率为1-uⁿ，两者均会随n增大而快速趋于1。此外，如果点到隐式曲线距离的计算是针对一组点{p₀，p₁，p₂，……}进行计算的话，在为第一个点p₀计算距离的过程中得到的任意一个非负值点和非正值点可以给组中的其它点p₁、p₂、……作为初始异号点使用。若碰到一个f(q)＝0的点，则这样的点既可作为非负值点又可作为非正值点。第二个办法是直接观察，对大部分具有特定形式的隐式函数f(x，y)，比如二元多项式、张量积B-样条函数，可根据函数f(p)＝f(x，y)在特殊方向(如坐标轴)上的单调性、渐进性经观察或取特定值就直接寻到非负值点、非正值点，甚至零值点。

参见图2，其为本发明较佳实施例中的圆与隐式曲线相交判断算法流程图，是图1步骤104的进一步细化。首先执行步骤201，取最小的正整数n，使圆周2ⁿ等分后的每段弧长小于容许误差ε，即n满足

$>\frac{2\mathrm{\pi r}}{2^n}<ϵ\le \frac{2\mathrm{\pi r}}{2^{n-1}}$>

其中r为所述圆的半径。

取圆周上的等分点为

$>q_i=(x_0+r\cos \frac{2\mathrm{\pi i}}{2^n},y_0+r\sin \frac{2\mathrm{\pi i}}{2^n}),(i=\mathrm{0,1},...,2^n-1)$>

其中r和p₀＝(x₀，y₀)分别为所述圆的半径和圆心，这些等分点用一个数组q[0..2ⁿ-1]表示。然后执行步骤202，置i的初值为0，置劣弧集A的初始状态为{[0，2ⁿ]}，即只有一个圆弧，它其实是整个圆周的“圆弧”。执行步骤203，计算i作为n位二进制数的逆序数，若i＝(b₁b₂……b_n)₂，则j＝(b_n……b₂b₁)₂。接着执行步骤204，判断等分点q[j]是否落在A的某个圆弧内，若不落在任一圆弧，执行步骤211与步骤212：i增1并判断i是否达到最大值2ⁿ，否则执行下面一系列步骤。首先是步骤205，计算函数值f(q[j])，然后是步骤206，判断该函数值是否与f(p₀)异号。若异号(包括f(q[j])＝0)，执行步骤207，返回判断结果“圆与隐式曲线相交”，算法结束；不然，f(q[j])与f(p₀)严格同号，执行步骤208，进行劣弧演化。劣弧演化是根据f(q[j])值与方向导数上界M排除f(q)肯定与f(p₀)同号的等分点q，对A中的圆弧进行删除、收缩和分裂处理，缩小劣弧集的覆盖范围。劣弧演化之后是步骤209，判断劣弧集是否为空，若为空，返回判断结果“圆与隐式曲线不相交”，算法结束，否则，执行步骤211与步骤212：i增1并判断i是否到达最大值2ⁿ。

若在步骤211与步骤212遇到i增1后达到最大值2ⁿ，则说明所有等分点处理完毕，并没有遇到异号点，于是执行步骤210，返回判断结果“圆与隐式曲线不相交”，算法结束。

通过步骤212→210结束算法的情况是不理想的，因为算法处理了所有等分点，效率比较低。实际上，本判断算法在很大的概率上不会从步骤212→210结束，而通过步骤206→207或者步骤209→210结束，在圆与隐式曲线相交的情况下会在O(1)时间内通过步骤206→207结束，这一点由本发明的等分点特殊处理顺序(见步骤203)所保证；在圆与隐式曲线不相交的情况，则会随着劣弧演化(步骤208)的进行，劣弧集的覆盖范围快速缩小，最后变空，于是通过步骤209→210结束。所以判断效率的期望值很高。

如图3所示，其为劣弧演化的较佳实施步骤，是图2步骤208的进一步细化。此时，f(q[j])与f(p₀)严格同号，首先执行步骤301，判断|f(q[j]|＞2Mr，其中：

r为圆的半径；

M为函数f(x，y)的方向导数的上界，即：

$>M=\underset{\Omega }{\text{max}}\left|\right|▿f(x,y)\left|\right|=\underset{\Omega }{\max }\left|\right|(\frac{\partial f(x,y)}{\partial x},\frac{\partial f(x,y)}{\partial y})\left|\right|<+\infty .$>

方向导数的上界M对各种典型的曲线拟合用的二元函数来说很容易估计。事实上，z＝f(x，y)可以看成是一个显式表示的空间曲面，所以，对二元多项式f(x，y)，M可根据多项式系数估计；对Bézier曲面、张量积B-样条、NURBS等参数曲面，M可根据基函数和控制参数估计(参见文献《计算机辅助几何设计》第二章2.3-2.7，王国瑾、汪国昭、郑建民著，高等教育出版社，2001)；对根据距离场估计而拟合得到的曲面，则更容易估计，因为准确的距离场具有方向导数之模处处为1的特点(参见文献《基于有向距离场的代数B-样条曲线重建》，李云夕、冯结青、金小刚，《软件学报》第18卷第9期第2306-2317页，2007)，所以M可以估计为一个比1稍大的数，比如1.5。

即使遇到很难估计的情况，M估计为

$>M=\frac{v_{\max }-v_{\min }}{D}\times 20$>

其中：

v_max为f(x，y)在区域Ω中的最大值；

v_min为f(x，y)在区域Ω中的最小值；

D为区域Ω的尺寸(圆形区域的直径，矩形区域的长宽平均值)。

如果f(x，y)的最大值与最小值还很难估计，则M取大数1000。

关于M的估计的两点说明：第一，M的估计只影响本发明算法的圆的等分点的采样密度，纵然估计不足，仅在隐式曲线的尖锐处的计算有可能存在偏差。可以想象，如果在所有等分点q_j，函数值f(q_j)与f(p₀)同号，据本发明误判为圆与隐式曲线不相交，而实际上相交，这说明隐式函数f(x，y)有一个狭长的异号部分通过两个相邻的等分点之间“伸入”圆内部，所以，曲线存在一个尖锐特征，而这样的隐式曲线对曲线拟合和曲线造型来说是不稳定的，在实际几何建模的技术研究中本来就需要避免。第二，假如做最保守的估计M(比如M＝1000)，只影响到本发明圆与隐式曲线不相交时的判断效率，相交情形的判断效率不受影响。

再回到图3，若步骤301判断结果为真，说明整个圆处在f(x，y)的同号区域中，则执行步骤302，置劣弧集A为空集，结束；否则，执行步骤303和步骤304，确定必定与f(p₀)同号的等分点的数目，步骤303是计算同号弧的弧度(单侧)，步骤304是计算同号等分点的个数k(单侧)，于是以q[j]为中心的两侧共有2k+1个同号等分点(包括q[j]本身)，标号从j-k到j+k。接着根据同号弧[j-k，j+k]对劣弧集进行演化处理，首先是执行步骤305，看是否在处理第一个等分点，在处理第一个等分点时，注意到同号弧覆盖q[0]两侧，所以执行步骤306，圆周经收缩后变为圆弧[k+1，2ⁿ-k-1]；处理以后的等分点时，执行步骤307，对劣弧集A中的每个圆弧[s，t]，分5种情形分别处理，如步骤308～312。5种情形如图4所示，列出了不同形状的同号弧对同一条劣弧的各种演化示意图，其中粗线弧401代表当前劣弧[s，t]，细线弧402代表同号弧[j-k，j+k]。

情形1(如图4(a))：j-k＞t或j+k＜s，说明同号弧与当前劣弧不交，保持[s，t]不变。

情形2(如图4(b))：j+k≥t≥j-k＞s，说明同号弧与当前劣弧尾部重叠，将[s，t]收缩为[s，j-k-1]。

情形3(如图4(c))：t＞j+k≥s≥j-k，说明同号弧与当前劣弧首部重叠，将[s，t]收缩为[j+k+1，t]。

情形4(如图4(d))：j+k≥t≥s≥j-k，说明同号弧覆盖了整个当前劣弧，将[s，t]从A中删除。

情形5(如图4(e))：t＞j+k≥j-k＞s，说明同号弧落在当前劣弧内部，将[s，t]分裂为[s，j-k-1]和[j+k+1，t]两个劣弧。

对A中每个劣弧处理完毕，结束。

图5为不相交情形的劣弧演化过程实例图，示出了圆与隐式曲线不相交情况下劣弧演化的分步动态变化过程，其中n＝4，等分点数是2n＝16，假定j＝0，1，2，……，15时同号弧[j-k，j+k]中的k值(根据f函数值确定)分别是3，3，4，3，2，1，0，1，1，0，0，1，1，2，3，3。

(a)i＝0，j＝0，计算f(q[j])得k＝3，收缩[0，16]，得A＝{[4，12]}；

(b)i＝1，j＝8，计算f(q[j])得k＝1，分裂[5，11]，得A＝{[4，6]，[10，12]}；

(c)i＝2，j＝4，计算f(q[j])得k＝2，删除[4，6]，得A＝{[10，12]}；

(d)i＝3，j＝12，计算f(q[j])得k＝1，收缩[10，12]，得A＝{[10，10]}；

(e)i＝4，j＝2，j不落在任一劣弧内；

(f)i＝5，j＝10，计算f(q[j])得k＝0，删除[10，10]，得返回结果“不相交”，结束。

按顺序逐点判断的方法需要16次f函数值计算，而采用劣弧演化只需要5次计算。

图6为相交情形的劣弧演化过程实例图，示出了圆与隐式曲线相交情况下劣弧演化的分步动态变化过程，其中n＝4，等分点数是2ⁿ＝16，假定j＝0，1，2，……，15时同号弧[j-k，j+k]中的k值(根据f函数值确定)分别是2，2，1，0，--，--，--，--，0，--，--，--，0，1，2，2。(--表示异号)

(a)i＝0，j＝0，计算f(q[j])得k＝2，收缩[0，16]，得A＝{[3，13]}；

(b)i＝1，j＝8，计算f(q[j])得k＝0，分裂[3，13]，得A＝{[3，7]，[9，13]}；

(c)i＝2，j＝4，计算f(q[j])得到异号点，返回结果“相交”，结束。

按顺序逐点判断的方法平均需要8次f函数值计算，而采用劣弧演化只需要3次计算。