CN102609987B - 一种通过计算零维三角多项式系统所有实根及其重数进行曲面绘制的方法和系统 - Google Patents

Abstract

一种通过计算零维三角多项式系统所有实根及其重数进行曲面绘制的方法和系统，包括步骤：确定曲面特征点，通过计算零维三角多项式系统所有实根及其重数得到曲面的特征点及其重数；生成曲面连接关系图，通过确定特征点之间的连接关系得到连接关系图；绘制曲面，在连接关系图的基础上直接绘制网格曲面或者通过其他光顺化方法绘制更光滑的曲面。

Description

一种通过计算零维三角多项式系统所有实根及其重数进行曲面绘制的方法和系统

技术领域

本发明涉及计算机三维建模技术，特别是涉及一种通过计算零维三角多项式系统所有实根及其重数进行曲面绘制的方法和系统。

背景技术

确定给定代数曲面的几何和连接关系，并用网格近似表示曲面是计算机图形学和几何建模的一个重要课题。曲面网格化可以正确地显示曲面，也可以用来展示曲面上的工程应用，例如有限元分析、计算机三维建模、计算机辅助设计、对象三维重建等应用。

如何绘制出给定的隐式方程确定的曲面是近年来颇为流行的一个课题。最流行的就是Marching Cube方法，将给定的空间不断细分为更小的立方体，直到每个小立方体中曲面的几何结构都能够被确定(参见LORENSEN，W.E.，AND CLINE，H.E.1987.Marching cubes：a high resolution 3d surfaceconstruction algorithm.In Proc.SIGGRAPH 1987，ACM Press.).Plantinga和Vegter提出了法线变化条件从而给出更好的网格化算法(参见PLANTINGA，S.，AND VEGTER，G.2004.Isotopic implicit surface meshing.In Proc.Symp.onGeometry Processing，251-260.PLANTINGA，S.，AND VEGTER，G.2007.Isotopic meshing of implicit surfaces.The Visual Computer 23，45-58.).Hart等提出了基于Morse理论的算法(参见HART，J.C.1998.Morse theory for implicitsurface modeling.In Mathematical Visualization，Springer-verlag，H.Hege and K.Polthier，Eds.，257-268.)，该算法引入一个变元来确定曲面的几何结构。这些算法都是纯粹的数值方法，直接对曲面网格化，速度较快，但是都只能处理无奇点的曲线或曲面，对有奇点的几何对象，可能会绘制出错误的图形。

近年有些学者提出符号数值混合算法，先用柱形代数方法将欧氏空间分为柱形胞腔，使给定曲面在每个胞腔中都不变号，从而得到一种计算曲面几何结构的新方法。这类方法的基本步骤都是先通过分解空间确定曲面的几何结构和连接关系，然后再将非奇部分网格化(参见D.S.Arnon，G.Collins，andS.Mccallum，1988.Cylindrical algebraic decomposition，iii：an adjacencyalgorithm for three-dimensional space.J.Symbolic Comput.5，1，2，163-187.J.S.Cheng，X.S.Gao，and M.Li，2005.Determine the topology of real algebraicsurfaces.In Mathematics of Surfaces，Springer-Verlag，121-146.S.Mccallum，and，G.E.Collins，2002.Local box adjacency algorithms for cylindrical algebraicdecompositions.J.Symbolic Comput.33，321-342.B.Mourrain，and J.P.Tecourt，2005.Computing the topology of real algebraic surfaces.In MEGA electronicproceedings.)。在这种方法中，如何进行空间的分解就显得尤为重要，空间分解一般是通过确定曲面的奇点，以及曲面在奇点处的连接关系来得到。曲面的奇点是一个零维三角多项式系统的实根，奇点连接着的分枝数，即该奇点的重数，就是对应实根的重数。

因此，求解零维三角多项式系统，即有有限个零点的三角多项式系统，是计算科学、工程等领域的一个基本问题。用数值的方法来求解零维三角多项式系统一般会因为数值表示的误差导致结果不正确，并且不能够计算重数。Lu等人提出一种隔离零维三角系统实根的方法(参见Z.Lu，B.He，Y.Luo and L.Pan，An algorithm of real root isolation for polynomial systems，in Proceedings ofInternational Workshop on Symbolic-Numeric Computation，Xi′an，China，Julyl9-21，94-107，2005.)，但是他们的方法没有终止判定条件，也不能处理重根。Collins等用区间算法和Descartes方法来计算(参见G.E.Collins，J.R.Johnson，and W.Krandick，Interval arithmetic in cylindrical algebraicdecom-position，J.Symb.Comput.，34：145-157，2002.)。Xia和Yang也提出了一种基于结式计算的算法(参见B.Xia and L.Yang，An algorithm for isolating thereal solutions of semi-algebraic systems，J.Symb.Comput.，34：461-477，2002.)，但是他们的方法在有些情况下会失效。

发明内容

基于现有技术的不足，本发明提供一种通过计算零维三角多项式系统所有实根及其重数进行曲面绘制的方法和系统，该方法通过计算零维三角多项式系统的所有实根及其重数得到曲面的特征点及其重数，从而确定特征点之间连接关系得到曲面连接关系图，最后绘制出空间结构准确的曲面图像。

一种通过计算零维三角多项式系统所有实根及其重数进行曲面绘制的方法，包括以下步骤：

1)确定曲面的特征点；

特征点包括曲面的独立奇点，奇异线上的点，与奇点邻近的非奇点。执行如下步骤：

(1)确定特征点满足的零维三角多项式系统；

(2)求解零维三角多项式系统；包括以下步骤：

I.对多项式系统进行预处理，通过线性变换将其变换到局部一般位置，包括以下步骤：

(a)计算实根的隔离界；

(b)估计根的上界，执行以下步骤：

(i)将多项式拆分成两个单项式项数最少的正系数多项式之和；

(ii)计算上界多项式和下界多下式；

(iii)计算上界多项式和下界多项式的根的上界；

(iv)取上界多项式和下界多项式的根的上界的最大值。

(c)确定线性变换，把多项式系统变换到局部一般位置；

II.将所有实根投影到一维空间中，通过计算结式来进行投影，得到一个单变元方程，该投影保持根的重数不变；

III.求计算单变元多项式方程的所有实根及其重数。

IV.通过线性变换将一元实根提升回原方程组的实根。

2)生成曲面连接关系图1

曲面的连接关系图包括三部分：特征点(点集合)，特征点之间的连接关系(边集合)，特征点构成的三角片与边的连接关系(面集合)。

3)绘制曲面。

步骤2)中得到的面集合已经是曲面的一个三角网格化，我们可以直接用他来绘制网格曲面，或者在此基础上通过其他光顺化方法绘制更光滑的曲面。

相应地，本发明提出了设计了一种通过计算零维三角多项式系统所有实根及其重数进行曲面绘制的系统。

目前，许多需要对象三维建模的领域都涉及到如何绘制隐式曲面图像的问题。本发明提出的技术方案，采用首先通过计算零维三角多项式系统得到曲面的特征点，然后通过确定特征点之间连接关系得到曲面连接关系图，最后通过连接关系图绘制曲面的方法，避免了纯数值方法绘制有奇点隐式曲面出错的工程难题，能够保证绘制出空间结构准确的立体曲面图像。本发明方法明确，结果鲁棒，可以用于隐式曲面的立体绘制。

附图说明

图1本发明的技术方案流程图；

图2确定曲面特征点流程图；

图3计算零维三角多项式系统所有实根流程图；

图4多项式系统预处理流程图；

图5估计根的上界流程图；

图6曲面连接关系图中确定边集合示意图；

图7曲面连接关系图中确定面集合示意图；

图8-10三维绘制效果图；

图11本发明的系统框图。

具体实施方式

一种通过计算零维三角多项式系统所有实根及其重数进行曲面绘制的方法和系统，首先，确定曲面的特征点；然后确定特征点处曲面的连接情况；最后，后绘制出曲面。具体流程参见图1。

下面具体介绍关键的实现细节：

1.确定曲面f(x₁，x₂，x₃)＝0的特征点，流程如图2所示。

特征点包括曲面的独立奇点，奇异线上的点，与奇点邻近的非奇点。执行如下步骤：

(1)确定特征点满足的零维三角多项式系统。

奇点包含在零维三角多项式系统

Σ₃＝{f₁(x₁)，f₂(x₁，x₂)，f₃(x₁，x₂，x₃)}  (1)

的实根中，其中

f₃(x₁，x₂，x₃)＝f(x₁，x₂，x₃)

$f_2(x_1,x_2)=\mathrm{SquareFreePart}\left({\mathrm{Resul}\tan t}_{x3}(f_3,\frac{\∂f_3}{\∂x_3})\right);---\left(2\right)$

$f_1\left(x_1\right)=\mathrm{SquareFreePart}\left({\mathrm{Resul}\tan t}_{x2}(f_2,\frac{\∂f_2}{\∂x_2})\right)$

奇异线上的点包含在零维三角系统

Σ₂＝{f₂(c₁，x₂)，f₃(c₁，x₂，x₃)}  (3)

的实根中；

奇点邻近的非奇点包含在零维三角系统

Σ₁＝{f₃(c₁，c₂，x₃)}  (4)

的实根中。

(2)求解零维三角多项式系统，流程如图3所示。

在计算实根的过程中，我们始终用一个包含该实根，端点为有理数的区间来近似表示实根，其中涉及到实根的计算都使用区间进行计算。对解零维三角多项式系统Σ_n＝{f₁(x₁)，f₂(x₁，x₂)，…，f_n(x₁，x₂，…，x_n)}，执行如下步骤：

I.进行预处理，通过线性变换将其变换到局部一般位置，流程如图4所示；

(a)计算零维三角多项式系统中第一个多项式组成的多项式系统Σ₁＝{f₁(x₁)}的所有实根Zero(Σ₁)；

(b)对i从2到n-1，令Σ_i＝{f₁(x₁)，f₂(x₁，x₂)，…，f_i(x₁，x₂，…，x_i)}，计算Σ_i的所有实根；

执行以下步骤：

(i)计算最倒数第二个方程关于x_i-1的实根的隔离界：

$r_{i-1}=\frac{1}{2}\underset{({\mathrm{\ξ}}_1,{\mathrm{\ξ}}_2,...,{\mathrm{\ξ}}_{i-2})\∈\mathrm{Zero}\left({\mathrm{\Σ}}_{i-2}\right)}{\min }\left\{\right|{\mathrm{\θ}}_1-{\mathrm{\θ}}_2|:{\mathrm{\θ}}_1,{\mathrm{\θ}}_2\∈\mathrm{Zero}\left(f_{i-1}({\mathrm{\ξ}}_1,{\mathrm{\ξ}}_2,...,{\mathrm{\ξ}}_{i-2},x_{i-1})\right),{\mathrm{\θ}}_1\≠{\mathrm{\θ}}_2\}.---\left(5\right)$

(ii)计算最后一个方程关于x_i的根的上界R_i，流程如图5所示：

设实数列(ξ₁，ξ₂，…，ξ_i-1)由一个包含它且端点为有理数的区间列[a₁，b₁]×[a₂，b₂]×…×[a_i-1，b_i-1]表示，将多项式拆分：

f_i(x₁，x₂，…，x_i)＝f_i ⁺-f_i ^-，  (6)

其中f_i ⁺，f_i ^-∈Z⁺[x₁，x₂，…，x_i]，是满足上式的单项式项数最少的正系数多项式。

定义上界多项式f_i ^u和下界多项式f_i ^d(x_i)：

f_i ^u(x_i)＝f_i ⁺(b₁，…，b_i-1，x_i)-f_i ^-(a₁，…，a_i-1，x_i)，

                                                 (7)

f_i ^d(x_i)＝f_i ⁺(a₁，…，a_i-1，x_i)-f_i ^-(b₁，…，b_i-1，x_i)，

设f_i ^u(x_i)形为

${f_i}^u\left(x_i\right)=c_d{x}_i^d+c_{d-1}{x}_i^{d-1}+...+c_0,---\left(8\right)$

计算f_i ^u(x_i)的根的上界

$\mathrm{Ru}=1+\underset{0\≤k\≤d-1}{\max }\left|\frac{c_k}{c_d}\right|,---\left(9\right)$

类似可以计算Rd。则f_i(ξ₁，ξ₂，…，ξ_i-1，x_i)的根的上界为

R＝max{Ru，Rd}。  (10)

最后计算f_i(x₁，x₂，…，x_i)关于x_i的根的上界

R_i＝max{R：R是f_i(ξ₁，ξ₂，…，ξ_i-1，x_i)根的上界，(ξ₁，ξ₂，…，ξ_i-1)∈Zero(Σ_i-1)}。  (11)

(iii)用线性变换将方程组Σ_i＝{f₁(x₁)，f₂(x₁，x₂)，…，f_i(x₁，x₂，…，x_i)}变换到局部一般位置：

对Σ_i做线性变换

得到新的多项式系统

${\mathrm{\Σ}}_i^{\′}=\{f_1(X_i^i-\frac{r_1}{R_2}{X}_{i-1}^i),...,f_i(X_i^i-\frac{r_1}{R_2}{X}_{i-1}^i,...,X_2^i-\frac{r_{i-1}}{R_i}{X}_1^i,X_1^i)\},---\left(13\right)$

Σ′_i位于局部一般位置，它的所有实根有不同的坐标，与Σ_i的实根一一对应，并且重数相同。

(iv)将所有实根投影到一维空间中，得到一个单变元方程：

计算结式列

的实根与Σ_i的实根一一对应，并且只差一个线性变换。

(v)求计算单变元多项式方程的实根

(vi)通过线性变换

${\mathrm{\ξ}}_1=\left(\underset{j=1}{\overset{i-1}{\mathrm{\Π}}}\frac{R_{j+1}}{r_j}\right)({\mathrm{\η}}_i-{\mathrm{\η}}_{i-1}),---\left(15\right)$

η_i∈Zero(T_i ⁱ)， ${\mathrm{\&eta;}}_{i-1}\&Element;\mathrm{Zero}\left(T_{i-1}^{i-1}\right),$ $|{\mathrm{\&eta;}}_i-{\mathrm{\&eta;}}_{i-1}|<\left(\underset{j=1}{\overset{i-2}{\mathrm{\&Pi;}}}\frac{r_j}{R_{j+1}}\right)r_{i-1}$

将一元实根η_i提升回原方程组的实根(ξ₁，ξ₂，…，ξ₁)。

(c)令i＝n，执行操作(b)(i)，(ii)，计算得到r_n-1和R_n；

(d)对Σ_n做线性变换

得到新的多项式系统

Σ′_n位于局部一般位置，它的所有实根有不同的坐标，与Σ_n的实根一一对应，并且重数相同。

II.将所有实根投影到一维空间中，得到一个单变元方程；

计算结式列

的实根与Σ_n的实根一一对应，只差一个线性变换，并且重数相同。

III.求计算单变元多项式方程的所有实根及其重数。

IV.通过线性变换

${\mathrm{\&xi;}}_n=\left(\underset{j=1}{\overset{n-1}{\mathrm{\&Pi;}}}\frac{R_{j+1}}{r_j}\right)({\mathrm{\&eta;}}_n-{\mathrm{\&eta;}}_{n-1}),---\left(19\right)$

${\mathrm{\&eta;}}_n\&Element;\mathrm{Zero}\left(T_n^n\right),$ ${\mathrm{\&eta;}}_{n-1}\&Element;\mathrm{Zero}\left(T_{n-1}^{n-1}\right),$ $|{\mathrm{\&eta;}}_n-{\mathrm{\&eta;}}_{n-1}|<\left(\underset{j=1}{\overset{n-2}{\mathrm{\&Pi;}}}\frac{r_j}{R_{j+1}}\right)r_{n-1}$

将一元实根η_n提升回原方程组的实根(ξ₁，ξ₂，…，ξ_n)。

这样就得到了Σ_n＝{f₁(x₁)，f₂(x₁，x₂)，…，f_n(x₁，x₂，…，x_n)}的所有实根：

Zero(Σ_n)＝{(ξ₁，ξ₂，…，ξ_n)∈Rⁿ|f₁(ξ₁)＝…＝f_n(ξ₁，…，ξ_n)＝0}  (20)

从而完成计算过程。

2.生成曲面连接关系图。

曲面的连接关系图包括三部分：特征点(点集合)，特征点之间的连接关系(边集合)，特征点构成的三角片与边的连接关系(面集合)。

点集合在第一步已经得到。

边集合包括独立奇点与奇异线上点的连接关系，奇异线上点与曲面上一般点的连接关系。对一个独立奇点，与他连接的奇异线的条数，可以通过计算邻接他的奇异线上的点的个数来确定；对一个奇异线上的点，与他连接的曲面上的点的个数，可以通过计算邻接他的曲面上一般点的个数来确定。如图6所示。

面集合是边集合中边与曲面的连接情况。可以通过计算连接到同一边两个端点的点的个数来确定。如图7所示。

3.绘制曲面。

步骤2中得到的面集合已经是曲面的一个三角网格化，我们可以直接用他来绘制网格曲面，或者在此基础上通过其他光顺化方法绘制更光滑的曲面，比如常用的MATLAB，OPENGL，MATH，MAYA等绘制软件平台均可以进行该三维绘制。效果如图8，图9，图10所示，相对于现有的三维绘制技术，本发明的技术方案绘制出了更加准确的空间立体结构，解决了奇点处及其邻近区域网格结构绘制不准确的技术问题。

相应地，本发明还设计了基于上述通过计算零维三角多项式系统所有实根及其重数进行曲面绘制的方法的系统，如图11所示，本系统包括：

隐式曲面方程生成装置；

曲面特征点生成装置，其通过计算零维三角多项式系统所有实根及其重数得到曲面的特征点及其重数；

曲面连接图生成装置，其通过确定特征点之间的连接关系得到连接关系图；

曲面绘制装置，其在连接关系图的基础上直接绘制网格曲面或者通过其他光顺化方法绘制更光滑的曲面；

显示装置，显示绘制得到的曲面图像。

Claims (5)

1.一种通过计算零维三角多项式系统所有实根及其重数进行曲面绘制的方法，包括以下步骤：

1)确定曲面的特征点；特征点包括曲面的独立奇点，奇异线上的点，与奇点邻近的非奇点；其包括如下步骤：

(1)确定特征点满足的零维三角多项式系统；

(2)计算零维三角多项式系统的所有实根及其重数；其包括以下步骤：

I.对零维三角多项式系统Σ_n＝{f₁(x₁)，f₂(x₁，x₂)，…，f_n(x₁，x₂，…，x_n)}进行预处理，通过线性变换将其变换到局部一般位置；

II.将所有实根投影到一维空间中，通过计算结式来进行投影，得到一个单变元方程，该投影保持根的重数不变；

III.计算单变元多项式方程的所有实根及其重数；

IV.通过线性变换将一元实根提升回原方程组的实根；

2)生成曲面连接关系图

曲面的连接关系图包括三部分：点集合，即步骤1)确定曲面的特征点集合；边集合，即特征点之间的连接关系；面集合，即特征点构成的三角片与边的连接关系；

3)绘制曲面

步骤2)中得到的面集合已经是曲面的一个三角网格化，以此直接绘制网格曲面，或者在此基础上通过其他光顺化方法绘制更光滑的曲面。

2.根据权利要求1所述的方法，其中，对零维三角多项式系统预处理包括：

1)计算实根的隔离界r_i，i＝1，2，…n-1；

2)估计根的上界R_i，i＝2，3，…n；

3)确定线性变换，把多项式系统变换到局部一般位置。

3.根据权利要求2所述的方法，其中，估计根的上界包括：

设实数列(ξ₁，ξ₂，…，ξ_i-1)由一个包含它且端点为有理数的区间列[a₁，b₁]×[a₂，b₂]×…×[a_i-1，b_i-1]表示，将多项式拆分：

f_i(x₁，x₂，…，x_i)＝f_i ⁺-f_i ^-，

其中f_i ⁺，f_i ^-∈Z⁺[x₁，x₂，…，x_i]，是满足上式的单项式项数最少的正系数多项式；

定义上界多项式f_i ^u和下界多项式f_i ^d(x_i)：

f_i ^u(x_i)＝f_i ⁺(b₁，…，b_i-1，x_i)-f_i ^-(a₁，…，a_i-1，x_i)，

f_i ^d(x_i)＝f_i ⁺(a₁，…，a_i-1，x_i)-f_i ^-(b₁，…，b_i-1，x_i)，

设f_i ^u(x_i)形为

${f_i}^u\left(x_i\right)=c_d{x}_i^d+c_{d-1}{x}_i^{d-1}+...+c_0,$

计算f_i ^u(x_i)的根的上界

$\mathrm{Ru}=1+\underset{0\&le;k\&le;d-1}{\max }\left|\frac{c_k}{c_d}\right|,$

运用同样的方法计算根的下界Rd；则f_i(ξ₁，ξ₂，…，ξ_i-1，x_i)的根的上界为

R＝max{Ru，Rd}

最后计算f_i(x₁，x₂，…，x_i)关于x_i的根的上界

R_i＝max{R：R是f_i(ξ₁，ξ₂，…，ξ_i-1，x_i)根的上界，(ξ₁，ξ₂，…，ξ_i-1)∈Zero(Σ_i-1)}。

4.根据权利要求2所述的方法，其中确定线性变换，把多项式系统变换到局部一般位置包括：

对Σ_n做线性变换

得到新的多项式系统

Σ′_n位于局部一般位置，它的所有实根有不同的坐标，与Σ_n的实根一一对应，并且重数相同。

5.根据权利要求1所述的方法，其中，通过线性变换将一元实根提升回原方程组的实根包括：

通过线性变换

${\mathrm{\&xi;}}_1={\mathrm{\&eta;}}_1,{\mathrm{\&xi;}}_i=\left(\underset{j=1}{\overset{i-1}{\mathrm{\&Pi;}}}\frac{R_{j+1}}{r_j}\right)({\mathrm{\&eta;}}_i-{\mathrm{\&eta;}}_{i-1}),$

η_i∈Zero(T_i ⁱ)， ${\mathrm{\&eta;}}_{i-1}\&Element;\mathrm{Zero}\left(T_{i-1}^{i-1}\right),$ $|{\mathrm{\&eta;}}_i-{\mathrm{\&eta;}}_{i-1}|<\left(\underset{j=1}{\overset{i-2}{\mathrm{\&Pi;}}}\frac{r_j}{R_{j+1}}\right)r_{i-1},$ i＝2，3，…，n

将一元实根η₁，η₂，…，η_n提升回原方程组的实根(ξ₁，ξ₂，…，ξ_n)；这样就得到了Σ_n＝{f₁(x₁)，f₂(x₁，x₂)，…，f_n(x₁，x₂，…，x_n)}的所有实根：

Zero(Σ_n)＝{(ξ₁，ξ₂，…，ξ_n)∈Rⁿ|f₁(ξ₁)＝…＝f_n(ξ₁，…，ξ_n)＝0}。