基于DP曲线构建PH样条

Construction of PH Splines Based on DP Curves

摘要: PH曲线是一类特殊的多项式参数曲线，其最显著的优点是弧长函数为多项式，其等距线可由兼容于CAD系统的有理多项式曲线表示。鉴于此，基于三次DP曲线，从平面PH曲线的定义出发，给出了三次DP曲线为PH曲线时其控制多边形满足的充分必要条件，得到了关于控制多边形的边长和夹角的几何特征条件，给出DP-PH曲线的定义。通过DP-PH控制多边形几何特征条件，给出可以先通过构造控制多边形进而构造出三次DP-PH曲线的几何构造方法。进一步分析了DP曲线和DP-PH曲线的误差。

Abstract: The PH curve is a special kind of polynomial parameter curve. The most significant advantage is that the arc length function is a polynomial. And the equidistant line can be represented by a rational polynomial curve, which is compatible with the CAD system. Based on the DP curve, a necessary and sufficient condition for a cubic plane DP curve to be a PH curve is obtained. The geometrical characteristics of the edge length and the angle of the control polygon are obtained and the definition of the DP-PH curve is obtained. The geometric feature condition of polygon is controlled by DP-PH. We describe the construction of a control polygon for a cubic DP-PH curve from geometric construction method, based on the procedure for a DP curve. Then the errors between DP curve and DP-PH curve are analyzed.

1. 引言

在计算机辅助几何设计中，非均匀有理B样条(NURBS)方法发展比较成熟，成为现代曲线曲面设计中最为广泛流行的技术，但是仍然存在不足：不能准确表示摆线、螺旋线、圆锥、圆弧曲线等工程问题中经常用到的超越曲线曲面。近年来，学者们不仅在非多项式空间构造出了新型曲线，而且在多项式空间也构造出了新型曲线。Delgado和Peña [1] 提出了新型参数曲线，称为DP曲线，该曲线不但在数值计算上具有稳定性、在算法上具有线性的计算时间复杂度，而且是由具有曲线保形性的全正基(NTP基)生成的 [2] [3] [4] [5]。

在CAD的很多领域涉及到曲线弧长和等距线的计算问题，例如铁路与公路的设计、机械零件的设计以及半智能机器人运动轨迹生成等，往往需要多项式曲线的弧长和等距线具备有理形式。为此，FAROUKI和SAKKALIS [6] 引入毕达哥拉斯速端(Pythagorean Hodograph, PH)曲线，具有多项式形式的弧长和有理形式的等距线。到目前为止，对PH曲线的研究过多地侧重于代数结构方面，在几何方面的研究成果很匮乏。众所周知，对于给定Bézier曲线的控制多边形，其相关边长和内角不依赖于坐标选择的固有内在几何参量，具体数据可由实际测量得。因此，从控制多边形的长度和角度来讨论PH曲线的几何性质无论在理论上还是实际应用中都具有重要意义。Farouki和Sakalis [6] 给出具有不同控制顶点的3次Bézier曲线为PH曲线分离形式的边长约束条件和角的约束条件。此外，五次PH曲线的几何特征条件也被确立，但对于边长和角度来说这个条件不是分离形式的。近年来，四次 [7]、五次 [8] 和六次 [9] PH曲线的充分必要的几何特征条件被确定。

对于基于三次DP曲线构造PH曲线的问题，本文给出了三次PH-DP曲线的定义，DP曲线成为PH曲线的充分必要条件，关于控制多边形的边长和夹角的几何特征条件，并且是边角完全分离的条件。通过DP-PH控制多边形几何特征条件，我们给出可以先通过构造控制多边形进而构造出三次DP-PH曲线的几何构造方法。最后，分析了DP曲线和DP-PH曲线的误差，给出数值例子。

2. 三次DP-PH曲线的构造

本章基于DP曲线，构造PH曲线(DP-PH曲线)，并给出三次DP曲线成为PH曲线的几何特征条件。

2.1. 三次DP曲线

给定控制点 $P_k$，对于 $\forall t\in \left[0,1\right]$，定义三次DP曲线：

$P\left(t\right)={\displaystyle \underset{i}{\overset{3}{\sum }}{p}_i}{D}_{i,3}(t)$

其中三次DP曲线基函数为：

$\left[D_{0,3}\left(t\right),D_{1,3}\left(t\right),D_{2,3}\left(t\right),D_{3,3}\left(t\right)\right]=\left[{\left(1-t\right)}^3,t\left(1-t\right)\left(2-t\right),t\left(1-t\right)\left(1+t\right),t^3\right]$

2.2. PH曲线

定义1 [6] ：给定平面参数多项式曲线 $q\left(t\right)=\left(x\left(t\right),y\left(t\right)\right)$，如果存在多项式 $\delta \left(t\right)$，使得，即 $\delta \left(t\right)=\sqrt{{x^{\prime }}^2\left(t\right)+{y^{\prime }}^2\left(t\right)}$ 是勾股数，称平面参数多项式曲线 $q\left(t\right)$ 为Pythagorean-hodograph曲线，简称PH曲线。

定理1 [6] ：设 $u\left(t\right),v\left(t\right),w\left(t\right)$ 为非常数的实多项式且 $\left(u\left(t\right),v\left(t\right)\right)=1$，则导数具有如下形式：

$x^{\prime }\left(t\right)=w\left(t\right)\left(u^2\left(t\right)-v^2\left(t\right)\right),y^{\prime }\left(t\right)=2w\left(t\right)u\left(t\right)v(t)$

的平面参数曲线 $q\left(t\right)=\left(x\left(t\right),y\left(t\right)\right)$ 为PH曲线。

2.3. DP-PH曲线几何特征条件

下面推导DP曲线成为PH曲线的条件。

给定两个线性多项式 $a\left(t\right)$ 和 $b\left(t\right)$ ：

$a\left(t\right)=a_0{D}_0^1\left(t\right)+a{D}_1^1(t)$

$b\left(t\right)=b_0{D}_0^1\left(t\right)+b{D}_1^1(t)$

假设 $a_0:a_1$ 和 $b_0:b_1$ 的比率不等，利用三次DP基函数 $\left[D_{0,3}\left(t\right),D_{1,3}\left(t\right),D_{2,3}\left(t\right),D_{3,3}\left(t\right)\right]=\left[{\left(1-t\right)}^3,t\left(1-t\right)\left(2-t\right),t\left(1-t\right)\left(1+t\right),t^3\right]$，求得：

$a{\left(t\right)}^2-b{\left(t\right)}^2=2\left(a_0^2-b_0^2\right)D_0^2\left(t\right)+3\left(a_0{a}_1-b_0{b}_1\right)D_1^2\left(t\right)+2\left(a_1^2-b_1^2\right)D_2^2(t)$

$2a\left(t\right)b\left(t\right)=4a_0{b}_0{D}_0^2\left(t\right)+3\left(a_0{b}_1+a_1{b}_0\right)D_1^2\left(t\right)+4a_1{b}_1{D}_2^2(t)$

$x\left(t\right)={\displaystyle {\int }_0^{t}a{\left(t\right)}^2-b{\left(t\right)}^2\text{d}t}={\displaystyle \underset{k=0}{\overset{3}{\sum }}{x}_k}{D}_k^3(t)$

$y\left(t\right)={\displaystyle {\int }_0^{t}2a\left(t\right)b\left(t\right)\text{d}t}={\displaystyle \underset{k=0}{\overset{3}{\sum }}{y}_k}{D}_k^3(t)$

所以PH-DP曲线的控制点为 $p_k=\left(x_k,y_k\right)$

$\{\begin{array}{l}{p}_1=p_0+\frac{1}{3}\left(2\left(a_0^2-b_0^2\right),4a_0{b}_0\right)\\ p_2=p_1+\frac{1}{3}\left(3\left(a_0{a}_1-b_0{b}_1\right),3\left(a_0{a}_1+b_0{b}_1\right)\right)\\ p_3=p_2+\frac{1}{3}\left(2\left(a_1^2-b_1^2\right),4a_1{b}_1\right)\end{array}$

其中 $p_0$ 是任意给定的。进一步的，给出上述条件公式(2.1)的几何解释，即DP-PH曲线几何特征条件。

定理2：对于任意给定的一个三次DP曲线 $q\left(t\right)$，其控制顶点 $p_k\left(k=0,1,2,3\right)$，控制多边形的各边长为 $L_1,L_2,L_3$， $d_{j,k}$ 表示 $p_k$ 和之间的距离， $\left(j\ne k\right)$， $L_1=d_{01},L_2=d_{12},L_3=d_{23},L=d_{02}$ 顶点 $p_1,p_2$ 对应的角为 ${\theta }_1,{\theta }_2$ 若满足条件

$2L=3\sqrt{L_1{L}_3}$ 和 ${\theta }_1={\theta }_2$

$q\left(t\right)$ 就是PH-DP曲线。

证明

若 $q\left(t\right)$ 满足PH曲线的条件公式(2.1)，可以得到：

$\{\begin{array}{l}{d}_{01}=\frac{2\left(a_0^2+b_0^2\right)}{3}\\ d_{12}=\sqrt{2\left(a_0^2{a}_1^2+b_0^2{b}_1^2\right)}\\ d_{23}=\frac{2\left(a_1^2+b_1^2\right)}{3}\end{array}$

Figure 1. The geometric form of the cubic DP-PH curve

图1. 三次DP-PH曲线控制多边形

如图1可以得出：

$\cos {\theta }_1=\frac{d_{01}^2+d_{12}^2-d_{02}^2}{2d_{12}{d}_{23}}$

$\cos {\theta }_2=\frac{d_{12}^2+d_{23}^2-d_{13}^2}{2d_{12}{d}_{23}}$

根据公式(2.1)给出 $d_{02},d_{13}$，则

$d_{02}^2=\frac{4}{9}\left(a_0^2+b_0^2\right)\left[9{\left(a_0+a_1\right)}^2+9{\left(b_0+b_1\right)}^2\right]$

$d_{13}^2=\frac{4}{9}\left(a_1^2+b_1^2\right)\left[9{\left(a_0+a_1\right)}^2+9{\left(b_0+b_1\right)}^2\right]$

得

$\cos {\theta }_1=\cos {\theta }_2=\frac{-4\left(a_0{a}_1+b_0{b}_1\right)}{3\sqrt{a_0^2{a}_1^2+b_0^2{b}_1^2}}$

可以推得 ${\theta }_1={\theta }_2$ 或 ${\theta }_2=2\pi -{\theta }_1$，进而比较与 $\sin {\theta }_1$ 和 $\sin {\theta }_2$，计算两个角的正弦值，即

$\sin {\theta }_1=\frac{\left(\Delta p_1\times \Delta p_0\right)\cdot z}{d_{12}{d}_{01}}$

$\sin {\theta }_2=\frac{\left(\Delta p_2\times \Delta p_1\right)\cdot z}{d_{23}{d}_{12}}$

其中z是与 $q\left(t\right)$ 平面正交的单位向量，将上述带入，可得到

$\sin {\theta }_1=\sin {\theta }_2=\frac{4\left(a_1{b}_0-a_0{b}_1\right)}{3\sqrt{a_0^2{a}_1^2+b_0^2{b}_1^2}}$

由此得出 ${\theta }_1={\theta }_2$。

并且根据条件公式(2.1)可得控制多边形边长之间的关系满足 $2L=3\sqrt{L_1{L}_3}$。

证毕。

3. 三次DP-PH曲线的几何构造方法

通过DP-PH控制多边形几何特征条件，我们给出可以先通过构造控制多边形进而构造出三次DP-PH曲线的几何构造方法。

给定始末控制顶点 $P_0,P_3$ 和任意点O，依次连接给定的三点，并且令 $\frac{\Vert O{P}_3\Vert }{\Vert P_{0}O\Vert }=\rho \ge 1$，角的范围为 $\pi >\angle P_{0}O{P}_3=\theta >0$。在边 $P_{0}O$ 和 $O{P}_3$ 分别取点 $P_1,P_2$。设单位向量 $e_1=O-P_0$， $e_2=P_3-O$。得到 $P_1,P_2$ 坐标为：

$\{\begin{array}{l}{P}_1=\lambda e_1+P_0\\ P_2=P_3-\mu e_2\end{array}$ (3.1)

$0\le \lambda ,\mu \le 1$。

以点O为原点， $P_{0}O$ 为x轴建立直角坐标系。令 $P_1=\left(1,0\right)$， $e_1=\left(-1,0\right)$， $e_2=\rho \left(\cos \theta ,\sin \theta \right)$，代入式(3.1)得到四个控制顶点的坐标为

$\{\begin{array}{l}{P}_0=\left(1,0\right)\\ P_1=\left(1-\lambda ,0\right)\\ P_2=\left(1-\mu \right)\rho \left(\cos \theta ,\sin \theta \right)\\ P_3=\rho \left(\cos \theta ,\sin \theta \right)\end{array}$ (3.2)

通过定理2中边 $P_{0}O$ 和 $O{P}_3$ 关系之比和推导出的控制多边形几何特征条件可以得出DP曲线控制顶点的中间两点 $D_1,D_2$，并得到如下关系式：

$\angle D_2{D}_1{P}_0=\angle P_3{D}_2{D}_1$ 和 $4{\left|D_1{D}_2\right|}^2=9\left|D_1{P}_0\right|\ast \left|P_3{D}_2\right|，$

由 $\frac{\Vert O{P}_3\Vert }{\Vert P_{0}O\Vert }=\rho \ge 1$ 得到 $\lambda ,\mu$ 表达关系式如下：

$1-\lambda =\rho \left(1-\mu \right)$ (3.3)

$\{\begin{array}{l}\left|D_1{P}_0\right|=\lambda \\ \left|P_3{D}_2\right|=\rho +\lambda +1\\ {\left|D_1{D}_2\right|}^2=2{\left(1-\lambda \right)}^2\left(1-\cos \theta \right)\end{array}$ (3.4)

由式(4)，得到关于参数 $\lambda$ 的方程

$\left(1+8\cos \theta \right){\lambda }^2+\left(7+9\rho -16\cos \theta \right)\lambda +\left(8\cos \theta -8\right)=0$ (3.5)

式(3.4)存在解

$\lambda =\frac{\left(16\cos \theta -7-9\rho \right)\pm \sqrt{81{\left(1+\rho \right)}^2-36\rho \left(8\cos \theta +1\right)}}{2\left(1+8\cos \theta \right)}$ (3.6)

令 $\theta =\frac{\pi }{3}$， $\rho =1$ 由式(3.3)得到 $\lambda =\mu =\frac{2}{5}$，把式(3.3)和(3.6)代入式(3.2)中，从而得出 $D_1=\left(\frac{3}{5},0\right)$， $D_2=\left(\frac{1}{5},\frac{\sqrt{3}}{5}\right)$，取点 $p_0,D_1,D_2,P_3$，画图如图2和图3：

Figure 2. DP curve

图2. DP曲线

Figure 3. DP-PH curve

图3. DP-PH曲线

4. 误差分析和举例

定义2：设 $P\left(t\right)$， $D\left(t\right)$ 是两条三次曲线，这两条曲线误差定义如下：

$\epsilon =\underset{t\in D}{\max }\sqrt{\langle P\left(t\right)-D\left(t\right),P\left(t\right)-D\left(t\right)\rangle }$ (4.1)

给定控制顶点 $P_i,D_i$

定义三次DP-PH曲线为：

$D\left(t\right)={\displaystyle \underset{i=0}{\overset{3}{\sum }}{D}_i{u}_{i,3}(t)}$

则DP曲线和DP-PH曲线的误差表示如下：

$\epsilon ={\displaystyle {\int }_0^a{\Vert {\displaystyle \underset{i=0}{\overset{3}{\sum }}{P}_i{u}_{i,3}\left(t\right)-{\displaystyle \underset{i=0}{\overset{3}{\sum }}{D}_i{u}_{i,3}\left(t\right)}}\Vert }^2\text{d}t}\le a\underset{i=0,1,2}{\max }\Vert P_i-D_i\Vert$

令 $n={\Vert {\displaystyle \underset{i=0}{\overset{3}{\sum }}{P}_i{u}_{i,3}\left(t\right)-{\displaystyle \underset{i=0}{\overset{3}{\sum }}{D}_i{u}_{i,3}\left(t\right)}}\Vert }^2$，则

$n=\left(t^2-2t^3+t^4\right)\left(\frac{{\left(3+\sqrt{3}\right)}^2}{400}-\frac{3}{5}+\frac{{\left(\sqrt{3}-3\right)}^2}{400}{t}^2\right)，$

从而

$\epsilon =\frac{6-3\sqrt{3}}{1400}{a}^7-\frac{2-\sqrt{3}}{200}{a}^6+\frac{27}{250}{a}^5-\frac{126+3\sqrt{3}}{400}{a}^4-\frac{19}{100}{a}^3$

$m=a\underset{t\in D}{\max }=\frac{3+\sqrt{8}}{60}a$

当 $a<0,\epsilon <0$，当 $0<a<1,\epsilon <m$，当 $a>1,\epsilon <m$，则得出两条曲线之间误差较小。

利用定理2中通过控制多边形几何特征条件求出的DP-PH曲线更逼近其控制多边形，具有更好的逼近效果。

基金项目

国家自然科学基金(61702244,61502217)；辽宁省教育厅项目(901132)。