- Email: [email protected]
Convergence analysis for B-spline geometric interpolation
Computers & Graphics 36 (2012) 884–891
Contents lists available at SciVerse ScienceDirect
Computers & Graphics journal homepage: www.elsevier.com/locate/cag
Computer Graphics in China
Convergence analysis for B-spline geometric interpolation Yunhui Xiong a, Guiqing Li b,n, Aihua Mao b a b
School of Science, South China University of Technology, Guangzhou, China School of Computer Science & Engineering, South China University of Technology, Guangzhou, China
a r t i c l e i n f o
a b s t r a c t
Article history: Received 10 December 2011 Received in revised form 10 July 2012 Accepted 10 July 2012 Available online 19 August 2012
In this paper, we propose a sufficient condition for the convergence of a geometric algorithm for interpolating a given polygon using non-uniform cubic B-splines. Geometric interpolation uses the given polygon as the initial shape of the control polygon of the B-spline and reduces the approximate error by iteratively updating the control points with the deviations from the corresponding interpolated vertices to their nearest footpoints on the current B-spline curve. The convergence condition is derived by employing a spectral radius estimation technique. The primary goal is to find for each control point a parametric interval within which the nearest footpoint should be confined such that the spectral radius of the error iteration matrix is smaller than 1. A convergent condition for the geometric interpolation of uniform B-splines can be derived as a special case of the new scheme. & 2012 Elsevier Ltd. All rights reserved.
Keywords: Geometric interpolation Convergence analysis Non-uniform B-splines Spectral analysis
1. Introduction The construction of a parametric curve passing through a set of scattered points is a fundamental task in computer aided geometric design due to its extensive applications in a variety of fields, such as economics, medicine, computer animation, manufacturing, and signal processing (see [1–6] for example). Though one may impose different requirements on generated curves depending on the application [1,7], the interpolatory curves should typically look pleasing and be as fair as possible. Interpolation problem based on B-spline curves can be formally described as follows: given a closed polygon P ¼ ðp1 ,p2 ,. . .,pn Þ and B-spline representation CðtÞ ¼
n X
Bi ðtÞp0i ,
i¼1
where Bi(t) and p0i are the non-uniform B-spline basis functions and the unknown control points, respectively, we need to determine p0i and t 1 o t 2 o ot n such that Cðt i Þ ¼ pi ,
i ¼ 1,2, . . ., n:
ð1Þ
Eq. (1) is actually a nonlinear under-constrained system with many solutions. In the literature, many approaches to obtaining an optimal solution are described. A straightforward method for obtaining an optimal solution is to specify {ti} explicitly and reduce Eq. (1) to a linear system. This n
Corresponding author. E-mail addresses: [email protected] (Y. Xiong), [email protected] (G. Li), [email protected] (A. Mao). 0097-8493/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.cag.2012.07.002
method usually leads to results with overshooting and undulating artifacts. Many parameterization methods have been proposed to remedy these problems [8–10]. Among these, arc-length parameterization and its variants are the most popular in practice. However, no solution is the best for all cases [11–13]. Some authors consider addressing the problem from the viewpoint of variation, but this method usually results in complicated optimizations [14]. In [11], Hoschek proposed an intrinsic parameterization to attack the curve fitting problem. Using this method, parameter values are computed implicitly by finding the nearest point from the data points to the fitting curve and are then employed to fit a new curve. The point nearest to a data point is called its footpoint. This approach iteratively refines the solution by repeating the process. Saux and Daniel [12] improved the efficiency of the approach, and Wang et al. [15] extended the approach with several new fitting errors. Unfortunately, these approaches degenerate if applied to interpolating problems instead of fitting tasks. Recently, Maekawa et al. [16] exploited a geometric interpolation algorithm to interpolate the vertices of a given mesh using subdivision surfaces. Their algorithm iteratively modifies each control vertex of the mesh using a deviation from the footpoint to its corresponding data point. In [16], two types of footpoints are discussed: the nearest point on the subdivision surface to the data point and the limit position of the corresponding control point on the surface to the data point (parametric distance deviation). As neither a large linear system nor the minimization of smooth energies is required, this approach is quite efficient and can generally produce high-quality results using the first type of correction due to its essence of intrinsic parameterization. In [17], Chen et al. extended this algorithm to the Catumll– Clark subdivision surfaces. Gofuku et al. [18] further studied the interpolation of tangents and normals using the geometric algorithm.
Y. Xiong et al. / Computers & Graphics 36 (2012) 884–891
However, the convergence analysis of the geometric interpolation algorithm remains open except in the special case of quadratic B-spline curves with the second type of deviation [15]. Cheng et al. showed that the algorithm is also convergent for Loop surfaces with the parametric distance deviation (the second type of deviation) [19]. Lin et al. introduced a class of iterative nonuniform B-spline curves [20] and presented a sufficient condition for the convergence of the geometric interpolation algorithm [21]. Carnicer et al. [22] further improved Lin’s proof. However, the convergence condition is too high level to derive a practical convergent scheme. In [23], Xiong et al. explored a sufficient condition that can guarantee the convergence of geometric interpolation based on a uniform B-spline. The limitation is that this condition does not work when using non-uniform B-splines due to the complexity of non-uniform B-spline basis functions. As an extended version of reference [23], this paper considers the case of non-uniform cubic B-spline curve interpolation. Our approach is also applicable to B-splines of degree 2–6 in which the maximum of the basis functions of the B-splines is greater than 0.5. The main contribution of this paper is the exploration of a sufficient condition that guarantees the convergence of the geometric interpolation algorithm. The convergent conditions for uniform B-spline geometric interpolation described in [23] is then derived as a special case of our new formulation. The remainder of this paper is organized as follows. Section 2 explains the mechanism of geometric interpolation using a simple case where only one control point of a given polygon is interpolated. Section 3 describes the proposed scheme in detail and then analyzes its convergence. In Section 4, several examples are given to examine the quality of interpolating curves using our new scheme. Finally, we conclude our paper in Section 5.
885
m Step 3 Vertex pm k in P is modified as follows
þ1 m pm ¼ pm k þ ek , k m m where the correction term em k ¼ pk C ðt k Þ . Step 4 Set m’m þ1, and repeat steps 2 and 3 until the correction term is less than a given threshold or the iterative step greater than a specified value. 0 Step 5 p0k ¼ pm k , and return P .
In Algorithm 1, the footpoint of pk on the current curve Cm(t) is required. Selecting nearest points as footpoints may yield a polynomial equation. If more than one nearest point is found, we choose the point whose parameter has the nearest value to uk. Then, we have the following results. Theorem 1. The convergence of Algorithm 1 can be guaranteed for a footpoint whose parametric value is constrained within ui1 r t m k r ui þ 1 , i.e., lim Hd ð pk ,C m ðtÞÞ ¼ 0, m-1
where Hd(A, B) is the Hausdorff distance from A to B. Proof. According to Algorithm 1, we have þ1 C m þ 1 ðtÞ ¼ ðpm pk ÞBk ðtÞ þ k
n X
pi Bi ðtÞ
i¼1 m ¼ ½pk C m ðt m k ÞBk ðtÞ þC ðtÞ, þ1 where {Bi(t)} are cubic B-spline basis functions. As C m þ 1 ðt m Þ is k mþ1 the nearest point to pk on the segment ui 1 rt rui þ 1 of C (t), þ1 :pk C m þ 1 ðt m Þ: r:pk C m þ 1 ðt m k Þ: k
2. Geometric interpolation to one control point with B-spline curves Given a polygon P ¼ ðp1 ,p2 ,. . .,pn Þ, we consider such a problem that constructs a non-uniform cubic B-spline C(t) defined by control polygon P 0 ¼ ðp1 , p2 ,. . .,p0k ,. . .,pn Þ and knot vector U ¼ ðu1 ,u0 ,u1 ,u2 ,. . .,uk ,. . .,un ,un þ 1 ,un þ 2 Þ, u1 ¼ un1 un ,u0 ¼ 0,un þ 1 ¼ un þ u1 , un þ 2 ¼ un þu2 , such that C(t) passes through the k-th vertex pk of polygon P. We modify the mechanism described in [23] slightly to address this single point interpolation problem using non-uniform B-splines. The new scheme also iteratively updates the position of the k-th control point and sustains other vertices unchanged. (See the pseudocode in Algorithm 1.)
m m m ¼ :pk ½pk C m ðt m k ÞBk ðt k ÞC ðt k Þ: m m ¼ :1Bk ðt m k Þ::pk C ðt k Þ:
¼
m Y
91Bk ðt jk Þ9:pk C 0 ðt 0k Þ::
j¼0
For
t jk A [ui 1,ui þ 1],
we can always guarantee 91Bk ðt jk Þ9 o1.
&
Remark 1. The convergent rate of Algorithm 1 is 5/6 for uniform cubic B-splines because 1/6rBk(t) r4/6 in this case. Remark 2. If the footpoint is sought on the whole Cm(t), Algorithm 1 is not guaranteed to converge [23]. 3. Geometric interpolation to closed polygons with B-spline curves
Algorithm 1. Constructing a B-spline curve interpolating a specified control point of the initial polygon Input
In this section, we consider solving Eq. (1) based on nonuniform cubic B-spline curves. Again, let
P ¼ ðp1 ,p2 ,. . .,pn Þ,
u1 ¼ un1 un ,u0 ¼ 0,un þ 1 ¼ un þu1 ,un þ 2 ¼ un þ u2 ,
U ¼ ðu1 ,u0 ,u1 ,u2 ,. . .,uk ,. . .,un ,un þ 1 ,un þ 2 Þ,
be the knot vector. Our geometric interpolation scheme works as described in Algorithm 2.
u1 ¼ un1 un ,u0 ¼ 0,un þ 1 ¼ un þ u1 ,un þ 2 ¼ un þ u2 : Output P 0 ¼ ðp1 ,p2 ,. . .,p0k ,. . .,pn Þ:
U ¼ ðu1 ,u0 ,u1 ,u2 ,. . .,uk ,. . .,un ,un þ 1 ,un þ 2 Þ,
Algorithm 2. The construction of the non-uniform B-spline curve interpolating all vertices of a given polygon. Input P ¼ ðp1 ,p2 ,. . .,pn Þ,
Step 1 Start from the given polygon: P 0 ¼P. Step 2 Assume we have gained P m ¼ ðp1 , ,pm k , ,pn Þ . Find the footpoint of pk on the curve Cm(t), which stands for the B-spline with control polygon Pm and knot vector U.
U ¼ ðu1 ,u0 ,u1 ,u2 ,. . .,uk ,. . .,un ,un þ 1 ,un þ 2 Þ, u1 ¼ un1 un ,u0 ¼ 0,un þ 1 ¼ un þu1 ,un þ 2 ¼ un þ u2 :
886
Y. Xiong et al. / Computers & Graphics 36 (2012) 884–891
Assume that Bi(u)(ui 1 otoui þ 1) achieves maximum at parameter value ti,max, namely, Bi ðt i,max Þ ¼ max Bi ðuÞ : t A ðui1 ,ui þ 1 Þ :
Output 0
P ¼
ðp01 ,p02 ,
,p0n Þ
0
Step 1 Start from the given polygon as the initial guess of P , namely, let P0 ¼P. m m Step 2 Assume we have obtained P m ¼ ðpm 1 ,p2 ,. . .,pn Þ: The control points of Pm þ 1 are calculated as follows m m þ1 pm ¼ pm i þ ðpi C ðt i ÞÞ i m
i ¼ 1,2,. . .,n
ð2Þ m
ðt m i Þ
where C is the footpoint of pi on C (t) within a given interval. (See Lemma 2.) Step 3 Set m’m þ1, and repeat Step 2 until specified thresholds are satisfied. 0 0 Step 4 P ¼Pm, and return P . From (2), we have C m þ 1 ðtÞ ¼
n X
þ1 pm Bj ðtÞ j
j¼1
¼
n X
m m ½pm j þðpj C ðt j ÞÞBj ðtÞ
j¼1
¼ C m ðtÞ þ
n X
½pj C m ðt m j ÞBj ðtÞ:
j¼1
It follows that n X
m ½pj C m ðt m j ÞBj ðt i Þ
m m ¼ ½B1 ðt m i Þ, ,Bi1 ðt i Þ,1Bi ðt i Þ, m m Bi þ 1 ðt m i Þ,. . .,Bn ðt i ÞD ,
i ¼ 1,2,. . .,n,
ð3Þ
where Dm is a column vector consisting of the correction terms m m m m T Dm ¼ ½p1 C m ðt m 1 Þ,p2 C ðt 2 Þ,. . .,pn C ðt n Þ
Eq. (3) has the following matrix form ¼ ðIBm ÞDm ,
D
Proof. We only need to show that there exists a real number e0A(0, 1) such that r(I Bm) r1 e0. As rðIBm Þ r:IBm : holds for arbitrary norm :: , our task is to find a concrete norm such that :IBm :r 1e0 holds under the given condition. Consider the row-sum norm of I–Bm 9 8 = < X :IBm :1 ¼ max 1Bi ðt m Bj ðt m i Þþ i Þ , ; i ¼ 1,...,n: 1 r j r n,j a i
where mþ1
Proof. The continuity of the basis functions implies that Bi ðtÞ Z 1=2 þc holds for some open interval containing ti,max . We can compute di,1 and di,2 by solving the polynomial equation Bi ðtÞ1=2c ¼ 0. & Bi ðt i,max Þ 41=2 is not always true. Fortunately, this is only the case for certain extreme situations in which the knot intervals are very non-uniform. To demonstrate this, we randomly generate 10,000 five-dimension knot vectors and then check the corresponding basis functions one by one for their maximum. Fig. 1 shows that only 25 examples include basis functions whose maxima (tagged with ‘ ’) are not greater than 0.5. For each unsuccessful case, we calculate the ratio max uj uj1 , i1r j r iþ 2 , min uj uj1 , i1 r j ri þ 2
Theorem 2. Algorithm 2 is convergent if the parametric value t m i of the footpoint of pi is confined within the interval [di,1,di,2] derived in Lemma 2 for i¼1, 2, y, n.
j¼1
mþ1
Lemma 2. For non-uniform cubic B-splines, there exist two real numbers di,1 and di,2 that satisfy ui1 r di,1 r t i,max r di,2 r ui þ 1 such that Bi ðtÞ Z 1=2 þc for an arbitrary tA[di,1,di,2], i¼1, 2, y, n.
and find that the ratio is greater than 16 for all these cases.
m m pi C m þ 1 ðt m i Þ ¼ pi C ðt i Þ
D
Let c be a user specified number satisfying Bi ðui,max Þ Z1=2 þ c if Bi ðui,max Þ 4 1=2 . We then have
¼ ½p1 C
mþ1
ðt m 1 Þ,p2 C
mþ1
ðt m 2 Þ,. . .,pn C
mþ1
T ðt m n Þ
and m
which implies that we should guarantee X 1Bi ðt m Bj ðt m i ¼ 1,2,. . .,n: i Þþ i Þ r 1e0 , 1 r j r n,j a i
m
Bm ¼ ðbij Þnn with bij ¼ Bj ðt m i Þ, Noting that
Noting the affine property of B-spline basis functions, the above inequality is equivalent to
þ1 :pj C m þ 1 ðt m Þ:r :pj C m þ 1 ðt m j Þ:, j
2ð1Bi ðt m i ÞÞ r1e0 : Namely,
:Dm þ 1 :r :D
mþ1
: r:IBm ::Dm : r r
m Y
:IBk ::D1 :,
k¼1
and the composite matrix following lemma.
m Q
:IBk : This leads to the
k¼1
Lemma 1. If there exist a positive constant e0 and an integer m0 such that the spectral radius r(I Bm) of I Bm satisfies r(I Bm)r1 e0 for all m4m0, then Algorithm 2 converges. Proof. It is obvious that lim :Dm þ 1 : ¼ 0 if the condition in m-1 Lemma 1 is satisfied. & Lin et al. derived a special case of Lemma 1 in [21]. This special case is not applicable yet because it does not address how to update Algorithm 2 such that r(I Bm)r1 e0 holds.
Bi ðt m i ÞZ
1 þ e0 : 2
ð4Þ
Set e0 ¼2c. According to Lemma 2, if tm i A ½di,1 , di,2 , then ð1 þ e0 Þ=2 ¼ 1=2 þc rBi ðt m & i Þ .
1 0.8 0.6 0.4 1 625 1249 1873 2497 3121 3745 4369 4993 5617 6241 6865 7489 8113 8737 9361 9985
we obtain
Fig. 1. The maximum distribution of the 10,000 basis functions.
Y. Xiong et al. / Computers & Graphics 36 (2012) 884–891
Remark 3. If Bi ðui,max Þ r1=2 appears occasionally we just set tm i ¼ ui to make the algorithm more robust. Corollary 1. Algorithm 2 is convergent for uniform cubic B-splines if the parameters t m i of the footpoints are confined to interval ½di,1 , di,2 , i¼1, 2, y, n, where di,1 ¼ i0:466178 and di,2 ¼iþ0.466178. Proof. See reference [23] for the proof. & In each iterative step, Algorithm 2 updates the control points simultaneously. Using recently renewed control points to evaluate the later control points results in Algorithm 3 is sometimes more efficient as demonstrated in Section 4. Algorithm 3. The construction of B-spline curve interpolating all vertices of an initial polygon using the instant updating mode. Step 1 Start from the given mesh as the initial control mesh, i.e., P 0 ¼P. þ1 Step 2 At the (mþ 1)-th iteration, P m þ 1 ¼ ðp1m þ 1 ,pm , ,pnm þ 1 Þ 2 þ1 mþ1 is renewed as follows. Assuming that p1m þ 1 ,pm ,. . .,p have 2 i1 been calculated, we first introduce an intermediate B-spline C m,i ðtÞ ¼
i1 X j¼1
þ1 pm Bj ðtÞ þ j
n X
pm j Bj ðtÞ,
j¼i
m,i and then find the footpoint (denoted by C m,i ðt m (t). i Þ) of pi on C mþ1 pi is finally modified by
m,i m pim þ 1 ¼ pm ðt i ÞÞ, i þ ðpi C
i ¼ 1,2,. . .,n,
ð5Þ
We will show the convergence of Algorithm 3 in Section 4 using numerical examples. 4. Experimental results and discussion We perform our experiments on a PC with a 2.40 GHz Intel Xeon E5620 CPU, 4 GB of RAM, and the Windows 7 OS. Five variants of our algorithms (4 for Algorithms 2 and 1 for Algorithm 3) summarized in Table 1 are implemented using MATLAB R2008a for non-uniform cubic B-splines. The distinction among these variants is that each seeks footpoints within differently sized parametric ranges. The footpoints (i.e., the nearest points) are obtained by finding the zero points of the derivative of the squared distance from the data point to the current B-spline using the MATLAB function roots(). We also conduct three types of experimental comparison in this section. The first one is used to investigate the efficiency of the variants, the second is used to show the spectra of the composite matrices numerically, and the third is used to compare the quality of the generated curves. All of the nonuniform knot vectors in the experiments are created using the chord length parameterization. To observe the efficiency of the variants, we generate 11 groups of polygon data. Each group includes 200 control polygons and 2200 polygons in total. The polygons in the same group have the same vertex number (5, 10, 15, 20, 30, 40, 50, 60, 70, 80, and 90 vertices). The positions of all of the polygonal vertices are randomly generated in the square [0,5] [0,5]. The convergent Table 1 Five variants of the interpolating algorithm. Range of Variant Variant Variant Variant Variant
0 1 2 3 4
tm i
(ui 2,ui þ 2) (ui 1,ui þ 1) (di,1,di,2) ui (di,1,di,2)
Updating mode Algorithm Algorithm Algorithm Algorithm Algorithm
2 2 2 2 3
887
criterion is that the maximal distance from data points to their footpoints is less than 10 6 or the iteration steps exceed 1000.
4.1. Convergence observation As the search range of variant 0 is far beyond our sufficient condition, many cases are expected to diverge. The experiment demonstrates that variant 0 fails to converge in the numerical sense for almost all examples when the control polygon has more than 40 vertices (see Fig. 2). Though variant 1 is numerically convergent for all of the above data, we find some abortive cases when increasing the number of control vertices to 300. In these cases, the distance error exhibits fluctuation around a small value (about 10 4).
4.2. Efficiency analysis Because variant 1 seeks the footpoints in a larger range than variants 2 and 4, its convergent rate is slightly higher than the convergent rates of variants 2 and 4. Variant 3 interpolates data points with fixed parametric values and therefore has the lowest convergent speed among the four schemes. This accords with the experimental results depicted in Fig. 3 in which the average iteration step for the polygons with the same vertex number is shown for the four schemes. In addition, the iteration step increases slowly for the variants other than variant 3. Nevertheless, the computational time is not completely consistent with the above order. Because variant 3 does not need to seek the footpoint, which is the most time-consuming operation in the algorithm, it is more efficient than variants 1, 2 and 4 for converging to the specified threshold, as shown in Fig. 4. It also shows that the time complexity of variants 1, 2 and 4 is very close. Variant 1 is slightly faster than the other two because of the higher convergent rate.
120% 100% 80% 60% 40% 20% 0%
5 10 15 20 30 40 50 60 70 80 90
Fig. 2. Percentage of those examples in which variant 0 converges against all examples (vertical axis: percent; horizontal axis: vertex number).
350 300 250 200 150 100 50 0
variant1 variant2 variant3 variant4 5 10 15 20 30 40 50 60 70 80 90
Fig. 3. The iteration step comparison among the four variants. vertical axis: the average times of iterator; horizontal axis: the vertex number of polygons.
888
Y. Xiong et al. / Computers & Graphics 36 (2012) 884–891
4.3. Spectral analysis To show the behavior of the composite matrices numerically, we have computed the spectral radius of the composite matrices for each round of iteration for every randomly generated polygon. As expected, the spectral radius for variant 0 increases for most of the experiments (see Fig. 5(a)). All experiments for variants 1, 2 and 3 are numerically convergent (Fig. 5(b–d)), and the
6 5 4
variant1
3
variant2
2
variant3
1
variant4
0
5 10 15 20 30 40 50 60 70 80 90
Fig. 4. The average time cost comparison among the four variants. vertical axis: the average time cost (second); horizontal axis: vertex number.
14
x 108
convergent rate of variant 1 appears to be slightly faster than that of variant 2. Again, variant 3 demands the maximal steps of iterations. 4.4. Quality comparison Though the range for seeking footpoints is confined to a small interval, variant 2 can still produce high quality curves. Figs. 6 and 7 illustrate that the curvature maps of curves generated by variant 2 are very close to those of variant 1. Variant 3, which is actually a common interpolating technique, generally results in curves with the worst quality. The other two experiments demonstrate an application of the algorithm in Chinese character generation. A Chinese character ‘‘Dragon’’ (Fig. 8) in Huawen Xingkai font and a Chinese character ‘‘Good’’ in Huawen Lishu font (Fig. 9) are used as examples here. Their contour polygons are extracted using the MATLAB function ginput(). We then apply the proposed scheme to construct the spline curves interpolating the character contours. Different simple polygons in the same character are handled separately in our implementation. Again, the results of variants 1 and 2 look very nice (as shown in Figs. 8(c and d) and 9(c and d)), but the results of variant 3 exhibit overshooting artifacts in some corners
1 0.9
12
0.8
10
0.7 0.6
8
0.5
6
0.4 0.3
4
0.2
2 0
0.1
0
5
10
15
20
25
30
35
40
45
50
0
1
1
0.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0
10
20
30
40
50
60
70
80
90
100
0
0
10
20
30
40
50
60
70
80
90
100
0
10
20
30
40
50
60
70
80
90
100
Fig. 5. The behavior of the spectral radius for the composite matrices of the 2200 non-uniform cubic B-splines. (The maximal iteration step is 1000.) (a, b, c, d) respectively correspond to variants 0, 1, 2, and 3. In the images, the horizontal and vertical axes stand for the spectral radius and the iteration step, respectively, and each curve plots the radius of a B-spline as a function of the iteration step.
Y. Xiong et al. / Computers & Graphics 36 (2012) 884–891
889
Fig. 6. Curvature maps. Top row: (a) variant 1, (b) variant 2, (c) variant 3, (d) the polygon to be interpolated. Bottom row: zoomed in images corresponding to the part of images marked by the dotted rectangles in the top row.
Fig. 7. Curvature maps. Top row: (a) variant 1, (b) variant 2, (c) variant 3, (d) the polygon to be interpolated. Bottom row: zoomed in images corresponding to the part of images marked by the dotted rectangles in the top row.
Fig. 8. Interpolation generation of the Chinese character ‘Dragon’: (a) original character; (b) polygon approximation of the character contour; (c, d, e) (middle row): the interpolation results of variants 1, 2, and 3 respectively; bottom row: zoomed in images corresponding to the part of images marked by the dotted rectangles in the middle row.
890
Y. Xiong et al. / Computers & Graphics 36 (2012) 884–891
Fig. 9. Interpolation generation of the Chinese character ‘Good’: (a) original character; (b) polygon approximation of the character contour; (c, d, e) (middle row): the interpolation results of variants 1, 2, and 3 respectively; bottom row: zoomed in images corresponding to the part of images marked by the dotted rectangles in the middle row.
of the characters (as shown in the zoom-in regions shown in Figs. 8e and 9e). We can conclude from aforementioned experiments that variant 2 can produce high-quality interpolating curves and can guarantee convergence except in some extreme cases.
5. Conclusions Geometric interpolation, which is an iterative method to construct control polygons or meshes, comes from Maekawa et al. [16], in which the convergence of the iterative method was not addressed. In this paper, a sufficient condition of convergence, which constrains the footpoint search to a local range, is presented for non-uniform B-splines geometric interpolation. Experiments demonstrate that the local method can still produce high-quality interpolating curves while enhancing the convergence speed greatly compared to the global method using non-uniform cubic B-splines. As for future work, it is interesting to give a theoretical condition in terms of interval lengths to guarantee the maximum of non-uniform B-spline bases greater than 0.5. For convenience, our discussion only involves the interpolation of closed polygons.
Bi ðtÞ
¼
The convergence analysis of the control polyline interpolation must be supplied with clamped conditions to apply the approach. Our goal is to develop a similar convergent condition for the geometric interpolation of parametric surfaces and subdivision surfaces. Another interesting problem is to investigate the mechanism that makes this method capable of producing smoother results than interpolation schemes specifying parameter values.
Acknowledgments We are delighted to thank the reviewers for their valuable comments. The work is partially supported by NSFC (60973084, 61003173), NSF of Guangdong (915106410-1000106), and Fundamental Research Funds for the Central Universities (2012zz0063, 2012zm0062).
Appendix A. Piecewise representation of base function Bi(t)
8 ðtui2 Þ3 > > ðui þ 1 ui2 Þðui ui2 Þðui1 ui2 Þ > > h i > > ðui þ 1 tÞðtui1 Þ ðui þ 2 tÞ ðtui2 Þ ðtui2 Þðui tÞ ðtui1 Þ2 > > > ðui þ 1 u Þ ðui ui2 Þðui ui1 Þ þ ðui þ 1 ui1 Þðui ui1 Þ þ ðui þ 2 ui1 Þ ðui þ 1 ui1 Þðui ui1 Þ > i2 < h i ðui þ 1 tÞ2 ðui þ 2 tÞ ðtui1 Þðui þ 1 tÞ ðui þ 2 tÞðtui Þ ðtui2 Þ > ðui þ 1 ui2 Þ ðui þ 1 ui Þðui þ 1 ui1 Þ þ ðui þ 2 ui1 Þ ðui þ 1 ui1 Þðui þ 1 ui Þ þ ðui þ 2 ui Þðui þ 1 ui Þ > > > > ðui þ 2 tÞ3 > > > > > ðui þ 2 ui þ 1 Þðui þ 2 ui Þðui þ 2 ui1 Þ : 0
t A ðui2 ,ui1 Þ t A ðui1 ,ui Þ t A ðui ,ui þ 1 Þ t A ðui þ 1 ,ui þ 2 Þ others
Y. Xiong et al. / Computers & Graphics 36 (2012) 884–891
References [1] Hagan P, West G. Interpolation methods for curve construction. Appl Math Finance 2006;13(2):89–129. [2] Jacob M, Blu T, Unser M. Efficient energies and algorithms for parametric snakes. IEEE Trans Image Process 2004;13(9):1231–44. [3] Ramamoorthi R, Barr AH. Fast construction of accurate quaternion splines. In: computer graphics Proceedings, annual conference series, ACM SIGGRAPH. Los Angeles; 1997. p. 287–92. [4] Nam S-H, Yang M-Y. A study on a generalized parametric interpolatior with real-time jerk-limited acceleration. Comput Aided Des 2004;36(1):27–36. [5] Wang F-C, Yang D. Nearly arc-length parameterized quintic-spline interpolation for precision machining. Comput Aided Des 1993;25(3):281–8. [6] Levin A. Interpolating nets of curves by smooth subdivision surfaces. In: computer graphics Proceedings, annual conference series, ACM SIGGRAPH. Los Angeles; 1999. p. 57–64. [7] Piegl LA, Rajab K, Smarodzinana V. Curve interpolation with directional constraints for engineering design. Eng Comput 2008;24:79–85. [8] Floater MS, Surazhsky T. Parameterization for curve interpolation. Stud Comput Math 2006;12:39–54. [9] Papaioannou SG, Patrikoussakis MM. Curve interpolation based on the canonical arc length parametrization. Comput Aided Des 2011;43:21–30. [10] Kouibia A, Pasadas M, Rodrı´guez ML. Optimization of parameters for curve interpolation by cubic splines. J Comput Appl Math 2011http://dx.doi.org/ 10.1016/j.cam.2011.03.013. [11] Hoschek J. Intrinsic parametrization for approximation. Comput Aided Geom Des 1988;5(2):27–31. [12] Saux E, Daniel M. An improved Hoschek intrinsic parametrization. Comput Aided Geom Des 2003;20(8–9):513–21.
891
[13] Pottmann H, Leopoldseder S. A concept for parametric surface fitting which avoids the parametrization problem. Comput Aided Geom Des 2003;20: 243–62. [14] Vassilev TI. Fair interpolation and approximation of B-splines by energy minimization and points insertion. Comput Aided Des 1996;28(9):753–60. [15] Wang W, Pottmann H, Liu Y. Fitting B-spline curves to point clouds by curvature-based squared distance minimization. ACM Trans Graphics 2006;25(2):214–38, http://dx.doi.org/10.1145/1138450.1138453. [16] Maekawa T, Matsumoto Y, Namiki K. Interpolation by geometric algorithm. Comput Aided Des 2007;39:313–32. [17] Chen Z, Luo X, Tan L, Ye B, Chen J. Progressive interpolation based on Catmull–Clark subdivision surfaces. Comput Graphics Forum 2008;27(7): 1823–7. [18] Gofuku S-i, Tamura S, Maekawa T. Point-tangent/point-normal B-spline curve interpolation by geometric algorithms. Comput Aided Des 2009;41:412–22. [19] Cheng F, Fan F, Lai S, Wang J, Yong J, Huang C. Progressive interpolation using Loop subdivision surfaces. In: Chen F, Juttler B, editors. GMP2008, vol. 4975. LNCS; 2008. p. 526–33. [20] Lin H, Wang G, Dong C. Constructing iterative non-uniform B-spline curve and surface to fit data points. Sci China Ser F 2004;47(3):315–31. [21] Lin H. The convergence of the geometric inter-polation algorithm. Comput Aided Des 2010;42(6):505–8. ˜ a JM. Progressive iteration approximation and the [22] Carnicer JM, Delgado J, Pen geometric algorithm. Comput Aided Des 2012;44(2):143–5. [23] Xiong Y, Li G, Mao A. Convergence of geometric interpolation using uniform b-splines. In: Proceedings of the computer-aided design and computer graphics. Jinan, China; 2011. p. 229–34.
Contents lists available at SciVerse ScienceDirect
Computers & Graphics journal homepage: www.elsevier.com/locate/cag
Computer Graphics in China
Convergence analysis for B-spline geometric interpolation Yunhui Xiong a, Guiqing Li b,n, Aihua Mao b a b
School of Science, South China University of Technology, Guangzhou, China School of Computer Science & Engineering, South China University of Technology, Guangzhou, China
a r t i c l e i n f o
a b s t r a c t
Article history: Received 10 December 2011 Received in revised form 10 July 2012 Accepted 10 July 2012 Available online 19 August 2012
In this paper, we propose a sufficient condition for the convergence of a geometric algorithm for interpolating a given polygon using non-uniform cubic B-splines. Geometric interpolation uses the given polygon as the initial shape of the control polygon of the B-spline and reduces the approximate error by iteratively updating the control points with the deviations from the corresponding interpolated vertices to their nearest footpoints on the current B-spline curve. The convergence condition is derived by employing a spectral radius estimation technique. The primary goal is to find for each control point a parametric interval within which the nearest footpoint should be confined such that the spectral radius of the error iteration matrix is smaller than 1. A convergent condition for the geometric interpolation of uniform B-splines can be derived as a special case of the new scheme. & 2012 Elsevier Ltd. All rights reserved.
Keywords: Geometric interpolation Convergence analysis Non-uniform B-splines Spectral analysis
1. Introduction The construction of a parametric curve passing through a set of scattered points is a fundamental task in computer aided geometric design due to its extensive applications in a variety of fields, such as economics, medicine, computer animation, manufacturing, and signal processing (see [1–6] for example). Though one may impose different requirements on generated curves depending on the application [1,7], the interpolatory curves should typically look pleasing and be as fair as possible. Interpolation problem based on B-spline curves can be formally described as follows: given a closed polygon P ¼ ðp1 ,p2 ,. . .,pn Þ and B-spline representation CðtÞ ¼
n X
Bi ðtÞp0i ,
i¼1
where Bi(t) and p0i are the non-uniform B-spline basis functions and the unknown control points, respectively, we need to determine p0i and t 1 o t 2 o ot n such that Cðt i Þ ¼ pi ,
i ¼ 1,2, . . ., n:
ð1Þ
Eq. (1) is actually a nonlinear under-constrained system with many solutions. In the literature, many approaches to obtaining an optimal solution are described. A straightforward method for obtaining an optimal solution is to specify {ti} explicitly and reduce Eq. (1) to a linear system. This n
Corresponding author. E-mail addresses: [email protected] (Y. Xiong), [email protected] (G. Li), [email protected] (A. Mao). 0097-8493/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.cag.2012.07.002
method usually leads to results with overshooting and undulating artifacts. Many parameterization methods have been proposed to remedy these problems [8–10]. Among these, arc-length parameterization and its variants are the most popular in practice. However, no solution is the best for all cases [11–13]. Some authors consider addressing the problem from the viewpoint of variation, but this method usually results in complicated optimizations [14]. In [11], Hoschek proposed an intrinsic parameterization to attack the curve fitting problem. Using this method, parameter values are computed implicitly by finding the nearest point from the data points to the fitting curve and are then employed to fit a new curve. The point nearest to a data point is called its footpoint. This approach iteratively refines the solution by repeating the process. Saux and Daniel [12] improved the efficiency of the approach, and Wang et al. [15] extended the approach with several new fitting errors. Unfortunately, these approaches degenerate if applied to interpolating problems instead of fitting tasks. Recently, Maekawa et al. [16] exploited a geometric interpolation algorithm to interpolate the vertices of a given mesh using subdivision surfaces. Their algorithm iteratively modifies each control vertex of the mesh using a deviation from the footpoint to its corresponding data point. In [16], two types of footpoints are discussed: the nearest point on the subdivision surface to the data point and the limit position of the corresponding control point on the surface to the data point (parametric distance deviation). As neither a large linear system nor the minimization of smooth energies is required, this approach is quite efficient and can generally produce high-quality results using the first type of correction due to its essence of intrinsic parameterization. In [17], Chen et al. extended this algorithm to the Catumll– Clark subdivision surfaces. Gofuku et al. [18] further studied the interpolation of tangents and normals using the geometric algorithm.
Y. Xiong et al. / Computers & Graphics 36 (2012) 884–891
However, the convergence analysis of the geometric interpolation algorithm remains open except in the special case of quadratic B-spline curves with the second type of deviation [15]. Cheng et al. showed that the algorithm is also convergent for Loop surfaces with the parametric distance deviation (the second type of deviation) [19]. Lin et al. introduced a class of iterative nonuniform B-spline curves [20] and presented a sufficient condition for the convergence of the geometric interpolation algorithm [21]. Carnicer et al. [22] further improved Lin’s proof. However, the convergence condition is too high level to derive a practical convergent scheme. In [23], Xiong et al. explored a sufficient condition that can guarantee the convergence of geometric interpolation based on a uniform B-spline. The limitation is that this condition does not work when using non-uniform B-splines due to the complexity of non-uniform B-spline basis functions. As an extended version of reference [23], this paper considers the case of non-uniform cubic B-spline curve interpolation. Our approach is also applicable to B-splines of degree 2–6 in which the maximum of the basis functions of the B-splines is greater than 0.5. The main contribution of this paper is the exploration of a sufficient condition that guarantees the convergence of the geometric interpolation algorithm. The convergent conditions for uniform B-spline geometric interpolation described in [23] is then derived as a special case of our new formulation. The remainder of this paper is organized as follows. Section 2 explains the mechanism of geometric interpolation using a simple case where only one control point of a given polygon is interpolated. Section 3 describes the proposed scheme in detail and then analyzes its convergence. In Section 4, several examples are given to examine the quality of interpolating curves using our new scheme. Finally, we conclude our paper in Section 5.
885
m Step 3 Vertex pm k in P is modified as follows
þ1 m pm ¼ pm k þ ek , k m m where the correction term em k ¼ pk C ðt k Þ . Step 4 Set m’m þ1, and repeat steps 2 and 3 until the correction term is less than a given threshold or the iterative step greater than a specified value. 0 Step 5 p0k ¼ pm k , and return P .
In Algorithm 1, the footpoint of pk on the current curve Cm(t) is required. Selecting nearest points as footpoints may yield a polynomial equation. If more than one nearest point is found, we choose the point whose parameter has the nearest value to uk. Then, we have the following results. Theorem 1. The convergence of Algorithm 1 can be guaranteed for a footpoint whose parametric value is constrained within ui1 r t m k r ui þ 1 , i.e., lim Hd ð pk ,C m ðtÞÞ ¼ 0, m-1
where Hd(A, B) is the Hausdorff distance from A to B. Proof. According to Algorithm 1, we have þ1 C m þ 1 ðtÞ ¼ ðpm pk ÞBk ðtÞ þ k
n X
pi Bi ðtÞ
i¼1 m ¼ ½pk C m ðt m k ÞBk ðtÞ þC ðtÞ, þ1 where {Bi(t)} are cubic B-spline basis functions. As C m þ 1 ðt m Þ is k mþ1 the nearest point to pk on the segment ui 1 rt rui þ 1 of C (t), þ1 :pk C m þ 1 ðt m Þ: r:pk C m þ 1 ðt m k Þ: k
2. Geometric interpolation to one control point with B-spline curves Given a polygon P ¼ ðp1 ,p2 ,. . .,pn Þ, we consider such a problem that constructs a non-uniform cubic B-spline C(t) defined by control polygon P 0 ¼ ðp1 , p2 ,. . .,p0k ,. . .,pn Þ and knot vector U ¼ ðu1 ,u0 ,u1 ,u2 ,. . .,uk ,. . .,un ,un þ 1 ,un þ 2 Þ, u1 ¼ un1 un ,u0 ¼ 0,un þ 1 ¼ un þ u1 , un þ 2 ¼ un þu2 , such that C(t) passes through the k-th vertex pk of polygon P. We modify the mechanism described in [23] slightly to address this single point interpolation problem using non-uniform B-splines. The new scheme also iteratively updates the position of the k-th control point and sustains other vertices unchanged. (See the pseudocode in Algorithm 1.)
m m m ¼ :pk ½pk C m ðt m k ÞBk ðt k ÞC ðt k Þ: m m ¼ :1Bk ðt m k Þ::pk C ðt k Þ:
¼
m Y
91Bk ðt jk Þ9:pk C 0 ðt 0k Þ::
j¼0
For
t jk A [ui 1,ui þ 1],
we can always guarantee 91Bk ðt jk Þ9 o1.
&
Remark 1. The convergent rate of Algorithm 1 is 5/6 for uniform cubic B-splines because 1/6rBk(t) r4/6 in this case. Remark 2. If the footpoint is sought on the whole Cm(t), Algorithm 1 is not guaranteed to converge [23]. 3. Geometric interpolation to closed polygons with B-spline curves
Algorithm 1. Constructing a B-spline curve interpolating a specified control point of the initial polygon Input
In this section, we consider solving Eq. (1) based on nonuniform cubic B-spline curves. Again, let
P ¼ ðp1 ,p2 ,. . .,pn Þ,
u1 ¼ un1 un ,u0 ¼ 0,un þ 1 ¼ un þu1 ,un þ 2 ¼ un þ u2 ,
U ¼ ðu1 ,u0 ,u1 ,u2 ,. . .,uk ,. . .,un ,un þ 1 ,un þ 2 Þ,
be the knot vector. Our geometric interpolation scheme works as described in Algorithm 2.
u1 ¼ un1 un ,u0 ¼ 0,un þ 1 ¼ un þ u1 ,un þ 2 ¼ un þ u2 : Output P 0 ¼ ðp1 ,p2 ,. . .,p0k ,. . .,pn Þ:
U ¼ ðu1 ,u0 ,u1 ,u2 ,. . .,uk ,. . .,un ,un þ 1 ,un þ 2 Þ,
Algorithm 2. The construction of the non-uniform B-spline curve interpolating all vertices of a given polygon. Input P ¼ ðp1 ,p2 ,. . .,pn Þ,
Step 1 Start from the given polygon: P 0 ¼P. Step 2 Assume we have gained P m ¼ ðp1 , ,pm k , ,pn Þ . Find the footpoint of pk on the curve Cm(t), which stands for the B-spline with control polygon Pm and knot vector U.
U ¼ ðu1 ,u0 ,u1 ,u2 ,. . .,uk ,. . .,un ,un þ 1 ,un þ 2 Þ, u1 ¼ un1 un ,u0 ¼ 0,un þ 1 ¼ un þu1 ,un þ 2 ¼ un þ u2 :
886
Y. Xiong et al. / Computers & Graphics 36 (2012) 884–891
Assume that Bi(u)(ui 1 otoui þ 1) achieves maximum at parameter value ti,max, namely, Bi ðt i,max Þ ¼ max Bi ðuÞ : t A ðui1 ,ui þ 1 Þ :
Output 0
P ¼
ðp01 ,p02 ,
,p0n Þ
0
Step 1 Start from the given polygon as the initial guess of P , namely, let P0 ¼P. m m Step 2 Assume we have obtained P m ¼ ðpm 1 ,p2 ,. . .,pn Þ: The control points of Pm þ 1 are calculated as follows m m þ1 pm ¼ pm i þ ðpi C ðt i ÞÞ i m
i ¼ 1,2,. . .,n
ð2Þ m
ðt m i Þ
where C is the footpoint of pi on C (t) within a given interval. (See Lemma 2.) Step 3 Set m’m þ1, and repeat Step 2 until specified thresholds are satisfied. 0 0 Step 4 P ¼Pm, and return P . From (2), we have C m þ 1 ðtÞ ¼
n X
þ1 pm Bj ðtÞ j
j¼1
¼
n X
m m ½pm j þðpj C ðt j ÞÞBj ðtÞ
j¼1
¼ C m ðtÞ þ
n X
½pj C m ðt m j ÞBj ðtÞ:
j¼1
It follows that n X
m ½pj C m ðt m j ÞBj ðt i Þ
m m ¼ ½B1 ðt m i Þ, ,Bi1 ðt i Þ,1Bi ðt i Þ, m m Bi þ 1 ðt m i Þ,. . .,Bn ðt i ÞD ,
i ¼ 1,2,. . .,n,
ð3Þ
where Dm is a column vector consisting of the correction terms m m m m T Dm ¼ ½p1 C m ðt m 1 Þ,p2 C ðt 2 Þ,. . .,pn C ðt n Þ
Eq. (3) has the following matrix form ¼ ðIBm ÞDm ,
D
Proof. We only need to show that there exists a real number e0A(0, 1) such that r(I Bm) r1 e0. As rðIBm Þ r:IBm : holds for arbitrary norm :: , our task is to find a concrete norm such that :IBm :r 1e0 holds under the given condition. Consider the row-sum norm of I–Bm 9 8 = < X :IBm :1 ¼ max 1Bi ðt m Bj ðt m i Þþ i Þ , ; i ¼ 1,...,n: 1 r j r n,j a i
where mþ1
Proof. The continuity of the basis functions implies that Bi ðtÞ Z 1=2 þc holds for some open interval containing ti,max . We can compute di,1 and di,2 by solving the polynomial equation Bi ðtÞ1=2c ¼ 0. & Bi ðt i,max Þ 41=2 is not always true. Fortunately, this is only the case for certain extreme situations in which the knot intervals are very non-uniform. To demonstrate this, we randomly generate 10,000 five-dimension knot vectors and then check the corresponding basis functions one by one for their maximum. Fig. 1 shows that only 25 examples include basis functions whose maxima (tagged with ‘ ’) are not greater than 0.5. For each unsuccessful case, we calculate the ratio max uj uj1 , i1r j r iþ 2 , min uj uj1 , i1 r j ri þ 2
Theorem 2. Algorithm 2 is convergent if the parametric value t m i of the footpoint of pi is confined within the interval [di,1,di,2] derived in Lemma 2 for i¼1, 2, y, n.
j¼1
mþ1
Lemma 2. For non-uniform cubic B-splines, there exist two real numbers di,1 and di,2 that satisfy ui1 r di,1 r t i,max r di,2 r ui þ 1 such that Bi ðtÞ Z 1=2 þc for an arbitrary tA[di,1,di,2], i¼1, 2, y, n.
and find that the ratio is greater than 16 for all these cases.
m m pi C m þ 1 ðt m i Þ ¼ pi C ðt i Þ
D
Let c be a user specified number satisfying Bi ðui,max Þ Z1=2 þ c if Bi ðui,max Þ 4 1=2 . We then have
¼ ½p1 C
mþ1
ðt m 1 Þ,p2 C
mþ1
ðt m 2 Þ,. . .,pn C
mþ1
T ðt m n Þ
and m
which implies that we should guarantee X 1Bi ðt m Bj ðt m i ¼ 1,2,. . .,n: i Þþ i Þ r 1e0 , 1 r j r n,j a i
m
Bm ¼ ðbij Þnn with bij ¼ Bj ðt m i Þ, Noting that
Noting the affine property of B-spline basis functions, the above inequality is equivalent to
þ1 :pj C m þ 1 ðt m Þ:r :pj C m þ 1 ðt m j Þ:, j
2ð1Bi ðt m i ÞÞ r1e0 : Namely,
:Dm þ 1 :r :D
mþ1
: r:IBm ::Dm : r r
m Y
:IBk ::D1 :,
k¼1
and the composite matrix following lemma.
m Q
:IBk : This leads to the
k¼1
Lemma 1. If there exist a positive constant e0 and an integer m0 such that the spectral radius r(I Bm) of I Bm satisfies r(I Bm)r1 e0 for all m4m0, then Algorithm 2 converges. Proof. It is obvious that lim :Dm þ 1 : ¼ 0 if the condition in m-1 Lemma 1 is satisfied. & Lin et al. derived a special case of Lemma 1 in [21]. This special case is not applicable yet because it does not address how to update Algorithm 2 such that r(I Bm)r1 e0 holds.
Bi ðt m i ÞZ
1 þ e0 : 2
ð4Þ
Set e0 ¼2c. According to Lemma 2, if tm i A ½di,1 , di,2 , then ð1 þ e0 Þ=2 ¼ 1=2 þc rBi ðt m & i Þ .
1 0.8 0.6 0.4 1 625 1249 1873 2497 3121 3745 4369 4993 5617 6241 6865 7489 8113 8737 9361 9985
we obtain
Fig. 1. The maximum distribution of the 10,000 basis functions.
Y. Xiong et al. / Computers & Graphics 36 (2012) 884–891
Remark 3. If Bi ðui,max Þ r1=2 appears occasionally we just set tm i ¼ ui to make the algorithm more robust. Corollary 1. Algorithm 2 is convergent for uniform cubic B-splines if the parameters t m i of the footpoints are confined to interval ½di,1 , di,2 , i¼1, 2, y, n, where di,1 ¼ i0:466178 and di,2 ¼iþ0.466178. Proof. See reference [23] for the proof. & In each iterative step, Algorithm 2 updates the control points simultaneously. Using recently renewed control points to evaluate the later control points results in Algorithm 3 is sometimes more efficient as demonstrated in Section 4. Algorithm 3. The construction of B-spline curve interpolating all vertices of an initial polygon using the instant updating mode. Step 1 Start from the given mesh as the initial control mesh, i.e., P 0 ¼P. þ1 Step 2 At the (mþ 1)-th iteration, P m þ 1 ¼ ðp1m þ 1 ,pm , ,pnm þ 1 Þ 2 þ1 mþ1 is renewed as follows. Assuming that p1m þ 1 ,pm ,. . .,p have 2 i1 been calculated, we first introduce an intermediate B-spline C m,i ðtÞ ¼
i1 X j¼1
þ1 pm Bj ðtÞ þ j
n X
pm j Bj ðtÞ,
j¼i
m,i and then find the footpoint (denoted by C m,i ðt m (t). i Þ) of pi on C mþ1 pi is finally modified by
m,i m pim þ 1 ¼ pm ðt i ÞÞ, i þ ðpi C
i ¼ 1,2,. . .,n,
ð5Þ
We will show the convergence of Algorithm 3 in Section 4 using numerical examples. 4. Experimental results and discussion We perform our experiments on a PC with a 2.40 GHz Intel Xeon E5620 CPU, 4 GB of RAM, and the Windows 7 OS. Five variants of our algorithms (4 for Algorithms 2 and 1 for Algorithm 3) summarized in Table 1 are implemented using MATLAB R2008a for non-uniform cubic B-splines. The distinction among these variants is that each seeks footpoints within differently sized parametric ranges. The footpoints (i.e., the nearest points) are obtained by finding the zero points of the derivative of the squared distance from the data point to the current B-spline using the MATLAB function roots(). We also conduct three types of experimental comparison in this section. The first one is used to investigate the efficiency of the variants, the second is used to show the spectra of the composite matrices numerically, and the third is used to compare the quality of the generated curves. All of the nonuniform knot vectors in the experiments are created using the chord length parameterization. To observe the efficiency of the variants, we generate 11 groups of polygon data. Each group includes 200 control polygons and 2200 polygons in total. The polygons in the same group have the same vertex number (5, 10, 15, 20, 30, 40, 50, 60, 70, 80, and 90 vertices). The positions of all of the polygonal vertices are randomly generated in the square [0,5] [0,5]. The convergent Table 1 Five variants of the interpolating algorithm. Range of Variant Variant Variant Variant Variant
0 1 2 3 4
tm i
(ui 2,ui þ 2) (ui 1,ui þ 1) (di,1,di,2) ui (di,1,di,2)
Updating mode Algorithm Algorithm Algorithm Algorithm Algorithm
2 2 2 2 3
887
criterion is that the maximal distance from data points to their footpoints is less than 10 6 or the iteration steps exceed 1000.
4.1. Convergence observation As the search range of variant 0 is far beyond our sufficient condition, many cases are expected to diverge. The experiment demonstrates that variant 0 fails to converge in the numerical sense for almost all examples when the control polygon has more than 40 vertices (see Fig. 2). Though variant 1 is numerically convergent for all of the above data, we find some abortive cases when increasing the number of control vertices to 300. In these cases, the distance error exhibits fluctuation around a small value (about 10 4).
4.2. Efficiency analysis Because variant 1 seeks the footpoints in a larger range than variants 2 and 4, its convergent rate is slightly higher than the convergent rates of variants 2 and 4. Variant 3 interpolates data points with fixed parametric values and therefore has the lowest convergent speed among the four schemes. This accords with the experimental results depicted in Fig. 3 in which the average iteration step for the polygons with the same vertex number is shown for the four schemes. In addition, the iteration step increases slowly for the variants other than variant 3. Nevertheless, the computational time is not completely consistent with the above order. Because variant 3 does not need to seek the footpoint, which is the most time-consuming operation in the algorithm, it is more efficient than variants 1, 2 and 4 for converging to the specified threshold, as shown in Fig. 4. It also shows that the time complexity of variants 1, 2 and 4 is very close. Variant 1 is slightly faster than the other two because of the higher convergent rate.
120% 100% 80% 60% 40% 20% 0%
5 10 15 20 30 40 50 60 70 80 90
Fig. 2. Percentage of those examples in which variant 0 converges against all examples (vertical axis: percent; horizontal axis: vertex number).
350 300 250 200 150 100 50 0
variant1 variant2 variant3 variant4 5 10 15 20 30 40 50 60 70 80 90
Fig. 3. The iteration step comparison among the four variants. vertical axis: the average times of iterator; horizontal axis: the vertex number of polygons.
888
Y. Xiong et al. / Computers & Graphics 36 (2012) 884–891
4.3. Spectral analysis To show the behavior of the composite matrices numerically, we have computed the spectral radius of the composite matrices for each round of iteration for every randomly generated polygon. As expected, the spectral radius for variant 0 increases for most of the experiments (see Fig. 5(a)). All experiments for variants 1, 2 and 3 are numerically convergent (Fig. 5(b–d)), and the
6 5 4
variant1
3
variant2
2
variant3
1
variant4
0
5 10 15 20 30 40 50 60 70 80 90
Fig. 4. The average time cost comparison among the four variants. vertical axis: the average time cost (second); horizontal axis: vertex number.
14
x 108
convergent rate of variant 1 appears to be slightly faster than that of variant 2. Again, variant 3 demands the maximal steps of iterations. 4.4. Quality comparison Though the range for seeking footpoints is confined to a small interval, variant 2 can still produce high quality curves. Figs. 6 and 7 illustrate that the curvature maps of curves generated by variant 2 are very close to those of variant 1. Variant 3, which is actually a common interpolating technique, generally results in curves with the worst quality. The other two experiments demonstrate an application of the algorithm in Chinese character generation. A Chinese character ‘‘Dragon’’ (Fig. 8) in Huawen Xingkai font and a Chinese character ‘‘Good’’ in Huawen Lishu font (Fig. 9) are used as examples here. Their contour polygons are extracted using the MATLAB function ginput(). We then apply the proposed scheme to construct the spline curves interpolating the character contours. Different simple polygons in the same character are handled separately in our implementation. Again, the results of variants 1 and 2 look very nice (as shown in Figs. 8(c and d) and 9(c and d)), but the results of variant 3 exhibit overshooting artifacts in some corners
1 0.9
12
0.8
10
0.7 0.6
8
0.5
6
0.4 0.3
4
0.2
2 0
0.1
0
5
10
15
20
25
30
35
40
45
50
0
1
1
0.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0
10
20
30
40
50
60
70
80
90
100
0
0
10
20
30
40
50
60
70
80
90
100
0
10
20
30
40
50
60
70
80
90
100
Fig. 5. The behavior of the spectral radius for the composite matrices of the 2200 non-uniform cubic B-splines. (The maximal iteration step is 1000.) (a, b, c, d) respectively correspond to variants 0, 1, 2, and 3. In the images, the horizontal and vertical axes stand for the spectral radius and the iteration step, respectively, and each curve plots the radius of a B-spline as a function of the iteration step.
Y. Xiong et al. / Computers & Graphics 36 (2012) 884–891
889
Fig. 6. Curvature maps. Top row: (a) variant 1, (b) variant 2, (c) variant 3, (d) the polygon to be interpolated. Bottom row: zoomed in images corresponding to the part of images marked by the dotted rectangles in the top row.
Fig. 7. Curvature maps. Top row: (a) variant 1, (b) variant 2, (c) variant 3, (d) the polygon to be interpolated. Bottom row: zoomed in images corresponding to the part of images marked by the dotted rectangles in the top row.
Fig. 8. Interpolation generation of the Chinese character ‘Dragon’: (a) original character; (b) polygon approximation of the character contour; (c, d, e) (middle row): the interpolation results of variants 1, 2, and 3 respectively; bottom row: zoomed in images corresponding to the part of images marked by the dotted rectangles in the middle row.
890
Y. Xiong et al. / Computers & Graphics 36 (2012) 884–891
Fig. 9. Interpolation generation of the Chinese character ‘Good’: (a) original character; (b) polygon approximation of the character contour; (c, d, e) (middle row): the interpolation results of variants 1, 2, and 3 respectively; bottom row: zoomed in images corresponding to the part of images marked by the dotted rectangles in the middle row.
of the characters (as shown in the zoom-in regions shown in Figs. 8e and 9e). We can conclude from aforementioned experiments that variant 2 can produce high-quality interpolating curves and can guarantee convergence except in some extreme cases.
5. Conclusions Geometric interpolation, which is an iterative method to construct control polygons or meshes, comes from Maekawa et al. [16], in which the convergence of the iterative method was not addressed. In this paper, a sufficient condition of convergence, which constrains the footpoint search to a local range, is presented for non-uniform B-splines geometric interpolation. Experiments demonstrate that the local method can still produce high-quality interpolating curves while enhancing the convergence speed greatly compared to the global method using non-uniform cubic B-splines. As for future work, it is interesting to give a theoretical condition in terms of interval lengths to guarantee the maximum of non-uniform B-spline bases greater than 0.5. For convenience, our discussion only involves the interpolation of closed polygons.
Bi ðtÞ
¼
The convergence analysis of the control polyline interpolation must be supplied with clamped conditions to apply the approach. Our goal is to develop a similar convergent condition for the geometric interpolation of parametric surfaces and subdivision surfaces. Another interesting problem is to investigate the mechanism that makes this method capable of producing smoother results than interpolation schemes specifying parameter values.
Acknowledgments We are delighted to thank the reviewers for their valuable comments. The work is partially supported by NSFC (60973084, 61003173), NSF of Guangdong (915106410-1000106), and Fundamental Research Funds for the Central Universities (2012zz0063, 2012zm0062).
Appendix A. Piecewise representation of base function Bi(t)
8 ðtui2 Þ3 > > ðui þ 1 ui2 Þðui ui2 Þðui1 ui2 Þ > > h i > > ðui þ 1 tÞðtui1 Þ ðui þ 2 tÞ ðtui2 Þ ðtui2 Þðui tÞ ðtui1 Þ2 > > > ðui þ 1 u Þ ðui ui2 Þðui ui1 Þ þ ðui þ 1 ui1 Þðui ui1 Þ þ ðui þ 2 ui1 Þ ðui þ 1 ui1 Þðui ui1 Þ > i2 < h i ðui þ 1 tÞ2 ðui þ 2 tÞ ðtui1 Þðui þ 1 tÞ ðui þ 2 tÞðtui Þ ðtui2 Þ > ðui þ 1 ui2 Þ ðui þ 1 ui Þðui þ 1 ui1 Þ þ ðui þ 2 ui1 Þ ðui þ 1 ui1 Þðui þ 1 ui Þ þ ðui þ 2 ui Þðui þ 1 ui Þ > > > > ðui þ 2 tÞ3 > > > > > ðui þ 2 ui þ 1 Þðui þ 2 ui Þðui þ 2 ui1 Þ : 0
t A ðui2 ,ui1 Þ t A ðui1 ,ui Þ t A ðui ,ui þ 1 Þ t A ðui þ 1 ,ui þ 2 Þ others
Y. Xiong et al. / Computers & Graphics 36 (2012) 884–891
References [1] Hagan P, West G. Interpolation methods for curve construction. Appl Math Finance 2006;13(2):89–129. [2] Jacob M, Blu T, Unser M. Efficient energies and algorithms for parametric snakes. IEEE Trans Image Process 2004;13(9):1231–44. [3] Ramamoorthi R, Barr AH. Fast construction of accurate quaternion splines. In: computer graphics Proceedings, annual conference series, ACM SIGGRAPH. Los Angeles; 1997. p. 287–92. [4] Nam S-H, Yang M-Y. A study on a generalized parametric interpolatior with real-time jerk-limited acceleration. Comput Aided Des 2004;36(1):27–36. [5] Wang F-C, Yang D. Nearly arc-length parameterized quintic-spline interpolation for precision machining. Comput Aided Des 1993;25(3):281–8. [6] Levin A. Interpolating nets of curves by smooth subdivision surfaces. In: computer graphics Proceedings, annual conference series, ACM SIGGRAPH. Los Angeles; 1999. p. 57–64. [7] Piegl LA, Rajab K, Smarodzinana V. Curve interpolation with directional constraints for engineering design. Eng Comput 2008;24:79–85. [8] Floater MS, Surazhsky T. Parameterization for curve interpolation. Stud Comput Math 2006;12:39–54. [9] Papaioannou SG, Patrikoussakis MM. Curve interpolation based on the canonical arc length parametrization. Comput Aided Des 2011;43:21–30. [10] Kouibia A, Pasadas M, Rodrı´guez ML. Optimization of parameters for curve interpolation by cubic splines. J Comput Appl Math 2011http://dx.doi.org/ 10.1016/j.cam.2011.03.013. [11] Hoschek J. Intrinsic parametrization for approximation. Comput Aided Geom Des 1988;5(2):27–31. [12] Saux E, Daniel M. An improved Hoschek intrinsic parametrization. Comput Aided Geom Des 2003;20(8–9):513–21.
891
[13] Pottmann H, Leopoldseder S. A concept for parametric surface fitting which avoids the parametrization problem. Comput Aided Geom Des 2003;20: 243–62. [14] Vassilev TI. Fair interpolation and approximation of B-splines by energy minimization and points insertion. Comput Aided Des 1996;28(9):753–60. [15] Wang W, Pottmann H, Liu Y. Fitting B-spline curves to point clouds by curvature-based squared distance minimization. ACM Trans Graphics 2006;25(2):214–38, http://dx.doi.org/10.1145/1138450.1138453. [16] Maekawa T, Matsumoto Y, Namiki K. Interpolation by geometric algorithm. Comput Aided Des 2007;39:313–32. [17] Chen Z, Luo X, Tan L, Ye B, Chen J. Progressive interpolation based on Catmull–Clark subdivision surfaces. Comput Graphics Forum 2008;27(7): 1823–7. [18] Gofuku S-i, Tamura S, Maekawa T. Point-tangent/point-normal B-spline curve interpolation by geometric algorithms. Comput Aided Des 2009;41:412–22. [19] Cheng F, Fan F, Lai S, Wang J, Yong J, Huang C. Progressive interpolation using Loop subdivision surfaces. In: Chen F, Juttler B, editors. GMP2008, vol. 4975. LNCS; 2008. p. 526–33. [20] Lin H, Wang G, Dong C. Constructing iterative non-uniform B-spline curve and surface to fit data points. Sci China Ser F 2004;47(3):315–31. [21] Lin H. The convergence of the geometric inter-polation algorithm. Comput Aided Des 2010;42(6):505–8. ˜ a JM. Progressive iteration approximation and the [22] Carnicer JM, Delgado J, Pen geometric algorithm. Comput Aided Des 2012;44(2):143–5. [23] Xiong Y, Li G, Mao A. Convergence of geometric interpolation using uniform b-splines. In: Proceedings of the computer-aided design and computer graphics. Jinan, China; 2011. p. 229–34.