Surface reconstruction via geodesic interpolation - Semantic Scholar
Computer-Aided Design 40 (2008) 480–492 www.elsevier.com/locate/cad
Surface reconstruction via geodesic interpolation N. Sprynski b , N. Szafran a , B. Lacolle a , L. Biard a,∗ a Laboratoire Jean Kuntzmann, Universit´e Joseph Fourier - Grenoble, France b CEA MINATEC, LETI/DCIS/SMOC, Grenoble, France
Received 24 May 2007; accepted 10 January 2008
Abstract This paper is concerned with reconstruction of numerical or real surfaces based on the knowledge of some geodesic curves on the surface. So, considering two regular 3D-curves f 0 (t) and f 1 (t), our purpose is to construct a surface which interpolates these two curves in such a way that these two curves are geodesics on this surface. This will be accomplished using Hermite interpolation. For a real surface, it will be shown that geodesics can be acquired using a ribbon of micro-sensors. c 2008 Elsevier Ltd. All rights reserved.
Keywords: Geodesic curves; Surface reconstruction; Hermite interpolation; Curvature minimization; Motion capture; Embedded micro-sensors
1. Introduction This paper is concerned with reconstruction of surfaces from curves running on this surface. Thus, given arbitrary curves on a surface, Coons methods [1,2,5,4] allow us to get a surface which interpolates the surface. Notice that Coons processes are mostly used to create/define 3D objects from some spatial curves. These 3D curves are designed/computed first. They provide the main characteristic features of the desired surface and model the essential shapes of the object. The Coons process, is then performed as a “filling” method to get the final surface. So, the surface is essentially modeled by the initial spatial curves. Thus, concerning a real 3D surface, the immediate question is how to get lines from this surface. Notice that we are concerned by the reconstruction problem of a real surface for which no parameterization is accessible. Then, in order to get a consistent reconstruction process, these lines must represent intrinsic characteristic features of the surface. Fortunately we are able to provide a positive answer to these two points. A natural way would be to acquire geodesic curves running on the surface. That can be achieved employing a ribbon of embedded micro-sensors (see Figs. 14 and 17) and it will be shown, in Section 5 that such a ribbon actually follows geodesic ∗ Corresponding author.
E-mail address: [email protected] (L. Biard). c 2008 Elsevier Ltd. All rights reserved. 0010-4485/$ - see front matter doi:10.1016/j.cad.2008.01.005
curves on the surface. This ribbon is inextensible and not “flexible” and has the same behaviour as a thin strip of paper, i.e., as a developable surface. Its mechanical properties are modeled in Section 5. The micro-sensors are able to give their own orientation in space from which are deduced data about the local tangency of the ribbon. This technology has been developed at the CEA-LETI, a research center of technology specialized in micro-electronics. Thus, from the tangential information of the ribbon, we can deduce its 3D shape. This reconstruction process is described in [9] and works for planar and spatial curves. Thus, we are now concerned with reconstruction of surfaces based on the knowledge of some geodesic curves on the surface. So, considering two regular 3D-curves f 0 (t) and f 1 (t), our purpose is to construct a surface which interpolates these two curves in such a way that these two curves are geodesics on this surface. This problem of surface reconstruction using only tangential orientations without any position data seems to have not yet been handled and the need to solve such a problem has become important in the field of motion capture. Notice that this is not an envelope problem in which the tangent lines are fixed. Here we only know the directions of these tangent lines and, moreover, we have the added constraint that the length between two consecutive contact points are given. This paper is organized as follows: Section 2 deals with fundamental properties of curves on surfaces. A solution to the geodesic interpolation problem is then presented in Section 3
481
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
with examples and an analysis of the algorithm. In Section 4 we consider numerical applications. Precisely, we consider numerical models of surfaces on which exact geodesic curves are evaluated with an algorithm issued from [7,8]. Then, the surface is reconstructed from these geodesics. Finally, in Section 5 we deal with reconstruction of real surfaces using the ribbon of micro-sensors. In particular, we show that such a ribbon follows geodesic curves on the surface. 2. Background on curves and surfaces In the following, all curves and surfaces are considered to be sufficiently smooth (at least of class C 3 ). The parameters of functions will sometimes be omitted when no confusion can occur and the dot will denote the inner product. For a regular 3D curve ϕ : t 7→ ϕ(t) we will denote by ϕ : s 7→ ϕ(s) its normal parameterization by arc length with the following properties kϕ 0 (s)k = 1 and ϕ 0 (s) ⊥ ϕ 00 (s) for each value of s. • With each point ϕ(s) of such a curve such that ϕ 00 (s) 6= 0, we associate the Serret–Frenet frame (τ (s), n(s), b(s)) where τ (s) = ϕ 0 (s), n(s) = ϕ 00 (s)/kϕ 00 (t)k and b(s) = τ (s)×n(s) are respectively the unit tangent vector, the unit normal vector and the binormal vector of the curve at point ϕ(s). The unit normal vector n(s) is also called the main normal. Considering the arc length, the Serret–Frenet frame satisfies the Serret–Frenet relations: τ (s) 0 k(s) 0 τ (s) d n(s) = −k(s) 0 θ (s) n(s) , ds b(s) 0 −θ (s) 0 b(s) where k(s) and θ (s) are the curvature and the torsion of the curve at point ϕ(s) defined by: k(s) = kϕ 00 (s)k
θ (s) =
and
det(ϕ 0 (s), ϕ 00 (s), ϕ 000 (s)) . (1) kϕ 00 (s)k2
The osculating plane of the curve at point ϕ(s) is defined by the two vectors τ (s) and n(s) and does not depend on the parameterization of the curve. If for some parameter s, the curvature k(s) = 0, then ϕ 00 (s) = 0, so that the normal vector n(s) and the osculating plane can not be defined. Such a point is called an inflection point of the curve. • For a regular oriented surface (u, v) 7→ F(u, v), we consider the normal vector N (u, v) and the Gaussian curvature K T (u, v) at point F(u, v) defined by N (u, v) =
Fu0 × Fv0 , kFu0 × Fv0 k
K T (u, v) = K 1 K 2 =
and
L = N · Fu002 ,
00 M = N · Fuv ,
g
which define the normal curvature kn (s), the geodesic curvature k g (s) and the geodesic torsion θg (s) by kn (s) =
dτ · N, ds
k g (s) =
dτ · V, ds
θg (s) =
dN · V. ds
(4)
• A regular curve f (s), running on the surface F is a geodesic curve of this surface if its geodesic curvature k g is identically zero. Thus, from k g ≡ 0 we deduce with the Darboux relations that dτ ds is collinear with the vector N for each value of the parameter, so that the Frenet and Darboux frames agree regardless of the orientation. Thus, since the osculating plane of a curve is unchanged under reparameterization, we get the following useful characterization of geodesic curves. A regular curve t 7→ f (t) running on a surface F(u, v) is a geodesic of F if its osculating plane at each non-inflection point f (t) is normal to the tangent plane of the surface F at point F(u(t), v(t)) = f (t). In the case of an isolated inflection point of the curve, we may encounter situations, see e.g. [3], where the osculating plane can not be defined by continuity as the limit planes approaching from the left and from the right side will disagree. However, in the given situation considering geodesic curves, it turn out later that the geodesic’s Frenet frame and the Darboux frame will agree for isolated inflection points, because, for geodesics, the main normal is collinear with respective surface normal. Therefore the continuity of the Darboux frame implies the continuity of the osculating plane in the corresponding Frenet frame. In the case of non-isolated inflection points, the vanishing of the curvature on a domain implies that the associated points belong to a straight line (so, a geodesic curve) running on the surface. Their Darboux frame can be defined on each point whereas the Frenet frame can not be defined at any point.
(2)
where K 1 and K 2 are the principal curvatures and where F = Fu0 · Fv0 ,
n
3. Geodesic interpolation
L N − M2 , EG − F2
E = Fu0 · Fu0 ,
elliptic points (K T > 0), hyperbolic points (K T < 0) and parabolic points (K T = 0). • Consider a curve f (s) = F(u(s), v(s)) running on the surface F. With each point f (s) we associate the Darboux frame (τ (s), V (s), N (s)) where τ (s) is the unit tangent vector of the curve f , N (s) is the unit normal vector of the surface F at point F(u(s), v(s)) = f (s) and V (s) = N (s) × τ (s). Considering the arc length, the Darboux frame satisfies the following Darboux relations: τ (s) 0 k g (s) kn (s) τ (s) d V (s) = −k g (s) 0 −θg (s) V (s) (3) ds N (s) −k (s) θ (s) 0 N (s)
G = Fv0 · Fv0 ; N = N · Fv002 ,
are the coefficients of the first and second fundamental form of the surface. The Gaussian curvature allows us to distinguish
Consider now two regular 3D curves f 0 (u) and f 1 (u), u ∈ [0, 1], such that f 0 (u) 6= f 1 (u) for each parameter u, and remember that our goal is to construct a surface which interpolates these two curves in such a way that these curves are geodesics on this surface. Precisely, let us denote by F(u, v), with (u, v) ∈ [0, 1]2 the interpolating surface. The surface F has to fulfill the two following properties.
482
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
Fig. 1. Geodesic interpolation of the two 3D-curves f 0 (u) and f 1 (u).
• Interpolating property: F(u, 0) = f 0 (u) and F(u, 1) = f 1 (u) for each u ∈ [0, 1]. • Geodesic property: for i ∈ {0, 1}, the osculating plane of the curve f i at each point f i (u) is normal to the tangent plane of the surface F at point F(u, i). Notice that by the following construction, each one of the two (geodesic) curves f i (u) will be an isoparametric curve of the surface F. 3.1. Main steps of the algorithm The algorithm can be described by the following steps (see also Fig. 1). 1. Tangent plane: for each one of the two curves f i (i = 0, 1), we have to determine the tangent plane Πi (u) of the surface F at each point F(u, i) = f i (u). From the above characterization of a geodesic, the tangent plane Πi (u) can be deduced from the osculating plane of the curve f i at each non-inflection point by Πi (u) = f i (u), f i0 (u), bi (u) = f i0 (u) × f i00 (u) . 2. Departure and arrival: we have to choose a starting vector T0 (u) in the tangent plane Π0 (u) at point F(u, 0) and a vector of arrival T1 (u) in the tangent plane Π1 (u) at point F(u, 1). The simplest choice would be to take Ti (u) = ± bi (u), i = 0, 1, according to the orientation of the binormal vectors bi (u). Nevertheless, since that T0 (u) has to be oriented towards point f 1 (u) and that −T1 (u) has to be oriented towards point f 0 (u), another natural choice would
be to take T0 (u) and T1 (u) as the normalized projections D0 (u) and D1 (u) of the vector f 1 (u) − f 0 (u) respectively on the tangent planes Π0 (u) and Π1 (u). More generally, we will consider the following starting and ending vectors, which allow us to get a more general family of interpolating surfaces, Ti (u) = cos[αi (u)]bi (u) + sin[αi (u)]τi (u),
i = 0, 1,
f i0 (u) k f i0 (u)k
with the unit tangent vectors τi (u) = and where the angles α0 and α1 are real functions of the parameter u. They measure the orientation of the vector Ti from vector bi and their influence will be analyzed in Section 3.3, criterion 4. 3. Construction of the surface F(u, v) by Hermite interpolation: F(u, v) = [H0 (v), H1 (v), H2 (v), H3 (v)] f 0 (u) λ0 (u)d(u)T0 (u) × λ1 (u)d(u)T1 (u) f 1 (u) with d(u) = k f 1 (u) − f 0 (u)k and H0 (v) = 2v 3 − 3v 2 + 1, H1 (v) = v 3 −2v 2 +v, H2 (v) = v 3 −v 2 , H3 (v) = 1− H0 (v), and where λ0 and λ1 are real functions of the parameter u and will be discussed in the following. 4. Smoothing: in this last step we have to find optimal parameter functions α0 (u), α1 (u), λ0 (u) and λ1 (u) which provide smooth interpolating surfaces. For this we will use local or global curvature constraints.
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
483
Fig. 2. Influence of parameters λ0 and λ1 .
Fig. 3. Two spatial curves f 0 (u) and f 1 (u) with non-zero torsion.
3.2. First examples Some examples with intuitive choices for T0 (u), T1 (u), λ0 (u), λ1 (u) are proposed here. In all these examples the parameter functions λ0 (u) and λ1 (u) are chosen as two constants λ0 and λ1 , whereas the angles α0 (u) and α1 (u) are identically zero. In Fig. 2, the same initial curves f 0 (u) and f 1 (u) are considered with different values of parameters λ0 and λ1 . These initial curves are in parallel planes and are quite similar except for size. Left, is represented the interpolating surface together with the starting and final vectors for Hermite interpolation: Ti (u) = Di (u), i = 0, 1. Middle and right, the influence of the real parameters λi on the “speed” of the interpolating curves can be observed, either at departure [for λ0 ] or at the arrival [for λ1 ]. Then, we consider in Fig. 3, two initial spatial curves with non-zero torsion. Left, the choice of parameters (λ0 , λ1 ) = (1, 1), corresponds to the classical conditions for Hermite interpolation as the function d(u), see step 3 of the algorithm in Section 3.1, already integrates the distance between the initial and final points. Right, the interpolating curves are stretched by shortening the initial and final tangent vectors. This interpolating method can lead to unwanted situations; essentially due to the inflection points of the initial curves and to the difficulty of finding a coherent orientation for the starting and ending vectors all along the curves. In Fig. 4, angles α0 (u) and α1 (u) are identically zero, and so the Hermite interpolating curves start with the binormal
Fig. 4. Inflection point on the initial curve.
vectors. It can be observed that, even if the osculating plane at the inflection point can be defined by continuity, the binormal vector does not really provide a good starting vector. According to the given initial curves, there is not always a good solution. Anyway, an approach will be proposed in Sections 3.3.1 and 3.3.2. 3.3. Analysis of the algorithm 3.3.1. Inflection points As said above, the difficulty of the first step is to deal with inflection points. As said previously at the end of Section 2, we will assume that in the context of geodesic interpolation, the osculating plane can be defined by continuity at each inflection point. Even in this case, the binormal vector usually changes its
484
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
Fig. 5. Left: Curvature vectors k(u) nE(u) and binormal vectors in bold. Right: Interpolating surface with adequate orientation for the starting vectors.
Fig. 6. Left: Binormal vectors. Right: Interpolating surface.
direction when passing through an inflection point (see Figs. 4 and 5(left)). In Fig. 5(right), we get a suitable interpolating surface, starting from the left curve of Fig. 5, by insuring continuity of the starting vectors, as proposed in Section 3.3.2. See also Fig. 6. From Section 2 (see also [3]) we deduce the following remarks which give information about the local behaviour of the interpolating surface in the neighbourhood of the initial geodesic curves. 1. From Darboux relations we deduce that in an inflection point M = f (u) of a curve running on a surface, the normal curvature kn and the geodesic curvature k g are zero. Thus this curve is locally geodesic in the neighbourhood of M, independent of the surface. Which shows that there is no uniqueness for our reconstruction problem close to an inflection point. 2. A curve running on an elliptic surface can not have inflection points. 3. On each point of an hyperbolic surface, there are exactly two directions where the normal curvature is zero on which a curve can admit an inflection point. From which we deduce. If a geodesic curve f running on a surface F has an inflection point in M = f (u), then M is either an hyperbolic point or
a parabolic point of the surface F. Furthermore, the normal curvature of the surface is zero in the direction f 0 (u). In practice, for numerical applications, we will apply a smoothing procedure to the binormal vector to avoid numerical oscillations near inflection points. See Figs. 12 and 13. 3.3.2. Starting and arrival vectors Even when no inflection occur, we have to be careful in the choice of the directions of the starting vector T0 (u) at point f 0 (u) and the arrival vector T1 (u) at point F1 (u). We have to take account of the rotation of the tangent plane Πi (u) about the tangent line ∆i (u) at point f i (u) along each curve f i . Precisely, the starting vector T0 (u) at point f 0 (u) is in one of the two half-planes of Π0 (u) separated by the line ∆0 (u) and the same is true for the vector of arrival T1 (u). Thus we get an intricate interpolating surface if the starting vector or the vector of arrival goes from one half-plane to the other. Notice that this situation also depends on the relative position of the two curves. Consider for example the two initial curves proposed in Fig. 6(left). The binormal vectors rotate around the curve. So that, these vectors are directed towards each other for some parts of the curve, and to the opposite direction for other parts of the curve. Thus, if for these latter vectors, we set Ti (u) = −bi (u) in order to go towards the other curve, we will get a “tearing” near the middle of the interpolating surface, similar to the one of Fig. 4. So, it seems better to insure continuity
485
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
Fig. 7. Interpolating surfaces obtained according to the four previous criteria with the same initial data.
(with parameter u) of the starting and arrival vectors, as in Fig. 6(right). But of course, if possible, these “pathological” situations have to be avoided in the choice of the initial spatial curves. 3.3.3. Optimization of parameters We are concerned here with the choice of the parameter functions α0 (u), α1 (u) and λ0 (u), λ1 (u). As proposed above, our approach consists in minimizing some energy of curvature of the isoparametric curves v 7→ F(u, ¯ v), for each u¯ ∈ [0, 1]. First, using formula (1), we deduce the curvature and the torsion for a standard parameterization t 7→ ϕ(t) of a 3D curve: k(t) =
kϕ 0 × ϕ 00 k (t) kϕ 0 k3
and
θ (t) =
det(ϕ 0 , ϕ 00 , ϕ 000 ) (t). (5) kϕ 0 × ϕ 00 k2
In this section, we will denote by k(u, ¯ v) and θ (u, ¯ v) the curvature and the torsion of the isoparametric curve F(u, ¯ v), v ∈ [0, 1]. • Criteria We have tested many criteria and of course the final interpolated surface depends strongly on the initial geodesic curves. Nevertheless, from our experiments we consider that the following criteria can be applied. Three of them are local criteria for the interpolated Hermite curves, i.e., for each value u¯ ∈ [0, 1], which allow us to define the parameter functions λ0 (u), λ1 (u) where α0 (u) ≡ 0, α1 (u) ≡ 0. The last one is a global criterion. Criterion 1: This criterion induces a constant speed parameterization. It provides good results for suitable initial conditions.
2 Z 1
d
d F(u,
dv. Min ¯ v)
λ0 ,λ1 0 dv dv Criterion 2: This criterion induces a constant curvature along interpolating curves. Results are reasonably good in
general. Nevertheless we found some divergent cases. 2 Z 1 d k(u, ¯ v) dv. Min λ0 ,λ1 0 dv Criterion 3: This criterion induces planar interpolating curves. Outside inflection points, results are suitable if initial data prohibit getting planar curves. Z 1 Min θ (u, ¯ v)2 dv. λ0 ,λ1 0
We suppose now that the parameters are affine or constant functions, i.e., α0 (u) = (1 − u)α00 + uα01 ,
λ0 (u) ≡ 1,
(1 − u)α10
λ1 (u) ≡ 1.
α1 (u) =
+ uα11 ,
Criterion 4: This global criterion induces planar interpolating curves. Results are suitable in general. ! Z 1 Z 1 2 θ (u, v) dv du. Min α00 ,α01 ,α10 ,α11 0
0
• Results The interpolating surfaces obtained according to the four previous criteria, with the same initial data, are presented in Fig. 7. These initial curves are segments of the two geodesic curves (2) and (3) of the Fig. 12a, and so are only known in a numerical way. Thus, as proposed above, tangent and binormal vectors are numerically evaluated and then smoothed. The last global criterion produces the best interpolating surface for these data. Notice that using the three local criteria, interpolating curves are optimized independently, which can generate oscillations, especially here in Fig. 7 with criterion 3. Furthermore, even though torsion is minimized, the approaches in the last two criteria are different. In the third criterion, we compute the
486
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
Fig. 8a. Geodesic curve on a revolution surface.
λi with fixed αi (locally) and then in the global criterion we compute the αi with fixed λi . 4. Numerical applications 4.1. Analytical geodesics In this section our purpose is to test the reconstruction method with simple surfaces on which we are able to compute exact geodesics. So, we consider surfaces of revolution on which geodesics can be determined using Clairaut’s relation (see [3]). But, notice that even with such elementary surfaces, the differential equations of geodesics are not so easy to integrate. As a first example we deal with a circular cone and then with the sphere. For a regular surface of revolution with the parameterization x = f (u) cos(v), y = f (u) sin(v), z = g(u), f > 0, a geodesic t 7→ G(t) is either a meridian curve v = C and u = u(t), or an “extremal” latitude circle v = v(t) and u = C (i.e., a latitude circle such that the tangent line of the orthogonal generating curve is parallel to the axis of revolution), or a curve which satisfies Clairaut’s relation: r cos(θ ) = C, where r is the distance to the axis of revolution and θ is the angle (h(t), G 0 (t)) of the geodesic with the latitude circle that intersects it (see Fig. 8a). 4.1.1. Geodesics on a cone We consider here a circular cone generated by the rotation of a straight line around the z-axis. If a is the angle of the line with the z-axis, we can write f (u) = u and g(u) = A u, with A = tan−1 (a) and the cone is defined by x = u cos v, y = u sin v, z = A u for u ≥ 0. Using the previous results, geodesic curves of this cone are parameterized, according to the angle v, by cos v u0 sin v , G(v) = cos ((v − v0 ) sin a) A where u 0 and v0 are the values of parameters such that the angle (h(v), G 0 (v)) is zero. That is, such that the geodesic curve is tangent to a latitude circle of the cone (see Fig. 8b). Using the previous parameterization, we define two geodesic curves on a cone in the following way. Either the two geodesic curves are aligned (their starting and ending points are on the
Fig. 8b. Geodesic on a circular cone.
same meridian line) as in Fig. 9, or they are shifted (they are not aligned) as in Fig. 10. For the reconstruction, in Figs. 9 and 10(left), we use the vectors Di (u) directed towards the other curve, as defined in step 2 of the algorithm in Section 3.1, whereas in Figs. 9 and 10(right), we use the binormal vectors. 4.1.2. Geodesics on the sphere Two shifted geodesic curves are now considered on a sphere, see Fig. 11. For the reconstruction, we use either the binormal vectors (left), or the vectors Di (u) (right). 4.2. Numerical geodesics We consider now parameterized surfaces on which geodesic curves are evaluated by employing an algorithm developed in [8]. So that we get only numerical data concerning these geodesics. Precisely, as shown in Fig. 12a, we evaluate four geodesic curves on a given Nurbs surface, using the previous cited algorithm. These geodesics are numbered from (1) to (4). Since only numerical data are available concerning these geodesic curves, the surface reconstruction is then carried out by means of a numerical evaluation of the first and the second derivatives. Then, from these first and second evaluated derivatives, the tangent and binormal vectors are computed, implying important oscillations near inflection points (see Fig. 13a). The three interpolating surfaces between geodesics (1) and (2), (2) and (3), (3) and (4) are reconstructed using these evaluated
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
487
Fig. 9. Reconstruction with two “aligned” geodesic curves.
Fig. 10. Reconstruction with two shifted geodesic curves.
Fig. 11. Two geodesic curves shifted on the sphere.
numerical binormal vectors before smoothing (Fig. 12b, left) and after smoothing (Fig. 12b, right). For each one of these four geodesic curves, the variations of the three cartesian components, labelled by x, y, z, of the binormal vectors are shown before smoothing (Fig. 13a) and after smoothing by means of least square method and cubic Bsplines (Fig. 13b).
5. Surface reconstruction via geodesic interpolation We will now experiment with this method on real surfaces. For this purpose, we need to get geodesics from the surface we have to reconstruct. Thus, we use Morphosense, which is a ribbon with embedded microsensors. These sensors are microaccelerometers (which give a static
488
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
measure that is the acceleration of the gravity, and then give orientation information with respect to the vertical) and micromagnetometers (which can provide angular information based on the earth’s magnetic field). The combination of this information gives the exact orientation of the sensor. So if we equip a ribbon with these two kinds of measuring devices, we will get the directions of the tangency on the sensors’ position. Let us note that we do not know the positions of the sensors nor their tangents, but only the direction of the tangents. In order to reconstruct spatial curves, we have equipped Morphosense with 16 3D-accelerometers alternating with 16 2D-magnetometers nearly every three centimeters. Morphosense is then able to provide its own shape (i.e., the spatial curve it represents). In order to get some geodesic curves from a surface, we will put Morphosense on this surface with a G 1 geometric contact (see Theoretical analysis, below, to see that the ribbon describes a geodesic of the surface). So the experimentation is as follows: we put the ribbon on the surface in some positions nearly parallel, and acquire the tangential informations. Then, the method for the surface reconstruction proceeds in two steps: first we have to reconstruct the spatial curves under tangential and length constraints (see below for a short overview), then we use the geodesic interpolation presented in this paper. 5.1. Reconstruction of spatial curves We resume here the process of curve reconstruction (see [9, 10] for a complete description). The data given by Morphosense are the separation of the sensors along the curve (the ribbon) and orientation of these sensors in space. Since we have length constraints (sensors separation), we work with the arc-length parameter for the reconstruction of the spatial curves. The derivative is thus a unit vector representing a curve lying on the unit sphere and the data are samples of this derivative at known parameters. Then, in order to interpolate this data by a curve staying on the sphere, we use cubic splines on the unit sphere, as proposed in [6].
Fig. 12a. Initial surface with 4 evaluated geodesics numbered from 1 to 4.
Thus, the method proceeds in two steps: first, we reconstruct the derivative function, and then we integrate this derivative to obtain the desired 3D curve. A numerical method is used for this second part. An example of such a spatial reconstruction is presented in Fig. 15. Left, is shown the original curve (thin) and the reconstructed one (bold). Right, is the derivative reconstructed on the sphere with values of sensors (dots). 5.2. Theoretical analysis A developable surface D is isometric to a plane, so that any geodesic curve γ on D can be associated with a straight line on a plane. Thus, considering any regular surface F, if we assume that the two surfaces D and F have a G 1 contact along a common curve γ and that γ is a geodesic curve on D, we will show that γ is a geodesic curve on the surface F. In practice, to ensure this G 1 contact we consider a ribbon D on the surface D, of width 2, centered on the curve γ . One can understand that even with a thin ribbon the contact is theoretical on hyperbolic points of the surface F. Then, using the fact that geodesic curves are analogous to straight lines in the plane, we will give a mechanical characterization of geodesics related with the G 1 installation of the ribbon on F. Precisely, Let (τ (s), V (s), N (s)) the Darboux frame at point γ (s) on the developable surface D . When installing the ribbon on the surface F, this one is deformed in an isometric way which only authorize local rotations around the vectors τ (s) and V (s) along the curve γ . The idea is that geodesics have to be straight, so that local rotations around the vector N (s) are forbidden. The curve γ on the ribbon D will be said N -rigid. Definition 1. A regular curve γ (s) running on a surface F is said to be N -rigid if its Darboux frame can not locally rotate (with parameter s) around the vector N (s). Proposition 1. Let γ (s) be a N -rigid regular curve running on a surface F. Then γ is a geodesic curve of the surface F. Proof. Since the curve γ is N -rigid, the Darboux frame (τ (s)+ dτ, V (s) + dV, N (s) + dN ) at point γ (s + ds) can only be deduced from the Darboux frame at point γ (s) by rotations around vectors τ (s) and V (s). Thus, considering elementary rotations around vector τ (s) by an angle dωx and around vector V (s) by an angle dω y , we can write τ (s) + dτ cos dω y 0 − sin dω y V (s) + dV = 0 1 0 N (s) + dN sin dω y 0 cos dω y
Fig. 12b. Reconstructed surfaces without (left) and with (right) smoothing of the binormal vectors.
489
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
Fig. 13a. Binormal vectors along the four evaluated geodesics.
Fig. 13b. Binormal vectors after smoothing.
1 0 × 0 cos dωx 0 sin dωx
0 τ (s) − sin dωx V (s) , cos dωx N (s)
which, assuming cos dωx ' 1 ' cos dω y , sin dωx ' dωx , sin dω y ' dω y , leads to dτ 0 0 −dω y τ (s) dV = 0 0 −dωx V (s) . (6) dN dω y dωx 0 N (s) Thus, comparing this result with Darboux relations (3), we get k g (s) ≡ 0 which proves that γ is a geodesic curve. Notice that we get exactly the same result if we first apply an elementary rotation around V (s). Proposition 2. Let γ (s) be a geodesic curve of the surface F. Then γ is N -rigid. Proof. The Darboux frame (τ (s)+dτ, V (s)+dV, N (s)+dN ) at point γ (s + ds) can be deduced from the Darboux frame at point γ (s) by 3 elementary rotations Rdωx , Rdω y , Rdωz
around vectors τ (s), V (s), N (s) by angles dωx , dω y , dωz . Thus, assuming cos dωx ' 1, cos dω y ' 1, cos dωz ' 1, sin dωx ' dωx , sin dω y ' dω y , sin dωz ' dωz , we have τ (s) + dτ τ (s) V (s) + dV = Rdωz ◦ Rdω y ◦ Rdωx V (s) , (7) N (s) N (s) + dN which leads to dτ 0 dV = −dωz dN dω y
−dωz 0 dωx
−dω y τ (s) −dωx V (s) . 0 N (s)
(8)
Now suppose that γ is a geodesic curve, that is k g (s) ≡ 0. Thus, from Darboux relations (3) we deduce that dωz = 0, which proves that the curve γ is N -rigid. Notice that the result does not depend on which order we apply the 3 elementary rotations. Proposition 3. Consider two regular oriented surfaces F1 and F2 with a G 1 contact along a common curve γ . Then, if γ is
490
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
Fig. 14. “Naked” Ribbon Morphosense with sensors and below, a real time demonstration.
the ribbon and F along the common curve γ ensures that the Darboux frames of γ on each of these two surfaces agree (regardless of the orientation). So that the curve γ is also N rigid on F. Then, from Proposition 1, we deduce that γ is a geodesic curve on F. Notice that Proposition 3 provides a more direct way to prove this result. Nevertheless, Propositions 1 and 2 connect the mechanical properties of the ribbon to its geometric behaviour. 5.3. Numerical results Fig. 15. Spatial curve reconstruction.
a geodesic curve on one surface, it is a geodesic curve on the other surface. Proof. Let us denote by (τi (s), Vi (s), Ni (s)), i = 1, 2, the Darboux frames of the curve γ on each one of the two surfaces Fi . Of course we have τ1 (s) = τ2 (s) for each value of the parameter s and since these two surfaces have a G 1 contact along the curve γ , we have N1 ≡ ±N2 . If necessary we can modify the orientation of one the two surfaces to get N1 ≡ N2 , so that the two Darboux frames are the same all along the curve γ . Finally, the result can be deduced from Darboux relations (3). From these results we are now able to verify that our ribbon actually follows geodesic curves on the surface F. The ribbon is inextensible and not “flexible”. It has the same behaviour as a thin strip of paper, i.e., as a developable surface. Consider a geodesic curve γ on the ribbon (this geodesic is isometric to a straight line of the plane). From Proposition 2, we deduce that this curve γ is N -rigid on the ribbon. Then, using the same arguments as in Proposition 3, the G 1 contact between
We first consider a numerical model (see Fig. 16, top) of a real surface (in Fig. 17). Five geodesic curves, numbered from 1 to 5, are evaluated and then the reconstruction is achieved via geodesic interpolation. In Fig. 16, each surface Si j is the interpolating surface between geodesics i and j and S12345 is the whole interpolating surface. Figure 17(left), shows the installation of the ribbon on the real surface, and (right) is the reconstructed surface via geodesic interpolation. 6. Conclusion In this paper it has been shown that numerical and real surfaces can be reconstructed from some geodesic running on this surface. These geodesics can be captured using a ribbon of micro-sensors so that the reconstruction process is really intrinsic. Even if this point is not at the heart of this paper, our experiments show data given by micro-sensors to be very sensitive. Thus, our process requires smoothing of this data. Nevertheless, the reconstruction is simple and fast and so could be suitable for deformation of surfaces.
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
Fig. 16. Reconstruction via geodesic interpolation.
Fig. 17. Reconstruction from ribbon’s data representing transversal geodesic curves.
491
492
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
Acknowledgments We are very grateful for the anonymous reviewers for their detailed constructive comments. References [1] Coons S. Surface patches and B-spline curves. In: Barnhill R, Riesenfeld R, editors. Computer aided geometric design. Academic Press; 1974. [2] Coons S. Surfaces for computer aided design, Technical report, M.I.T. 1964. Available as AD 663 504 from the National Technical Information Service. Springfield, VA 22161. [3] Do Carmo MP. Differential geometry of curves and surfaces. PrenticeHall; 1976. [4] Farin G. Curves and surfaces for CAGD, 5th ed. Academic Press; 2002. [5] Gordon W. An operator calculus for surface and volume modeling. IEEE
Computer Graphics and Applications 1983;3:18–22. [6] Nielson GM. ν-quaternion splines for the smooth interpolation of orientations. IEEE Transactions on Visualization and Computer Graphics 2004;10(2):224–9. [7] Pham-Trong V. (2001) D´etermination g´eom´etrique de chemins g´eod´esiques sur des surfaces de subdivision. Th`ese de l’universit´e Joseph Fourier, Math´ematiques Appliqu´ees. Laboratoire LMC-IMAG. 28 Septembre 2001. [8] Pham-Trong V, Biard L, Szafran N. Pseudo-geodesics on threedimensional surfaces and pseudo-geodesic meshes. Numerical Algorithms 2001;26:305–15. [9] Sprynski N, Lacolle B, Biard L, David D. In: Chenin P, Lyche T, Schumaker LL, editors. Curve and surface reconstruction via tangential informations, curve and surface design, Avignon 2006. Nashboro Press; 2007. p. 254–63. [10] Sprynski N. Reconstruction de courbes et surfaces a` partir de donn´ees tangentielles, Laboratoires CEA/LETI, LJK. Th`ese de l’universit´e Joseph Fourier, Grenoble. 5 Juillet 2007.
Surface reconstruction via geodesic interpolation N. Sprynski b , N. Szafran a , B. Lacolle a , L. Biard a,∗ a Laboratoire Jean Kuntzmann, Universit´e Joseph Fourier - Grenoble, France b CEA MINATEC, LETI/DCIS/SMOC, Grenoble, France
Received 24 May 2007; accepted 10 January 2008
Abstract This paper is concerned with reconstruction of numerical or real surfaces based on the knowledge of some geodesic curves on the surface. So, considering two regular 3D-curves f 0 (t) and f 1 (t), our purpose is to construct a surface which interpolates these two curves in such a way that these two curves are geodesics on this surface. This will be accomplished using Hermite interpolation. For a real surface, it will be shown that geodesics can be acquired using a ribbon of micro-sensors. c 2008 Elsevier Ltd. All rights reserved.
Keywords: Geodesic curves; Surface reconstruction; Hermite interpolation; Curvature minimization; Motion capture; Embedded micro-sensors
1. Introduction This paper is concerned with reconstruction of surfaces from curves running on this surface. Thus, given arbitrary curves on a surface, Coons methods [1,2,5,4] allow us to get a surface which interpolates the surface. Notice that Coons processes are mostly used to create/define 3D objects from some spatial curves. These 3D curves are designed/computed first. They provide the main characteristic features of the desired surface and model the essential shapes of the object. The Coons process, is then performed as a “filling” method to get the final surface. So, the surface is essentially modeled by the initial spatial curves. Thus, concerning a real 3D surface, the immediate question is how to get lines from this surface. Notice that we are concerned by the reconstruction problem of a real surface for which no parameterization is accessible. Then, in order to get a consistent reconstruction process, these lines must represent intrinsic characteristic features of the surface. Fortunately we are able to provide a positive answer to these two points. A natural way would be to acquire geodesic curves running on the surface. That can be achieved employing a ribbon of embedded micro-sensors (see Figs. 14 and 17) and it will be shown, in Section 5 that such a ribbon actually follows geodesic ∗ Corresponding author.
E-mail address: [email protected] (L. Biard). c 2008 Elsevier Ltd. All rights reserved. 0010-4485/$ - see front matter doi:10.1016/j.cad.2008.01.005
curves on the surface. This ribbon is inextensible and not “flexible” and has the same behaviour as a thin strip of paper, i.e., as a developable surface. Its mechanical properties are modeled in Section 5. The micro-sensors are able to give their own orientation in space from which are deduced data about the local tangency of the ribbon. This technology has been developed at the CEA-LETI, a research center of technology specialized in micro-electronics. Thus, from the tangential information of the ribbon, we can deduce its 3D shape. This reconstruction process is described in [9] and works for planar and spatial curves. Thus, we are now concerned with reconstruction of surfaces based on the knowledge of some geodesic curves on the surface. So, considering two regular 3D-curves f 0 (t) and f 1 (t), our purpose is to construct a surface which interpolates these two curves in such a way that these two curves are geodesics on this surface. This problem of surface reconstruction using only tangential orientations without any position data seems to have not yet been handled and the need to solve such a problem has become important in the field of motion capture. Notice that this is not an envelope problem in which the tangent lines are fixed. Here we only know the directions of these tangent lines and, moreover, we have the added constraint that the length between two consecutive contact points are given. This paper is organized as follows: Section 2 deals with fundamental properties of curves on surfaces. A solution to the geodesic interpolation problem is then presented in Section 3
481
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
with examples and an analysis of the algorithm. In Section 4 we consider numerical applications. Precisely, we consider numerical models of surfaces on which exact geodesic curves are evaluated with an algorithm issued from [7,8]. Then, the surface is reconstructed from these geodesics. Finally, in Section 5 we deal with reconstruction of real surfaces using the ribbon of micro-sensors. In particular, we show that such a ribbon follows geodesic curves on the surface. 2. Background on curves and surfaces In the following, all curves and surfaces are considered to be sufficiently smooth (at least of class C 3 ). The parameters of functions will sometimes be omitted when no confusion can occur and the dot will denote the inner product. For a regular 3D curve ϕ : t 7→ ϕ(t) we will denote by ϕ : s 7→ ϕ(s) its normal parameterization by arc length with the following properties kϕ 0 (s)k = 1 and ϕ 0 (s) ⊥ ϕ 00 (s) for each value of s. • With each point ϕ(s) of such a curve such that ϕ 00 (s) 6= 0, we associate the Serret–Frenet frame (τ (s), n(s), b(s)) where τ (s) = ϕ 0 (s), n(s) = ϕ 00 (s)/kϕ 00 (t)k and b(s) = τ (s)×n(s) are respectively the unit tangent vector, the unit normal vector and the binormal vector of the curve at point ϕ(s). The unit normal vector n(s) is also called the main normal. Considering the arc length, the Serret–Frenet frame satisfies the Serret–Frenet relations: τ (s) 0 k(s) 0 τ (s) d n(s) = −k(s) 0 θ (s) n(s) , ds b(s) 0 −θ (s) 0 b(s) where k(s) and θ (s) are the curvature and the torsion of the curve at point ϕ(s) defined by: k(s) = kϕ 00 (s)k
θ (s) =
and
det(ϕ 0 (s), ϕ 00 (s), ϕ 000 (s)) . (1) kϕ 00 (s)k2
The osculating plane of the curve at point ϕ(s) is defined by the two vectors τ (s) and n(s) and does not depend on the parameterization of the curve. If for some parameter s, the curvature k(s) = 0, then ϕ 00 (s) = 0, so that the normal vector n(s) and the osculating plane can not be defined. Such a point is called an inflection point of the curve. • For a regular oriented surface (u, v) 7→ F(u, v), we consider the normal vector N (u, v) and the Gaussian curvature K T (u, v) at point F(u, v) defined by N (u, v) =
Fu0 × Fv0 , kFu0 × Fv0 k
K T (u, v) = K 1 K 2 =
and
L = N · Fu002 ,
00 M = N · Fuv ,
g
which define the normal curvature kn (s), the geodesic curvature k g (s) and the geodesic torsion θg (s) by kn (s) =
dτ · N, ds
k g (s) =
dτ · V, ds
θg (s) =
dN · V. ds
(4)
• A regular curve f (s), running on the surface F is a geodesic curve of this surface if its geodesic curvature k g is identically zero. Thus, from k g ≡ 0 we deduce with the Darboux relations that dτ ds is collinear with the vector N for each value of the parameter, so that the Frenet and Darboux frames agree regardless of the orientation. Thus, since the osculating plane of a curve is unchanged under reparameterization, we get the following useful characterization of geodesic curves. A regular curve t 7→ f (t) running on a surface F(u, v) is a geodesic of F if its osculating plane at each non-inflection point f (t) is normal to the tangent plane of the surface F at point F(u(t), v(t)) = f (t). In the case of an isolated inflection point of the curve, we may encounter situations, see e.g. [3], where the osculating plane can not be defined by continuity as the limit planes approaching from the left and from the right side will disagree. However, in the given situation considering geodesic curves, it turn out later that the geodesic’s Frenet frame and the Darboux frame will agree for isolated inflection points, because, for geodesics, the main normal is collinear with respective surface normal. Therefore the continuity of the Darboux frame implies the continuity of the osculating plane in the corresponding Frenet frame. In the case of non-isolated inflection points, the vanishing of the curvature on a domain implies that the associated points belong to a straight line (so, a geodesic curve) running on the surface. Their Darboux frame can be defined on each point whereas the Frenet frame can not be defined at any point.
(2)
where K 1 and K 2 are the principal curvatures and where F = Fu0 · Fv0 ,
n
3. Geodesic interpolation
L N − M2 , EG − F2
E = Fu0 · Fu0 ,
elliptic points (K T > 0), hyperbolic points (K T < 0) and parabolic points (K T = 0). • Consider a curve f (s) = F(u(s), v(s)) running on the surface F. With each point f (s) we associate the Darboux frame (τ (s), V (s), N (s)) where τ (s) is the unit tangent vector of the curve f , N (s) is the unit normal vector of the surface F at point F(u(s), v(s)) = f (s) and V (s) = N (s) × τ (s). Considering the arc length, the Darboux frame satisfies the following Darboux relations: τ (s) 0 k g (s) kn (s) τ (s) d V (s) = −k g (s) 0 −θg (s) V (s) (3) ds N (s) −k (s) θ (s) 0 N (s)
G = Fv0 · Fv0 ; N = N · Fv002 ,
are the coefficients of the first and second fundamental form of the surface. The Gaussian curvature allows us to distinguish
Consider now two regular 3D curves f 0 (u) and f 1 (u), u ∈ [0, 1], such that f 0 (u) 6= f 1 (u) for each parameter u, and remember that our goal is to construct a surface which interpolates these two curves in such a way that these curves are geodesics on this surface. Precisely, let us denote by F(u, v), with (u, v) ∈ [0, 1]2 the interpolating surface. The surface F has to fulfill the two following properties.
482
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
Fig. 1. Geodesic interpolation of the two 3D-curves f 0 (u) and f 1 (u).
• Interpolating property: F(u, 0) = f 0 (u) and F(u, 1) = f 1 (u) for each u ∈ [0, 1]. • Geodesic property: for i ∈ {0, 1}, the osculating plane of the curve f i at each point f i (u) is normal to the tangent plane of the surface F at point F(u, i). Notice that by the following construction, each one of the two (geodesic) curves f i (u) will be an isoparametric curve of the surface F. 3.1. Main steps of the algorithm The algorithm can be described by the following steps (see also Fig. 1). 1. Tangent plane: for each one of the two curves f i (i = 0, 1), we have to determine the tangent plane Πi (u) of the surface F at each point F(u, i) = f i (u). From the above characterization of a geodesic, the tangent plane Πi (u) can be deduced from the osculating plane of the curve f i at each non-inflection point by Πi (u) = f i (u), f i0 (u), bi (u) = f i0 (u) × f i00 (u) . 2. Departure and arrival: we have to choose a starting vector T0 (u) in the tangent plane Π0 (u) at point F(u, 0) and a vector of arrival T1 (u) in the tangent plane Π1 (u) at point F(u, 1). The simplest choice would be to take Ti (u) = ± bi (u), i = 0, 1, according to the orientation of the binormal vectors bi (u). Nevertheless, since that T0 (u) has to be oriented towards point f 1 (u) and that −T1 (u) has to be oriented towards point f 0 (u), another natural choice would
be to take T0 (u) and T1 (u) as the normalized projections D0 (u) and D1 (u) of the vector f 1 (u) − f 0 (u) respectively on the tangent planes Π0 (u) and Π1 (u). More generally, we will consider the following starting and ending vectors, which allow us to get a more general family of interpolating surfaces, Ti (u) = cos[αi (u)]bi (u) + sin[αi (u)]τi (u),
i = 0, 1,
f i0 (u) k f i0 (u)k
with the unit tangent vectors τi (u) = and where the angles α0 and α1 are real functions of the parameter u. They measure the orientation of the vector Ti from vector bi and their influence will be analyzed in Section 3.3, criterion 4. 3. Construction of the surface F(u, v) by Hermite interpolation: F(u, v) = [H0 (v), H1 (v), H2 (v), H3 (v)] f 0 (u) λ0 (u)d(u)T0 (u) × λ1 (u)d(u)T1 (u) f 1 (u) with d(u) = k f 1 (u) − f 0 (u)k and H0 (v) = 2v 3 − 3v 2 + 1, H1 (v) = v 3 −2v 2 +v, H2 (v) = v 3 −v 2 , H3 (v) = 1− H0 (v), and where λ0 and λ1 are real functions of the parameter u and will be discussed in the following. 4. Smoothing: in this last step we have to find optimal parameter functions α0 (u), α1 (u), λ0 (u) and λ1 (u) which provide smooth interpolating surfaces. For this we will use local or global curvature constraints.
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
483
Fig. 2. Influence of parameters λ0 and λ1 .
Fig. 3. Two spatial curves f 0 (u) and f 1 (u) with non-zero torsion.
3.2. First examples Some examples with intuitive choices for T0 (u), T1 (u), λ0 (u), λ1 (u) are proposed here. In all these examples the parameter functions λ0 (u) and λ1 (u) are chosen as two constants λ0 and λ1 , whereas the angles α0 (u) and α1 (u) are identically zero. In Fig. 2, the same initial curves f 0 (u) and f 1 (u) are considered with different values of parameters λ0 and λ1 . These initial curves are in parallel planes and are quite similar except for size. Left, is represented the interpolating surface together with the starting and final vectors for Hermite interpolation: Ti (u) = Di (u), i = 0, 1. Middle and right, the influence of the real parameters λi on the “speed” of the interpolating curves can be observed, either at departure [for λ0 ] or at the arrival [for λ1 ]. Then, we consider in Fig. 3, two initial spatial curves with non-zero torsion. Left, the choice of parameters (λ0 , λ1 ) = (1, 1), corresponds to the classical conditions for Hermite interpolation as the function d(u), see step 3 of the algorithm in Section 3.1, already integrates the distance between the initial and final points. Right, the interpolating curves are stretched by shortening the initial and final tangent vectors. This interpolating method can lead to unwanted situations; essentially due to the inflection points of the initial curves and to the difficulty of finding a coherent orientation for the starting and ending vectors all along the curves. In Fig. 4, angles α0 (u) and α1 (u) are identically zero, and so the Hermite interpolating curves start with the binormal
Fig. 4. Inflection point on the initial curve.
vectors. It can be observed that, even if the osculating plane at the inflection point can be defined by continuity, the binormal vector does not really provide a good starting vector. According to the given initial curves, there is not always a good solution. Anyway, an approach will be proposed in Sections 3.3.1 and 3.3.2. 3.3. Analysis of the algorithm 3.3.1. Inflection points As said above, the difficulty of the first step is to deal with inflection points. As said previously at the end of Section 2, we will assume that in the context of geodesic interpolation, the osculating plane can be defined by continuity at each inflection point. Even in this case, the binormal vector usually changes its
484
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
Fig. 5. Left: Curvature vectors k(u) nE(u) and binormal vectors in bold. Right: Interpolating surface with adequate orientation for the starting vectors.
Fig. 6. Left: Binormal vectors. Right: Interpolating surface.
direction when passing through an inflection point (see Figs. 4 and 5(left)). In Fig. 5(right), we get a suitable interpolating surface, starting from the left curve of Fig. 5, by insuring continuity of the starting vectors, as proposed in Section 3.3.2. See also Fig. 6. From Section 2 (see also [3]) we deduce the following remarks which give information about the local behaviour of the interpolating surface in the neighbourhood of the initial geodesic curves. 1. From Darboux relations we deduce that in an inflection point M = f (u) of a curve running on a surface, the normal curvature kn and the geodesic curvature k g are zero. Thus this curve is locally geodesic in the neighbourhood of M, independent of the surface. Which shows that there is no uniqueness for our reconstruction problem close to an inflection point. 2. A curve running on an elliptic surface can not have inflection points. 3. On each point of an hyperbolic surface, there are exactly two directions where the normal curvature is zero on which a curve can admit an inflection point. From which we deduce. If a geodesic curve f running on a surface F has an inflection point in M = f (u), then M is either an hyperbolic point or
a parabolic point of the surface F. Furthermore, the normal curvature of the surface is zero in the direction f 0 (u). In practice, for numerical applications, we will apply a smoothing procedure to the binormal vector to avoid numerical oscillations near inflection points. See Figs. 12 and 13. 3.3.2. Starting and arrival vectors Even when no inflection occur, we have to be careful in the choice of the directions of the starting vector T0 (u) at point f 0 (u) and the arrival vector T1 (u) at point F1 (u). We have to take account of the rotation of the tangent plane Πi (u) about the tangent line ∆i (u) at point f i (u) along each curve f i . Precisely, the starting vector T0 (u) at point f 0 (u) is in one of the two half-planes of Π0 (u) separated by the line ∆0 (u) and the same is true for the vector of arrival T1 (u). Thus we get an intricate interpolating surface if the starting vector or the vector of arrival goes from one half-plane to the other. Notice that this situation also depends on the relative position of the two curves. Consider for example the two initial curves proposed in Fig. 6(left). The binormal vectors rotate around the curve. So that, these vectors are directed towards each other for some parts of the curve, and to the opposite direction for other parts of the curve. Thus, if for these latter vectors, we set Ti (u) = −bi (u) in order to go towards the other curve, we will get a “tearing” near the middle of the interpolating surface, similar to the one of Fig. 4. So, it seems better to insure continuity
485
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
Fig. 7. Interpolating surfaces obtained according to the four previous criteria with the same initial data.
(with parameter u) of the starting and arrival vectors, as in Fig. 6(right). But of course, if possible, these “pathological” situations have to be avoided in the choice of the initial spatial curves. 3.3.3. Optimization of parameters We are concerned here with the choice of the parameter functions α0 (u), α1 (u) and λ0 (u), λ1 (u). As proposed above, our approach consists in minimizing some energy of curvature of the isoparametric curves v 7→ F(u, ¯ v), for each u¯ ∈ [0, 1]. First, using formula (1), we deduce the curvature and the torsion for a standard parameterization t 7→ ϕ(t) of a 3D curve: k(t) =
kϕ 0 × ϕ 00 k (t) kϕ 0 k3
and
θ (t) =
det(ϕ 0 , ϕ 00 , ϕ 000 ) (t). (5) kϕ 0 × ϕ 00 k2
In this section, we will denote by k(u, ¯ v) and θ (u, ¯ v) the curvature and the torsion of the isoparametric curve F(u, ¯ v), v ∈ [0, 1]. • Criteria We have tested many criteria and of course the final interpolated surface depends strongly on the initial geodesic curves. Nevertheless, from our experiments we consider that the following criteria can be applied. Three of them are local criteria for the interpolated Hermite curves, i.e., for each value u¯ ∈ [0, 1], which allow us to define the parameter functions λ0 (u), λ1 (u) where α0 (u) ≡ 0, α1 (u) ≡ 0. The last one is a global criterion. Criterion 1: This criterion induces a constant speed parameterization. It provides good results for suitable initial conditions.
2 Z 1
d
d F(u,
dv. Min ¯ v)
λ0 ,λ1 0 dv dv Criterion 2: This criterion induces a constant curvature along interpolating curves. Results are reasonably good in
general. Nevertheless we found some divergent cases. 2 Z 1 d k(u, ¯ v) dv. Min λ0 ,λ1 0 dv Criterion 3: This criterion induces planar interpolating curves. Outside inflection points, results are suitable if initial data prohibit getting planar curves. Z 1 Min θ (u, ¯ v)2 dv. λ0 ,λ1 0
We suppose now that the parameters are affine or constant functions, i.e., α0 (u) = (1 − u)α00 + uα01 ,
λ0 (u) ≡ 1,
(1 − u)α10
λ1 (u) ≡ 1.
α1 (u) =
+ uα11 ,
Criterion 4: This global criterion induces planar interpolating curves. Results are suitable in general. ! Z 1 Z 1 2 θ (u, v) dv du. Min α00 ,α01 ,α10 ,α11 0
0
• Results The interpolating surfaces obtained according to the four previous criteria, with the same initial data, are presented in Fig. 7. These initial curves are segments of the two geodesic curves (2) and (3) of the Fig. 12a, and so are only known in a numerical way. Thus, as proposed above, tangent and binormal vectors are numerically evaluated and then smoothed. The last global criterion produces the best interpolating surface for these data. Notice that using the three local criteria, interpolating curves are optimized independently, which can generate oscillations, especially here in Fig. 7 with criterion 3. Furthermore, even though torsion is minimized, the approaches in the last two criteria are different. In the third criterion, we compute the
486
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
Fig. 8a. Geodesic curve on a revolution surface.
λi with fixed αi (locally) and then in the global criterion we compute the αi with fixed λi . 4. Numerical applications 4.1. Analytical geodesics In this section our purpose is to test the reconstruction method with simple surfaces on which we are able to compute exact geodesics. So, we consider surfaces of revolution on which geodesics can be determined using Clairaut’s relation (see [3]). But, notice that even with such elementary surfaces, the differential equations of geodesics are not so easy to integrate. As a first example we deal with a circular cone and then with the sphere. For a regular surface of revolution with the parameterization x = f (u) cos(v), y = f (u) sin(v), z = g(u), f > 0, a geodesic t 7→ G(t) is either a meridian curve v = C and u = u(t), or an “extremal” latitude circle v = v(t) and u = C (i.e., a latitude circle such that the tangent line of the orthogonal generating curve is parallel to the axis of revolution), or a curve which satisfies Clairaut’s relation: r cos(θ ) = C, where r is the distance to the axis of revolution and θ is the angle (h(t), G 0 (t)) of the geodesic with the latitude circle that intersects it (see Fig. 8a). 4.1.1. Geodesics on a cone We consider here a circular cone generated by the rotation of a straight line around the z-axis. If a is the angle of the line with the z-axis, we can write f (u) = u and g(u) = A u, with A = tan−1 (a) and the cone is defined by x = u cos v, y = u sin v, z = A u for u ≥ 0. Using the previous results, geodesic curves of this cone are parameterized, according to the angle v, by cos v u0 sin v , G(v) = cos ((v − v0 ) sin a) A where u 0 and v0 are the values of parameters such that the angle (h(v), G 0 (v)) is zero. That is, such that the geodesic curve is tangent to a latitude circle of the cone (see Fig. 8b). Using the previous parameterization, we define two geodesic curves on a cone in the following way. Either the two geodesic curves are aligned (their starting and ending points are on the
Fig. 8b. Geodesic on a circular cone.
same meridian line) as in Fig. 9, or they are shifted (they are not aligned) as in Fig. 10. For the reconstruction, in Figs. 9 and 10(left), we use the vectors Di (u) directed towards the other curve, as defined in step 2 of the algorithm in Section 3.1, whereas in Figs. 9 and 10(right), we use the binormal vectors. 4.1.2. Geodesics on the sphere Two shifted geodesic curves are now considered on a sphere, see Fig. 11. For the reconstruction, we use either the binormal vectors (left), or the vectors Di (u) (right). 4.2. Numerical geodesics We consider now parameterized surfaces on which geodesic curves are evaluated by employing an algorithm developed in [8]. So that we get only numerical data concerning these geodesics. Precisely, as shown in Fig. 12a, we evaluate four geodesic curves on a given Nurbs surface, using the previous cited algorithm. These geodesics are numbered from (1) to (4). Since only numerical data are available concerning these geodesic curves, the surface reconstruction is then carried out by means of a numerical evaluation of the first and the second derivatives. Then, from these first and second evaluated derivatives, the tangent and binormal vectors are computed, implying important oscillations near inflection points (see Fig. 13a). The three interpolating surfaces between geodesics (1) and (2), (2) and (3), (3) and (4) are reconstructed using these evaluated
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
487
Fig. 9. Reconstruction with two “aligned” geodesic curves.
Fig. 10. Reconstruction with two shifted geodesic curves.
Fig. 11. Two geodesic curves shifted on the sphere.
numerical binormal vectors before smoothing (Fig. 12b, left) and after smoothing (Fig. 12b, right). For each one of these four geodesic curves, the variations of the three cartesian components, labelled by x, y, z, of the binormal vectors are shown before smoothing (Fig. 13a) and after smoothing by means of least square method and cubic Bsplines (Fig. 13b).
5. Surface reconstruction via geodesic interpolation We will now experiment with this method on real surfaces. For this purpose, we need to get geodesics from the surface we have to reconstruct. Thus, we use Morphosense, which is a ribbon with embedded microsensors. These sensors are microaccelerometers (which give a static
488
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
measure that is the acceleration of the gravity, and then give orientation information with respect to the vertical) and micromagnetometers (which can provide angular information based on the earth’s magnetic field). The combination of this information gives the exact orientation of the sensor. So if we equip a ribbon with these two kinds of measuring devices, we will get the directions of the tangency on the sensors’ position. Let us note that we do not know the positions of the sensors nor their tangents, but only the direction of the tangents. In order to reconstruct spatial curves, we have equipped Morphosense with 16 3D-accelerometers alternating with 16 2D-magnetometers nearly every three centimeters. Morphosense is then able to provide its own shape (i.e., the spatial curve it represents). In order to get some geodesic curves from a surface, we will put Morphosense on this surface with a G 1 geometric contact (see Theoretical analysis, below, to see that the ribbon describes a geodesic of the surface). So the experimentation is as follows: we put the ribbon on the surface in some positions nearly parallel, and acquire the tangential informations. Then, the method for the surface reconstruction proceeds in two steps: first we have to reconstruct the spatial curves under tangential and length constraints (see below for a short overview), then we use the geodesic interpolation presented in this paper. 5.1. Reconstruction of spatial curves We resume here the process of curve reconstruction (see [9, 10] for a complete description). The data given by Morphosense are the separation of the sensors along the curve (the ribbon) and orientation of these sensors in space. Since we have length constraints (sensors separation), we work with the arc-length parameter for the reconstruction of the spatial curves. The derivative is thus a unit vector representing a curve lying on the unit sphere and the data are samples of this derivative at known parameters. Then, in order to interpolate this data by a curve staying on the sphere, we use cubic splines on the unit sphere, as proposed in [6].
Fig. 12a. Initial surface with 4 evaluated geodesics numbered from 1 to 4.
Thus, the method proceeds in two steps: first, we reconstruct the derivative function, and then we integrate this derivative to obtain the desired 3D curve. A numerical method is used for this second part. An example of such a spatial reconstruction is presented in Fig. 15. Left, is shown the original curve (thin) and the reconstructed one (bold). Right, is the derivative reconstructed on the sphere with values of sensors (dots). 5.2. Theoretical analysis A developable surface D is isometric to a plane, so that any geodesic curve γ on D can be associated with a straight line on a plane. Thus, considering any regular surface F, if we assume that the two surfaces D and F have a G 1 contact along a common curve γ and that γ is a geodesic curve on D, we will show that γ is a geodesic curve on the surface F. In practice, to ensure this G 1 contact we consider a ribbon D on the surface D, of width 2, centered on the curve γ . One can understand that even with a thin ribbon the contact is theoretical on hyperbolic points of the surface F. Then, using the fact that geodesic curves are analogous to straight lines in the plane, we will give a mechanical characterization of geodesics related with the G 1 installation of the ribbon on F. Precisely, Let (τ (s), V (s), N (s)) the Darboux frame at point γ (s) on the developable surface D . When installing the ribbon on the surface F, this one is deformed in an isometric way which only authorize local rotations around the vectors τ (s) and V (s) along the curve γ . The idea is that geodesics have to be straight, so that local rotations around the vector N (s) are forbidden. The curve γ on the ribbon D will be said N -rigid. Definition 1. A regular curve γ (s) running on a surface F is said to be N -rigid if its Darboux frame can not locally rotate (with parameter s) around the vector N (s). Proposition 1. Let γ (s) be a N -rigid regular curve running on a surface F. Then γ is a geodesic curve of the surface F. Proof. Since the curve γ is N -rigid, the Darboux frame (τ (s)+ dτ, V (s) + dV, N (s) + dN ) at point γ (s + ds) can only be deduced from the Darboux frame at point γ (s) by rotations around vectors τ (s) and V (s). Thus, considering elementary rotations around vector τ (s) by an angle dωx and around vector V (s) by an angle dω y , we can write τ (s) + dτ cos dω y 0 − sin dω y V (s) + dV = 0 1 0 N (s) + dN sin dω y 0 cos dω y
Fig. 12b. Reconstructed surfaces without (left) and with (right) smoothing of the binormal vectors.
489
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
Fig. 13a. Binormal vectors along the four evaluated geodesics.
Fig. 13b. Binormal vectors after smoothing.
1 0 × 0 cos dωx 0 sin dωx
0 τ (s) − sin dωx V (s) , cos dωx N (s)
which, assuming cos dωx ' 1 ' cos dω y , sin dωx ' dωx , sin dω y ' dω y , leads to dτ 0 0 −dω y τ (s) dV = 0 0 −dωx V (s) . (6) dN dω y dωx 0 N (s) Thus, comparing this result with Darboux relations (3), we get k g (s) ≡ 0 which proves that γ is a geodesic curve. Notice that we get exactly the same result if we first apply an elementary rotation around V (s). Proposition 2. Let γ (s) be a geodesic curve of the surface F. Then γ is N -rigid. Proof. The Darboux frame (τ (s)+dτ, V (s)+dV, N (s)+dN ) at point γ (s + ds) can be deduced from the Darboux frame at point γ (s) by 3 elementary rotations Rdωx , Rdω y , Rdωz
around vectors τ (s), V (s), N (s) by angles dωx , dω y , dωz . Thus, assuming cos dωx ' 1, cos dω y ' 1, cos dωz ' 1, sin dωx ' dωx , sin dω y ' dω y , sin dωz ' dωz , we have τ (s) + dτ τ (s) V (s) + dV = Rdωz ◦ Rdω y ◦ Rdωx V (s) , (7) N (s) N (s) + dN which leads to dτ 0 dV = −dωz dN dω y
−dωz 0 dωx
−dω y τ (s) −dωx V (s) . 0 N (s)
(8)
Now suppose that γ is a geodesic curve, that is k g (s) ≡ 0. Thus, from Darboux relations (3) we deduce that dωz = 0, which proves that the curve γ is N -rigid. Notice that the result does not depend on which order we apply the 3 elementary rotations. Proposition 3. Consider two regular oriented surfaces F1 and F2 with a G 1 contact along a common curve γ . Then, if γ is
490
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
Fig. 14. “Naked” Ribbon Morphosense with sensors and below, a real time demonstration.
the ribbon and F along the common curve γ ensures that the Darboux frames of γ on each of these two surfaces agree (regardless of the orientation). So that the curve γ is also N rigid on F. Then, from Proposition 1, we deduce that γ is a geodesic curve on F. Notice that Proposition 3 provides a more direct way to prove this result. Nevertheless, Propositions 1 and 2 connect the mechanical properties of the ribbon to its geometric behaviour. 5.3. Numerical results Fig. 15. Spatial curve reconstruction.
a geodesic curve on one surface, it is a geodesic curve on the other surface. Proof. Let us denote by (τi (s), Vi (s), Ni (s)), i = 1, 2, the Darboux frames of the curve γ on each one of the two surfaces Fi . Of course we have τ1 (s) = τ2 (s) for each value of the parameter s and since these two surfaces have a G 1 contact along the curve γ , we have N1 ≡ ±N2 . If necessary we can modify the orientation of one the two surfaces to get N1 ≡ N2 , so that the two Darboux frames are the same all along the curve γ . Finally, the result can be deduced from Darboux relations (3). From these results we are now able to verify that our ribbon actually follows geodesic curves on the surface F. The ribbon is inextensible and not “flexible”. It has the same behaviour as a thin strip of paper, i.e., as a developable surface. Consider a geodesic curve γ on the ribbon (this geodesic is isometric to a straight line of the plane). From Proposition 2, we deduce that this curve γ is N -rigid on the ribbon. Then, using the same arguments as in Proposition 3, the G 1 contact between
We first consider a numerical model (see Fig. 16, top) of a real surface (in Fig. 17). Five geodesic curves, numbered from 1 to 5, are evaluated and then the reconstruction is achieved via geodesic interpolation. In Fig. 16, each surface Si j is the interpolating surface between geodesics i and j and S12345 is the whole interpolating surface. Figure 17(left), shows the installation of the ribbon on the real surface, and (right) is the reconstructed surface via geodesic interpolation. 6. Conclusion In this paper it has been shown that numerical and real surfaces can be reconstructed from some geodesic running on this surface. These geodesics can be captured using a ribbon of micro-sensors so that the reconstruction process is really intrinsic. Even if this point is not at the heart of this paper, our experiments show data given by micro-sensors to be very sensitive. Thus, our process requires smoothing of this data. Nevertheless, the reconstruction is simple and fast and so could be suitable for deformation of surfaces.
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
Fig. 16. Reconstruction via geodesic interpolation.
Fig. 17. Reconstruction from ribbon’s data representing transversal geodesic curves.
491
492
N. Sprynski et al. / Computer-Aided Design 40 (2008) 480–492
Acknowledgments We are very grateful for the anonymous reviewers for their detailed constructive comments. References [1] Coons S. Surface patches and B-spline curves. In: Barnhill R, Riesenfeld R, editors. Computer aided geometric design. Academic Press; 1974. [2] Coons S. Surfaces for computer aided design, Technical report, M.I.T. 1964. Available as AD 663 504 from the National Technical Information Service. Springfield, VA 22161. [3] Do Carmo MP. Differential geometry of curves and surfaces. PrenticeHall; 1976. [4] Farin G. Curves and surfaces for CAGD, 5th ed. Academic Press; 2002. [5] Gordon W. An operator calculus for surface and volume modeling. IEEE
Computer Graphics and Applications 1983;3:18–22. [6] Nielson GM. ν-quaternion splines for the smooth interpolation of orientations. IEEE Transactions on Visualization and Computer Graphics 2004;10(2):224–9. [7] Pham-Trong V. (2001) D´etermination g´eom´etrique de chemins g´eod´esiques sur des surfaces de subdivision. Th`ese de l’universit´e Joseph Fourier, Math´ematiques Appliqu´ees. Laboratoire LMC-IMAG. 28 Septembre 2001. [8] Pham-Trong V, Biard L, Szafran N. Pseudo-geodesics on threedimensional surfaces and pseudo-geodesic meshes. Numerical Algorithms 2001;26:305–15. [9] Sprynski N, Lacolle B, Biard L, David D. In: Chenin P, Lyche T, Schumaker LL, editors. Curve and surface reconstruction via tangential informations, curve and surface design, Avignon 2006. Nashboro Press; 2007. p. 254–63. [10] Sprynski N. Reconstruction de courbes et surfaces a` partir de donn´ees tangentielles, Laboratoires CEA/LETI, LJK. Th`ese de l’universit´e Joseph Fourier, Grenoble. 5 Juillet 2007.
