Interpolatory Subdivision Curves via Diffusion of Normals Yutaka Ohtake
Alexander Belyaev
Hans-Peter Seidel
Computer Graphics Group, Max-Planck-Institut f¨ur Informatik Stuhlsatzenhausweg 85, 66123 Saarbr¨ucken, Germany E-mails: fohtake,belyaev,
[email protected] Phone: +49 681 9325-408 Fax: +49 681 9325-499
Abstract
that our approach often produces to a better interpolation than uniform and nonuniform “4-point interpolatory subdivision” schemes [4, 5, 6].
In this paper, we propose a new interpolatory subdivision scheme for generating nice-looking curvature-continuous curves of round shapes. The scheme is based on a diffusion of normals. Given a subdivided polyline, the new polyline vertices inserted at the the splitting step are updated in order to fit diffused (averaged with appropriate weights) normals. Although the resulting interpolatory subdivision scheme is non-stationary, nonlinear, and nonuniform from the traditional point of view, the scheme is easy to implement because the same simple geometric procedure for generating new vertices is used at each subdivision step. According to our experiments, the scheme is robust and demonstrate very good convergence properties. Keywords: interpolatory subdivision, diffusion of normals.
Figure 1. Interpolatory subdivision urves produ es by a method developed in this paper. Left: a spiral. Right: from ve unevenly pla ed points to an approximation of a ir le with maximum radius error = 3.5%.
1 Introduction Discrete data interpolation dates back to ancient times [14]. Nevertheless until now solving a “simple” problem of a fast and robust interpolation of sparse point data on a plane by a smooth nice-looking curve remains to be an active research area because of numerous CAD/CAM applications. The interpolatory subdivision techniques developed recently [3, 4, 5, 10, 11, 12] constitute a very promising approach to the problem. Recently shape processing via diffusion of normals has become a promising research direction in the geometric modeling field [1, 7, 8, 9, 15, 16, 19, 20]. In this paper, we combine a simple subdivision technique with a diffusion of normals for generating nice-looking round-shaped interpolating curves. Fig. 1 demonstrates the power of our approach. In next sections, we will also show
Subdivision. The subdivision approach for shape design is based on modeling smooth shapes as a limit of successively refined polygonal shapes (see, for example, [18, 22] and references therein). Consider a polyline (: : : ; pj1 ; pj0 ; pj 1 ; : : :), where the upper index j enumerates the subdivision steps. The standard subdivision procedure consists of two steps: the splitting step and the averaging step. The splitting step inserts the midpoints in the segments of the polyline
q2 j
On a leave from University of Aizu, Aizu-Wakamatsu 965-8580 Japan
k
1
=
p
j k
;
q2 j
k+1
=
1 2
p
j k
+
p
j k+1
– apply local averaging to the normals and produce modified (smoothed) normals;
The averaging step computes a weighted average of several successive points
p +1 = j
X
ai qk+i
– update the positions of midpoints (odd vertices) by fitting the polyline to the field of modified normals.
j
k
i
This process is demonstrated in Fig. 3. In the classical sense the resulting interpolatory subdivision scheme is non-stationary, nonlinear, and nonuniform. However its practical implementation is easy because the scheme uses the same simple geometric procedure for generating new vertices at each subdivision step.
For example, if we consider the averaging mask (a 1 ; a0 ; a1 ) =
1 4
(1; 2; 1)
the polyline (: : : ; pj1 ; pj0 ; pj 1 ; : : :) converges to a cubic Bspline as j ! 1 [13].
2 Subdivion via Diffusion of Normals
Interpolatory subdivision. Linear interpolatory subdivision schemes are given by
(
p
j +1 k
=
q
j
P i
Basic Idea Consider an oriented polyline (: : : ; pj1 ; pj0 ; pj 1 ; : : :), where the upper index j enumerates the subdivision steps, and following the standard subdivision procedure insert the midpoints in the segments of the polyline
if k is even
k
ai qk+i j
if k is odd
The uniform 4-point interpolatory subdivision [4, 5] is defined by
q2 j
k
(a 2; a 1 ; a0 ; a1 ; a2 ) =
1 16
(
2; 5; 10; 5;
2):
=
p
j k
q2 j
;
k+1
1
=
2
p
j k
+
p
j
k+1
The even vertices remain at the same positions
p2+1 = q2 j
Fig. 2 demonstrates generating an interpolatory curve via the 4-point interpolatory subdivision.
j
k
k
:
The odd vertices are updated as follows. Let njk be the unit vector orthogonal to the segment qj2k qj2k+2 and consistent with the orientation of the polyline. Let us attach the normal vectors njk at the midpoints qj2k+1 . For the segments qj2k qj2k+1 and qj2k+1 qj2k+2 we compute unit vectors mj2k and mj2k+1 , respectively, using an averaging procedure, see the next section and left image of Fig. 4 for details. Then the position for the vertex pj2+1 is determined such that k+1 j j the vectors m2k and m2k+1 approximate the normals of the segments pj2+1 qj and pj2+1 qj , respectively, in a k+1 2k k+1 2k+2 certain least square sense.
Figure 2. From left to right: generating an interpolatory urve by subdivision. The uniform 4-point interpolatory subdivision s heme is used. The subdivision schemes described above are linear because new positions of vertices are obtained by linear averaging, uniform since the same averaging mask is used everywhere along the curve, and stationary because the same mask is used for each iteration of the subdivision process. More complex, nonuniform and non-stationary, subdivision schemes can be found, for example, in [21].
Averaging Normals Denote by (v) 2 [ ; ℄ the angle between a vector v and a prescribed fixed direction. We compute the vectors mj2k and mj2k+1 by linear averaging of angles
Basic idea of the proposed scheme. Our approach is based on averaging certain the polyline normals instead of averaging the polyline vertex positions. Our basic idea is conceptually simple. One iteration of our subdivision scheme consists of the following steps. Given a polyline, we
(mj2k )
(m2k+1 ) j
= =
w1 (njk ) + w2 (njk 1 ) j j w3 (nk ) + w4 (nk+1 );
(1)
where wi are weights: w1 + w2 = w3 + w4 = 1. We choose the weights equal to those used in the Chaikin corner cutting procedure:
– double the number of vertices by inserting the midpoints;
(w1 ; w2 ) = (w3 ; w4 ) =
2
1 4
(3; 1):
(2)
Figure 3. From left to right by zigzag: inserting the midpoints, averaging normals, updating positions of the inserted midpoints a
ording to new normals, and so on. The s heme qui kly produ es ne polygonal approximations of a smooth urve interpolating the initial polyline verti es.
w2
n k-1
m2k n k
q 2k
Let us consider pj2+1 as unknown x and let n1 and n2 k+1 be the unit normals of the segments xqj2k and xqj2k+2 , respectively. An error function E (x) measuring how close the pairs n1 ; n2 and mj2k ; mj2k+1 are can be defined as
w3
w1
q2k+1
m2k+1
q 2k+2
w4
m2k
n k+1
nk
E (x) = l1
n k-1 w2
kn
1
m2 j
k
k
2
+ l2
kn
m2 j
2
k+1
( n k , n k-1 )
p2+1+1 = argmin E (x): j
x
k
The gradient of error function E (x) can be easily com-
!
k
j
k
= =
R w2 6 (nk ; nk j
j
n
R w4 6 (nk ; nk+1 ) j
j
1)
j
k
n
j k
!
puted. Let us denote vectors qj2k x and qj2k+2 x by d1 and d2 , respectively, and let vector Pm (d) be the projection of d onto the perpendicular to m. Then
In practice, we implement (1) by j
2
where l1 and l2 are the lengths of segments xqj2k and xqj2k+2 , respectively, as seen in the left image of Fig. 5. is determined by Now pj2+1 k+1
Figure 4. The ve tors m2k and m2k+1 are omputed by averaging normals nk 1 , nk , and nk+1 .
8 < m2 : m2 +1
k;
rE (x)
;
where R(') is the 2 2 matrix representing the rotation in the counter-clockwise direction on the angle ', and 6 (v1 ; v2 ) 2 [ ; ℄ is the angle between two unit vectors v1 to v2 measured by
=
2
!
m2 (d1 ) kPm2 (d1 )k ! (d ) Pm d2 2 +1 2 kd2 k kPm2 +1 (d2 )k ;
d1 kd1 k
P
j
k
j
k
j
+2
k
j
k
The right image of Fig. 5 explains a geometric meaning of rE (x). Now E (x) can be effectively minimized. In our current implementation, the conjugate gradient method [17] is used. To obtain more evenly sampled curves the initial guess is set to the midpoint qj2k+1 of the segment.
6 (v1 ; v2 ) = sign (det(v1 ; v2 )) ar
os(v1 v2 ) Fig. 4 demonstrates the above averaging procedure.
Updating Positions of Odd Vertices
Dealing with Open Polylines
Now we want to place the vertex p2k+1 such that the vectors mj2k and mj2k+1 approximate the normals of the qj and pj2+1 qj , respectively. segments pj2+1 k+1 2k k+1 2k+2 j +1
So far we have considered only the closed curve case. Open curves can be treated in a standard way. Given an initial polyline, first we put two dummy points on the extended 3
m2k
E(x)
m2k
n1
n2
m2k+1
Pm2k(d1 )
x q 2k
l1
l2
q 2k+2
q 2k
d1
x
Figure 5. Left: error fun tion E (x). Right: geometri meaning of rE (x).
π θ(n)
π θ(n)
0
end segments of the polyline. Next the dummy points are connected and a closed polyline is generated. Then our interpolatory subdivision scheme is applied. Finally, the dummy part is removed. Fig. 6 demonstrates this process.
−π
Figure 7. Top: two urves interpolating the same ontrol points and generated with the uniform weights (2) (left image) and non-uniform weights (3) (right image). Bottom: plots of (n(s)), the angle between the urve normal n(s) and a xed dire tion, versus ar length s for both the urves. θ(m2k ) θ(n k )
q 2k-2
(w3 ; w4 )
=
1 ) (ek + 2 ek 1 ; ek ) +1 ) (ek + 2 ek+1 ; ek );
w4 : w3
e k-1
q 2k
ek
q 2k+2
e k+1
q 2k+4
Figure 8. Nonuniform weighting (3) allows us to generate urvature- ontinuous urves.
The subdivision scheme developed in the three previous sections may produce a curve with curvature discontinuities, as shown in the left image of Fig. 7. To generate curvature-continuous curves we use a non-uniform corner cutting scheme [6] instead of Chaikin’s uniform weights (2). That new non-uniform weights lead to generating curvature-continuous curves, as it is demonstrated in the right image of Fig. 7. Let ek denote the length of the segment qj2k qj2k+2 . Following [6] we use the following weights =
θ(n k+1 )
arclength parameter
Non-uniform Weighting for Curvature Continuity
(w1 ; w2 )
θ(m2k+2 )
w1 : w2
Figure 6. Left: losed subdivision urve produ ed by
losed polyline. Middle: open urve subdivision urve is generated from the same ontrol points. Right: inserting dummy points redu es the losed urve ase to the open one.
1 2(ek +ek 1 2(ek +ek
s
−π
θ(n k-1 )
(
0
s
Let us a control point where we want to have a corner with a pair of parameters B ; F 2 [0; 1℄, where subscripts B and F are used for “Backward” and “Forward”, respectively. To compute the averaged normals associated with the segments containing the control points where corners are created, the following special weights are used.
(w b1 ; wb2 )
(w b3 ; wb4 )
= =
F (1; 0) + (1
B (1; 0) + (1
F ) (w1 ; w2 )
B ) (w3 ; w4 )
(4)
Using these weights produces sharp corners if F;B = 1 and creates rounded parts if F;B = 0. Fig. 9 demonstrates how F;B = 0 control sharpness of the subdivision curve at control points.
(3)
in (1). Fig. 8 presents the geometric meaning of the above non-uniform weights.
Generating Corners
3 Discussion
In order to be able to generate subdivision curves with corners we adapt to an idea presented in [2].
We have presented a new interpolatory subdivision scheme for generating nice-looking curvature-continuous 4
Figure 9. Using (4) with various F;B allows us to reate urves with/without orners from the same ontrol points. The numbers are values of F;B used to generate the urves (zero values are not shown).
Figure 11. Curves interpolating ontrol points spa ed irregularly. Top-left: the uniform 4-point subdivision s hemes was used. Top-right: the non-uniform 4-point subdivision s heme was used. Bottom: the method developed in the paper was applied. Bottomleft: uniform weights (2) were used. Bottom-right: non-uniform weights (3) were used.
curves of round shapes. A simple Java applet demonstrating the scheme can be found at 1 . Our method is very useful if there is a need to generate curves containing parts approximating circular arcs, since it is difficult using the classical 4-point interpolatory subdivision schemes [4, 5, 6]. Fig. 10 shows advantages of our subdivision method in generating circular shapes.
of the control polygon intersects a segment of the polygon, as seen in Fig. 12.
Figure 10. The subdivision urves interpolating the verti es of a regular triangle. Left: the 4-point subdivision s heme was used. Right: the subdivision s heme developed in this paper (uniform weight) was used; the maximum radius error = 0.34%.
Figure 12. A subdivision urve generated a
ording to the s heme des ribed in the paper bifur ates when a vertex of the ontrol polygon rosses a segment of the polygon. The idea of using normals for subdivision purposes seems to be very promising. In future we hope to extend our approach to 3D.
Fig. 11 demonstrate robustness of our subdivision scheme against uneven distribution of control points. Our research on subdivision via averaging normals is at the beginning stage. A rigorous mathematical treatment of the proposed interpolatory subdivision scheme is needed. We have also noticed a small drawback of our current implementation: the subdivision curve bifurcates when a vertex
References [1] A. G. Belyaev, Yu. Ohtake, and K. Abe. Detection of ridges and ravines on range images and triangular meshes. In Vision Geometry IX, Proc. SPIE 4117, pages 146–154, 2000.
1 http://www.mpi-sb.mpg.de/˜ohtake/subdiv
5
[2] T. DeRose, M. Kass, and T. Truong. Subdivision surfaces in character animation. In Proceedings of ACM SIGGRAPH 1998, pages 85–94, 1998. [3] G. Deslauriers and S. Dubic. Symmetric iterative interpolation processes. Constr. Approx., 5:49–68, 1989. [4] S. Dubic. Interpolation through an iterative scheme. Journal of Mathematical Analysis and Applications, 114:185–204, 1986. [5] N. Dyn, D. Levin, and J. Gregory. A 4-point interpolatory subdivision scheme for curve design. Computer-Aided Geometric Design, 4:257–268, 1987. [6] J. Gregory and R. Qu. Non-uniform corner cutting. Computer-Aided Geometric Design, 13(8):763–772, 1996. [7] S. Karbacher and G. H¨ausler. New approach for modeling and smoothing of scattered 3D data. In Three-Dimensional Image Capture and Applications, Proc. SPIE 3313, pages 115–125, San Jose, California, January 1998. [8] S. Karbacher, S. Seeger, and G. H¨ausler. A non-linear subdivision scheme for triangle meshes. In Vision, Modeling and Visualization 2000, pages 163–170, Saarbr¨ucken, Germany, November 2000. [9] S. Karbacher, S. Seeger, and G. H¨ausler. Refining triangle meshes by non-linear subdivision. In Proceedings of Third International Conference on 3-D Digital Imaging amd Modeling, pages 270–277. IEEE Computer Society, 2001. [10] L. Kobbelt. Iterative Erzeugung glatter Interpolaten. PhD thesis, Universit¨at Karlsruhe, 1994. [11] L. Kobbelt. A variational approach to subdivision. Computer Aided Geometric Design, 13:743–761, 1996. [12] L. Kobbelt and P. Schr¨oder. A multiresolution framework for variational subdivision. ACM Transactions on Graphics, 17(4):209–237, 1998. [13] J. M. Lane and Riesenfeld R. F. A theoretical development for the computer generation and display of piecewise polynomial surfaces. IEEE Trans. Pattern Analysis and Machine Intelligence, 2(1):35–46, 1980. [14] E. Meijering. A chronology of interpolation: From ancient astronomy to modern signal and image processing. Proceedings of the IEEE, 90:319–342, 2002. [15] Yu. Ohtake, A. G. Belyaev, and I. A. Bogaevski. Mesh regularization and adaptive smoothing. Computer-Aided Design, 33(11):789–800, 2001. [16] Yu. Ohtake, A. G. Belyaev, and H.-P. Seidel. Mesh smoothing by adaptive and anisotropic Gaussian filter. In Vision, Modeling, and Visualization 2002, pages 203–210, Erlangen, Germany, November 2002. [17] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipies in C: The Art of Scientific Computing. Cambridge University Press, 1993. [18] E. J. Stollnitz, T. D. DeRose, and D. H. Salesin. Wavelets for Computer Graphics. Morgan Kaufmann Publishers, 1996. [19] T. Tasdizen, R. Whitaker, P. Burchard, and S. Osher. Geometric surface smoothing via anisotropic diffusion of normals. In Proceedings of IEEE Visualization 2002, pages 125–132, November 2002. [20] G. Taubin. Linear anisotropic mesh filtering. IBM Research Report RC22213 (W0110-051), IBM, October 2001. [21] J. Warren and H. Weimer. Subdivision Methods for Geometric Design. Academic Press, 2002. [22] D. Zorin, P. Schr¨oder, T. DeRose, L. Kobbelt, A. Levin, and W. Sweldens. Subdivision for modeling and animation. In SIGGRAPH 2000 Course Notes, July 2000.
6