发布时间: 2020-07-16
摘要点击次数:
全文下载次数:
DOI: 10.11834/jig.190534
2020 | Volume 25 | Number 7

计算机图形学

插值给定数据点的四次PH曲线构造

方林聪¹, 阳诚砖¹, 邸文钰², 刘芳¹

1. 浙江财经大学信息管理与人工智能学院, 杭州 310018;

2. 同济大学经济与管理学院, 上海 200082

收稿日期: 2019-11-04; 修回日期: 2019-12-27; 预印本日期: 2020-01-03

基金项目: 浙江省自然科学基金项目（LY18F020023，LQ19F020003，LY19F020011）

第一作者简介: 方林聪, 1982年生, 男, 副教授, 硕士生导师, 主要研究方向为计算机辅助设计与图形学、几何设计与计算。E-mail:lincongfang@gmail.com;
阳诚砖, 男, 讲师, 主要研究方向为计算机视觉、模式识别、人工智能。E-mail:chengzhuanyang@zufe.edu.cn;
邸文钰, 女, 本科生, 主要研究方向为数据可视化、数据挖掘与可视分析。E-mail:1753566@tongji.edu.cn;
刘芳, 女, 讲师, 主要研究方向为数据可视化、数据挖掘与可视分析。E-mail:maggie_liufang@126.com.

中图法分类号: TP391.41

文献标识码: A

文章编号: 1006-8961(2020)07-1473-08

摘要

目的 PH（Pythagorean hodograph）曲线由于具备有理等距曲线、弧长可精确计算等优良的几何性质，广泛应用于数控加工和路径规划等方面。曲线插值是曲线构造的主要手段之一，虽然对PH曲线的Hermite插值方法进行了广泛研究，但插值给定数据点的构造方法仍有待突破，为推广四次PH曲线的应用范围，提出了一种新的四次PH曲线的3点插值问题解决方法。方法从四次PH曲线的代数充分必要条件出发，在该曲线的Bézier控制多边形中引入辅助控制顶点，指出其中实参数的几何意义，该实参数可作为形状调节因子对构造曲线进行交互。对给定的3个平面型值点进行参数化确定相应的参数值；通过对四次PH曲线一阶导数积分得到曲线的显式表达，其中包含一个待定复常量，将给定的约束点代入曲线的显式表达式得到关于待定复常量的一元二次复方程，求解该复方程并反求Bézier控制顶点得到符合约束条件的四次PH曲线。结果实验对通过构造插值给定数据点的四次PH曲线进行比较，当形状调节因此改变时，曲线形状可进行有效交互。每次交互得到两条四次PH曲线，通过弧长、弯曲能量、绝对旋转数的计算得到最优曲线，并构造得到PH曲线的等距线。结论本文方法给定的形状调节参数具有明确的代数意义和几何意义，本文方法易于实现，可有效进行交互。

关键词

计算机辅助设计; Bézier曲线; 控制多边形; 等距曲线; 四次PH曲线; 插值

Construction of quartic PH curves via interpolating given points

Fang Lincong¹, Yang Chengzhuan¹, Di Wenyu², Liu Fang¹

1. School of Information Management and Artificial Intelligence, Zhejiang University of Finance and Economics, Hangzhou 310018, China;

2. School of Economics and Management, Tongji University, Shanghai 200082, China

Supported by: Natural Science Foundation of Zhejiang Province, China (LY18F020023, LQ19F020003, LY19F020011)

Abstract

Objective The problem of interpolating three distinct planar points using quartic Pythagorean hodograph (PH) curves is studied. PH curves comprise an important class of polynomial curves that form a mathematical foundation of most current computer-aided design (CAD) tools. By incorporating special algebraic structures into their hodograph curves, PH curves exhibit many advanced properties over ordinary polynomial parametric curves. These properties include polynomial arc-length functions and rational offset curves. Thus, PH curves are considered a sophisticated solution for a variety of difficult issues arising in applications (e.g., tool paths) in the fields of computer numerically controlled machining and real-time motion control. For example, the arc-length of a PH curve can be computed without numerical integration, accelerating algorithms for numerically controlled machining. The offsets of a PH curve can also be represented exactly rather than being approximated in CAD systems. Thus, analyzing and manipulating PH curves are of considerable practical value in CAD and other applications. PH curves can be represented as widely used Bézier curves. The most intuitive and efficient method for constructing PH curves is by adjusting the control points of Bézier curves under conditions that guarantee PH properties. Therefore, a variety of methods for identifying PH curves are developed. Another important application of PH curves is the geometric construction of these curves. Considerable work has been conducted on Hermite interpolation with different degrees of PH curves. However, methods for interpolating three or more planar points have been rarely studied. Method The necessary and sufficient condition for a planar curve to be a PH curve is a form of a product of complex polynomials, and a Bézier curve and its first derivative are Bernstein polynomials, which are a form of the sum of Bernstein basis functions. We derive a system of complex nonlinear equations by considering the compatibility of the two forms. The geometric meanings of the coefficients are then introduced by presenting several auxiliary points for their Bézier control polygons. To construct a quartic PH curve that interpolates any given three planar points, the first and last points are used as the two endpoints of a Bézier curve. The second point is parameterized by computing the chord lengths by connecting three given points. A complex unknown should be solved considering the integration of the first derivative. The compatibility of complex systems provides a quadratic complex equation with a real coefficient. Thus, in accordance with the fundamental theorem of algebra, two quartic PH curves satisfy any given conditions. A user may interactively construct a series of quartic PH curves by specifying a real coefficient. Result The method is implemented using MATLAB. A maximum of two families of quartic PH curves can be constructed for any given three points. Moreover, arc-lengths, bending energy, and absolute rotation numbers can be computed to select the best solution. Curves with low energy and/or an absolute rotation number can be generally regarded as the best solution because curves with a large bending energy and/or absolute rotation numbers are typically self-intersected. Examples show that the shape can be interactively adjusted by changing a real coefficient, determining the parameter value of the cusp. Lastly, the offsets of the constructed quartic PH curves are shown. Conclusion The proposed method can efficiently construct quartic PH curves for any given three planar points. Only a quadratic complex equation is required to be solved. Thus, the method is robust and efficient. Future studies may consider other applications of the proposed method, e.g., data interpolation using quartic PH splines.

Key words

computer aided design (CAD); Bézier curve; control polygon; offset curve; quartic PH curve; interpolation

0 引言

等距曲线在路桥设计、加工路径规划、数控机器设计等计算机辅助设计(computer aided design，CAD)中有着广泛应用。CAD系统中主要的造型工具是多项式曲线曲面，具有Bézier控制多边形功能，给设计交互造型带来了极大便利。此外，Bézier控制多边形具有凸包性、几何不变性、变差缩减等优良的几何性质，进一步为多项式曲线曲面的几何计算提供了高效算法。

然而，多项式曲线的单位法向量在计算时包含求根运算，导致等距曲线通常不具有有理形式，为此研究者提出了许多等距线的逼近算法。当参数曲线的一阶导数的欧氏范数为一个实多项式(满足Pythagorean条件)时，该曲线的等距线自然可以表示为有理多项式形式，这类曲线称为PH(Pythagorean hodograph)曲线(Farouki和Sakkalis，1990)。由于这种特殊性质，使得PH曲线的弧长、曲率等几何量可被精确计算，而无须使用数值积分运算。同时由于PH曲线及其等距曲线能够很好地与现有CAD系统兼容，等距曲线在CAD系统中也可精确表示。这些性质不仅提升了许多几何计算中的算法效率(Farouki和Neff，1995)，而且为各种实际应用提供了良好的解决方案，特别是传统数控加工中的切、削、磨的刀具路径规划等。

PH曲线是特殊的多项式参数曲线，给定一个多项式曲线及其Bézier控制多边形，能否有效判断该曲线为PH曲线，在逆向工程中具有重要作用(Farouki等，2015)。研究者对PH曲线的判别方法进行了深入研究。早在PH曲线提出时，Farouki和Sakkalis(1990)就给出了三次PH曲线的几何判别方法，指出若三次PH曲线的控制多边形满足两内角相等，且第2条边是首末两边的等比中项，则该曲线为三次PH曲线。Wang和Fang(2009)采用平面曲线的复数表示方法，给出了四次PH曲线的几何判别方法，继而推广到五次、六次、七次PH曲线的几何特征研究(Fang和Wang，2018；Farouki，2008；Zheng等，2016)。此外，Farouki等人(2015)从代数学的角度提出了一种PH曲线的代数判别方法，并给出了在数控实际应用中的可靠计算精度范围。

另一方面，PH曲线的构造方法也在不断改进，为使用者对PH曲线进行交互设计提供了有效工具，困扰PH曲线的Hermite插值问题也有了许多解决办法。Farouki和Neff(1995)针对五次PH曲线的$C^{1}$插值问题，指出对于任意给定的$C^{1}$端点条件，总存在4个可行解，即可以构造4条符合条件的五次PH曲线。Wang等人(2017)通过对4个可行解中的曲线形状的讨论，提出了选取最优解的标准与方法。雍俊海和郑文(2005)从几何角度给出了五次PH曲线的$C^{1}$ Hermite插值问题的一种求解方法。Meek和Walton(1997)针对给定$G^{1}$连续的端点约束条件则存在两条符合条件的三次PH曲线给出了相应的构造方法，Byrtus和Bastl(2010)在此基础上深入分析了插值解的数量与质量。Jüttler(2001)针对七次PH曲线则可满足更高连续阶的端点约束条件，提出了一类正则七次PH曲线在$G^{2}$[$C^{1}$]端点条件下的构造方法，并应用七次PH曲线对任意给定曲线进行逼近。Kong等人(2008)对于特殊的PH曲线的$C^{1}$ Hermite插值方法，采用曲线的复数表示，对PH曲线的导矢曲线零点在复数域中分布情况进行了研究。Farouki(2014)考虑到在数控加工中经常出现的直角过渡曲线设计问题，研究了奇数次PH曲线在$G^{2}$连续的端点约束条件下的构造方法。此外，偶数次的PH曲线也用于曲线构造。Wang和Fang(2009)提出了四次PH曲线$G^{1}$条件下的几何构造方法。郭宇等人(2019)基于Möbius变换引入参数，利用复分析的方法构造了四次有理抛物-PH曲线的$C^{2}$ Hermite插值。方林聪和汪国昭(2014)以及王慧等人(2016)分别讨论了六次PH曲线在$C^{1}$和$G^{2}$端点条件下的构造方法以及圆弧的逼近问题。随着研究的不断深入，PH曲线的相关理论推广到了高维曲线(Farouki等，2002；Sakkalis和Farouki，2012)、有理曲线(Pottmann，1995；Farouki和Šír，2011)以及Minkowski空间(Moon，1999；Kosinka和Jüttler，2006)。更多的内容以及实例可以参考Farouki等人(2008)以及Kosinka和Lávićka(2014)的相关论述。总之，PH曲线Hermite插值方法的研究成果比较丰富，然而PH样条的构造方法尚未深入研究。

数据插值在几何造型中的应用非常广泛，设计者希望通过给定的型值点构造相应的几何模型。本文以PH曲线的代数充要条件为基础，提出了四次PH曲线的插值构造方法，并讨论了均匀参数化、弦长参数化、弧长参数化方法在插值效果上的差异。通过比较曲线的弧长、弯曲能量、绝对旋转数等几何量来选择最优曲线。

1 四次PH曲线

定理1 一条平面参数曲线是PH曲线当且仅当其一阶导数有如下形式的因式分解

$ {P^\prime }(t) = w(t){[u(t) + {\rm{i}}v(t)]^2} $

式中，$u(t)$，$v(t)$，$w(t)$为实多项式，且$u(t)$，$v(t)$互素。

若四次平面参数曲线$P (t)$为PH曲线，则有$u(t)$，$v(t)$，$w(t)$均为一次实多项式。因此四次PH曲线的一阶导数有如下形式

$ {P^\prime }(t) = [a(1 - t) + t]{[{z_0}(1 - t) + {z_1}t]^2} $

(1)

式中，$a∈ {\bf R} $，$z_{0}$，$z_{1}∈ C $, $z_{0}≠kz_{1}$，$∀k∈ {\bf R} $。注意到$P'\left( {\frac{a}{{a - 1}}} \right) = 0$，因此实系数$a$确定了四次PH曲线的奇异点的参数值。

令$P _{i}$，$i=0, …, n$，为Bézier曲线$P (t)$的控制顶点，则Bézier曲线的一阶导数可表示为

$ {P^\prime }(t) = 3\sum\limits_{i = 0}^2 {B_i^2} (t)\Delta {P_i} $

(2)

式中，Δ$P _{i}=P _{i+1}－P _{i}$，$i=0, 1, 2$。通过比较Bernstein多项式系数，得

$ \left\{ \begin{array}{l} 4\Delta {P_0} = az_0^2\\ 12\Delta {P_1} = z_0^2 + 2a{z_0}{z_1}\\ 12\Delta {P_2} = az_1^2 + 2{z_0}{z_1}\\ 4\Delta {P_3} = z_1^2 \end{array} \right. $

(3)

不妨令

$ \left\{ \begin{array}{l} {Q_0} = {P_1} + \frac{{z_0^2}}{{12}}\\ {Q_1} = {P_3} + \frac{{z_1^2}}{{12}} \end{array} \right. $

于是有

$ {a = \frac{{\Delta {P_0}}}{{3({Q_0} - {P_1})}} = \frac{{{P_2} - {Q_0}}}{{{Q_1} - {P_2}}} = \frac{{3({P_3} - {Q_1})}}{{\Delta {P_3}}}} $

$ {{\rm{arg}}\frac{{{P_2} - {Q_0}}}{{\Delta {P_0}}} = {\rm{arg}}\frac{{\Delta {P_3}}}{{{Q_1} - {P_2}}} = {\rm{arg}}\frac{{{z_1}}}{{{z_0}}}} $

即直线$Q _{0}Q _{1}$为直线$P _{0}P _{1}$与$P _{3}P _{4}$所成角的角平分线，且经过控制顶点$P _{2}$，因此$a$的几何意义是明确的，如图 1所示。

图 1 四次PH曲线的Bézier控制多边形

Fig. 1 Bézier control polygon of a quartic PH curve

给定一条平面参数曲线，令$k$为曲线的曲率，则该曲线的弧长$S$、弯曲能量$E$、绝对旋转数$R_{\rm abs}$等曲线度量计算方法如下

$ {S = \int {\rm{d}} s = \int_0^1 {\sqrt {{x^{\prime 2}}(t) + {y^{\prime 2}}(t)} } {\rm{d}}t} $

$ {E = \int {{k^2}} {\rm{d}}s = \int_0^1 {{k^2}} \left\| {{P^\prime }(t)} \right\|{\rm{d}}t} $

$ {{R_{{\rm{abs}}}} = \frac{1}{{2\pi }}\int | k|{\rm{d}}s = \frac{1}{{2\pi }}\int_0^1 | k|{\kern 1pt} {\kern 1pt} {\kern 1pt} \left\| {{P^\prime }(t)} \right\|{\rm{d}}t} $

2 三点插值

考虑给定平面上3个型值点$\tilde P _{i}$, $i=0, 1, 2$，要求构造一条三次PH曲线$P (t)$，$t∈[0, 1]$, 插值这些数据点，其中有$P (0)=\tilde P _{0}=P _{0}$，$P (1)=\tilde P _{2}=P _{4}$。记$\tilde P _{1}$点的参数值为$t_{1}$，即$P (t_{1})=\tilde P _{1}$。这里$t_{1}$可以通过弦长参数化方法计算得到，即

$ {t_1} = \frac{{\left\| {{{\tilde P}_1} - {{\tilde P}_0}} \right\|}}{{\left\| {{{\tilde P}_1} - {{\tilde P}_0}} \right\| + \left\| {{{\tilde P}_2} - {{\tilde P}_1}} \right\|}} $

(4)

由于$P ′(t)$具有式(1)的形式，因此$P (t)$可通过对式(1)积分得到，即

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} P(t) = \int {{P^\prime }} (t){\rm{d}}t = \\ \int {[a(1 - t) + t]{{[{z_0}(1 - t) + {z_1}t]}^2}{\rm{d}}t = } \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rho _0}(t)z_0^2 + {\rho _1}(t){z_0}{z_1} + {\rho _2}(t)z_1^2 + Z \end{array} $

式中，$Z $为待定的复常数，其余实多项式系数$ρ_{0}(t)$，$ρ_{1}(t)$，$ρ_{2}(t)$的表达式为

$ \begin{array}{l} {\rho _0}(t) = - \frac{a}{4}{(1 - t)^4} + \frac{1}{2}{(1 - t)^2}{t^2} + \frac{1}{3}{t^3} - \frac{1}{4}{t^4}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {{\rho _1}(t) = a{{(1 - t)}^2}{t^2} + \frac{2}{3}(1 - t){t^3} - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{2}a{t^4} + \frac{1}{6}{t^4} + \frac{2}{3}a{t^3}}\\ {{\rho _2}(t) = \frac{1}{3}a(1 - t){t^3} + \frac{1}{4}{t^4} + \frac{1}{{12}}a{t^4}} \end{array} \end{array} $

(5)

记$f_{i}=ρ_{i}(t_{1})$，$i=0, 1, 2$。由于$P (0)=\tilde P _{0}=P _{0}$，$P (t_{1})=\tilde P _{1}$，$P (1)=\tilde P _{2}=P _{4}$，可以得到方程组

$ \left\{ {\begin{array}{*{20}{l}} { - \frac{a}{4}z_0^2 + Z = {{\tilde P}_0}}\\ {{f_0}z_0^2 + {f_1}{z_0}{z_1} + {f_2}z_1^2 + Z = {{\tilde P}_1}}\\ {\frac{1}{{12}}z_0^2 + \frac{{1 + a}}{6}{z_0}{z_1} + \frac{{3 + a}}{{12}}z_1^2 + Z = {{\tilde P}_2}} \end{array}} \right. $

在该方程组中，将$z^{2}_{0}$，$z_{0}z_{1}$，$z^{2}_{1}$作为未知数，可求得

$ \left\{ {\begin{array}{*{20}{l}} {z_0^2 = \frac{4}{a}(Z - {{\tilde P}_0})}\\ {{z_0}{z_1} = \frac{{{f_2}({g_4} - {g_3}Z) - {g_1}({g_5} - {g_2}Z)}}{{{f_2}{g_0} - {f_1}{g_1}}}}\\ {z_1^2 = \frac{{{g_0}({g_5} - {g_2}Z) - {f_1}({g_4} - {g_3}Z)}}{{{f_2}{g_0} - {f_1}{g_1}}}} \end{array}} \right. $

(6)

式中，

$ \left\{ {\begin{array}{*{20}{l}} {{g_0} = \frac{{1 + a}}{6};{g_1} = \frac{{3 + a}}{{12}};}\\ {{g_2} = \frac{{4{f_0}}}{a};{g_3} = \frac{1}{{3a}};}\\ {{g_4} = {{\tilde P}_2} + \frac{1}{{3a}}{{\tilde P}_0};{g_5} = {{\tilde P}_1} + \frac{{4{f_0}}}{a}{{\tilde P}_0}} \end{array}} \right. $

(7)

进一步考虑$(z_{0}z_{1})^{2}=z^{2}_{0}z^{2}_{1}$, 得到关于$Z $的一元二次复方程

$ A{Z^2} + BZ + C = 0 $

(8)

式中，

$ \left\{ {\begin{array}{*{20}{l}} {A = {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{({g_1}{g_2} - {f_2}{g_3})}^2} - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{4}{a}({f_2}{g_0} - {f_1}{g_1})({f_1}{g_3} - {g_0}{g_2})}\\ {B = {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 2({g_1}{g_2} - {f_2}{g_3})({f_2}{g_4} - {g_1}{g_5}) - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{4}{a}({f_2}{g_0} - {f_1}{g_1})[{g_0}{g_5} - {f_1}{g_4} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ({f_1}{g_3} - {g_0}{g_2}){{\tilde P}_0}]}\\ {C = {\kern 1pt} {\kern 1pt} {\kern 1pt} {{({f_2}{g_4} - {g_1}{g_5})}^2} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{4}{a}({f_2}{g_0} - {f_1}{g_1})({g_0}{g_5} - {f_1}{g_4}){{\tilde P}_0}} \end{array}} \right. $

(9)

由代数基本定理知，当给定实数$a$时，方程式(8)存在两个复数解，即存在两条插值四次PH曲线。综上所述，给定实数平面上任意3个型值点$\tilde P _{i}$，$i=0, 1, 2$，以及实数$a$(实数$a$具有明确的代数意义和几何意义)，则插值的四次PH曲线可以用如下算法构造得到。

输入：复数表示的3个型值点$\tilde P _{i}$，$i=0, 1, 2$，一个实参数$a$。

输出：复数表示的控制顶点$P _{i}$，$i=0, …, 4$。

1) 根据式(4)计算$\tilde P _{1}$的参数值$t_{1}$。

2) 根据式(5)计算$f_{i}$，$i=0, 1, 2$。

3) 根据式(7)计算$g_{i}$，$i=0, …, 5$。

4) 将$f_{i}$，$g_{j}$，$i=0, 1, 2$，$j=0, …, 5$代入式(9)得到$A $，$B $，$C $。

5) 求解方程式(8)得到$Z $的两个解。

6) 通过式(6)计算得到$z^{2}_{0}$，$z_{0}z_{1}$，$z^{2}_{1}$。

7) 已知$P _{0}=\tilde P_{0}$，$P _{4}=\tilde P_{2}$，按照式(3)计算得$P _{i}$，$i=1, 2, 3$。

3 数值实例

图 2—图 4是应用本文算法构造的两个四次PH曲线实例。在图 2的实例中，给定型值点$\tilde P _{0}$=0，$\tilde P _{1}=3.5+2$i$，$\tilde P _{2}=6$，给出了$a$=1, …, 10时构造的四次PH曲线。为比较得到的两族曲线，表 1给出了相应的曲线几何量，包括曲线的弧长、弯曲能量、绝对旋转数等。注意到对于相同的参数$a$，得到的插值曲线具有相同的弧长，但图 2(a)中的曲线具有较小的弯曲能量和旋转数，曲线无自交；图 2(b)中的曲线具有较大的弯曲能量和旋转数，表现为自交曲线，实际应用中通常不采用。图 3给出了当$a=2$时，本文方法构造的曲线及其等距曲线，图中实线表示PH曲线，虚线表示等距曲线。

图 2 第1个实例的两族四次PH曲线

Fig. 2 The first example of two families of quartic PH curves constructed ((a) curves with smaller bending energy and absolute rotation numbers are not self-intersected; (b) curves with bigger bending energy and absolute rotation numbers are self-intersected)

图 3 $a=2$时第1个实例的两条四次PH曲线及其等距曲线

Fig. 3 Two quartic PH curves with their offsets of the first example when $a$ = 2 ((a) none self-intersected curve with its offsets; (b) self-intersected curve with its offsets)

图 4 第2个实例的两族四次PH曲线

Fig. 4 The second example of two families of quartic PH curves constructed ((a) curves with smaller bending energy and absolute rotation numbers are not self-intersected; (b) curves with bigger bending energy and absolute rotation numbers are self-intersected)

表 1 图 2中两族曲线几何量比较
Table 1 Comparison of curves' geometric measures in Fig. 2

下载CSV

a	图 2(a)			图 2(b)
a	弧长	弯曲能量	旋转数	弧长	弯曲能量	旋转数
1	7.617 7	0.630 5	0.342 8	7.617 7	57.839 6	0.874 0
2	7.683 2	0.654 2	0.350 8	7.683 2	53.881 9	0.873 0
3	7.746 9	0.692 4	0.360 6	7.746 9	50.319 1	0.872 5
4	7.801 1	0.729 2	0.369 2	7.801 1	47.660 8	0.872 2
5	7.846 6	0.762 0	0.376 3	7.846 6	45.663 5	0.872 0
6	7.884 9	0.790 7	0.382 3	7.884 9	44.122 0	0.871 8
7	7.917 6	0.815 6	0.387 3	7.917 6	42.900 4	0.871 6
8	7.945 8	0.837 4	0.391 5	7.945 8	41.909 8	0.871 5
9	7.970 2	0.856 5	0.395 2	7.970 2	41.091 1	0.871 4
10	7.991 7	0.873 4	0.398 3	7.991 7	40.403 2	0.871 2

在图 4的实例中，给定型值点$\tilde P _{0}=0$，$\tilde P _{1}=3.5+10$i$，$\tilde P _{2}=6$，给出了$a=1, …, 10$时构造的四次PH曲线。为比较得到的两族曲线，表 2给出了相应的曲线几何量，包括曲线的弧长、弯曲能量、旋转数等。注意到对于相同的参数$a$，得到的插值曲线具有相同的弧长，但图 4(a)中的曲线具有较小的弯曲能量和旋转数，曲线无自交；图 4(b)中的曲线具有较大的弯曲能量和旋转数，表现为自交曲线，通常不被实际应用所采用。图 5给出了当$a=2$时，本文方法构造的曲线及其等距曲线，图中PH曲线是实线表示的，而等距曲线是虚线表示的。

图 5 $a=2$时第2个实例的两条四次PH曲线及其等距曲线

Fig. 5 Two quartic PH curves with their offsets of the second example where $a$ = 2 ((a) none self-intersected curve with its offsets; (b) self-intersected curve with its offsets)

表 2 图 4中两族曲线几何量比较
Table 2 Comparison of curves' geometric measures in Fig. 4

下载CSV

a	图 4(a)			图 4(b)
a	弧长	弯曲能量	旋转数	弧长	弯曲能量	旋转数
1	23.888 3	0.819 9	0.591 9	23.888 3	2.348 1	0.732 5
2	24.111 9	0.812 9	0.595 5	24.111 9	2.289 5	0.734 0
3	24.397 9	0.805 9	0.600 4	24.397 9	2.221 4	0.736 2
4	24.647 9	0.800 9	0.604 5	24.647 9	2.168 4	0.738 0
5	24.857 1	0.797 3	0.607 8	24.857 1	2.127 9	0.739 4
6	25.031 9	0.794 6	0.610 5	25.031 9	2.096 5	0.740 6
7	25.179 3	0.792 6	0.612 7	25.179 3	2.071 5	0.741 6
8	25.304 8	0.791 1	0.614 6	25.304 8	2.051 2	0.742 4
9	25.412 7	0.789 9	0.616 1	25.412 7	2.034 5	0.743 0
10	25.506 5	0.788 9	0.617 5	25.506 5	2.020 4	0.743 6

4 结论

本文指出给定平面3个型值点的条件下，依次插值这3个点的四次PH曲线的构造问题可转化为一个带有实参数的一元二次复方程的求解问题。因此，由代数基本定理可知，当实参数给定时，存在不超过两条四次PH曲线满足插值条件。本文通过对四次PH曲线代数充分必要条件的深入分析，给出了一元二次复方程的具体形式，进一步提出了详细的曲线构造算法。

实验结果表明，本文方法可有效构造四次PH曲线，且实参数可作为形状调节因子对曲线进行交互构造。通过对弧长、弯曲能量、绝对旋转数的计算，从两个可行解中可有效排除自交曲线。因此，本文的方法易于实现，应用灵活。

然而，由于本文方法未考虑端点的切向及高阶导数的约束条件，因此无法直接用于构造四次PH样条曲线。本文工作继续推广，可对任意给定的曲线使用四次PH曲线进行插值逼近，今后也可进一步讨论PH样条构造的可能性。

参考文献

Byrtus M, Bastl B. 2010. G¹ Hermite interpolation by PH cubics revisited. Computer Aided Geometric Design, 27(8): 622-630 [DOI:10.1016/j.cagd.2010.06.004]

Fang L C, Wang G Z. 2014. C¹ Hermite interpolation using sextic PH curves. Scientia Sinica Mathematica, 44(7): 799-804 (方林聪, 汪国昭. 2014. 六次PH曲线C¹ Hermite插值. 中国科学:数学, 44(7): 799-804) [DOI:10.1360/N012013-00124]

Fang L C, Wang G Z. 2018. Geometric characteristics of planar quintic Pythagorean-hodograph curves. Journal of Computational and Applied Mathematics, 330: 117-127 [DOI:10.1016/j.cam.2017.08.014]

Farouki R T. 2008. Pythagorean-Hodograph Curves: Algebra and Geometry Inseparable. Berlin: Springer

Farouki R T. 2014. Construction of G² rounded corners with Pythagorean-hodograph curves. Computer Aided Geometric Design, 31(2): 127-139 [DOI:10.1016/j.cagd.2014.02.002]

Farouki R T, Al-Kandari M, Sakkalis T. 2002. Hermite interpolation by rotation-invariant spatial Pythagorean-hodograph curves. Advances in Computational Mathematics, 17(4): 369-383 [DOI:10.1023/A:1016280811626]

Farouki R T, Giannelli C, Sestini A. 2015. Identification and "reverse engineering" of Pythagorean-hodograph curves. Computer Aided Geometric Design, 34: 21-36 [DOI:10.1016/j.cagd.2015.04.001]

Farouki R T, Jüttler B, Manni C. 2008. Pythagorean-hodograph curves and related topics. Computer Aided Geometric Design, 25(4/5): 203-204 [DOI:10.1016/j.cagd.2008.01.001]

Farouki R T, Neff C A. 1995. Hermite interpolation by Pythagorean hodograph quintics. Mathematics of Computation, 64(212): 1589-1609 [DOI:10.1090/S0025-5718-1995-1308452-6]

Farouki R T, Sakkalis T. 1990. Pythagorean hodographs. IBM Journal of Research and Development, 34(5): 736-752 [DOI:10.1147/rd.345.0736]

Farouki R T, Šír Z. 2011. Rational Pythagorean-hodograph space curves. Computer Aided Geometric Design, 28(2): 75-88 [DOI:10.1016/j.cagd.2011.01.002]

Guo Y, Jiang P, Wang J M, Liu Z. 2019. C² Hermite interpolation based on quartic rational parabolic-PH curves by using Möbius transformation. Journal of Image and Graphics, 24(1): 96-102 (郭宇, 江平, 王剑敏, 刘植. 2019. Möbius变换下四次有理抛物-PH曲线的C²Hermite插值. 中国图象图形学报, 24(1): 96-102) [DOI:10.11834/jig.180318]

Jüttler B. 2001. Hermite interpolation by Pythagorean hodograph curves of degree seven. Mathematics of Computation, 70(235): 1089-1111 [DOI:10.1090/S0025-5718-00-01288-6]

Kong J H, Jeong S P, Lee S, Kim G I. 2008. C¹ Hermite interpolation with simple planar PH curves by speed reparametrization. Computer Aided Geometric Design, 25(4/5): 214-229 [DOI:10.1016/j.cagd.2007.11.006]

Kosinka J, Jüttler B. 2006. C¹Hermite interpolation by Minkowski Pythagorean hodograph cubics. Computer Aided Geometric Design, 23(5): 401-418 [DOI:10.1016/j.cagd.2006.01.004]

Kosinka J, Lávička M. 2014. Pythagorean hodograph curves:a survey of recent advances. Journal for Geometry and Graphics, 18(1): 23-43

Meek D S, Walton D J. 1997. Geometric Hermite interpolation with Tschirnhausen cubics. Journal of Computational and Applied Mathematics, 81(2): 299-309 [DOI:10.1016/S0377-0427(97)00066-6]

Moon H P. 1999. Minkowski Pythagorean hodographs. Computer Aided Geometric Design, 16(8): 739-753 [DOI:10.1016/S0167-8396(99)00016-3]

Pottmann H. 1995. Rational curves and surfaces with rational offsets. Computer Aided Geometric Design, 12(2): 175-192 [DOI:10.1016/0167-8396(94)00008-G]

Sakkalis T, Farouki R T. 2012. Pythagorean-hodograph curves in Euclidean spaces of dimension greater than 3. Journal of Computational and Applied Mathematics, 236(17): 4375-4382 [DOI:10.1016/j.cam.2012.04.002]

Wang G Z, Fang L C. 2009. On control polygons of quartic Pythagorean-hodograph curves. Computer Aided Geometric Design, 26(9): 1006-1015 [DOI:10.1016/j.cagd.2009.08.003]

Wang H, Zhu C G, Li C Y. 2016. G² Hermite interpolation by Pythagorean hodograph of degree six. Journal of Graphics, 37(2): 155-165 (王慧, 朱春钢, 李彩云. 2016. 六次PH曲线G² Hermite插值. 图学学报, 37(2): 155-165) [DOI:10.11996/JG.j.2095-302X.2016020155]

Wang H, Zhu C G, Li C Y. 2017. Identification of planar sextic Pythagorean-hodograph curves. Journal of Mathematical Research with Applications, 37(1): 59-72

Yong J H, Zheng W. 2005. Geometric method for Hermite interpolation by a class of PH quintics. Journal of Computer-Aided Design and Computer Graphics, 17(5): 990-995 (雍俊海, 郑文. 2005. 一类五次PH曲线Hermite插值的几何方法. 计算机辅助设计与图形学学报, 17(5): 990-995) [DOI:10.3321/j.issn:1003-9775.2005.05.019]

Zheng Z H, Wang G Z, Yang P. 2016. On control polygons of Pythagorean hodograph septic curves. Journal of Computational and Applied Mathematics, 296: 212-227 [DOI:10.1016/j.cam.2015.09.006]