A research note on design of fair surfaces over irregular domains using data-dependent triangulation

Available online at www.sciencedirect.com

Applied Mathematical Modelling 32 (2008) 2172–2195 www.elsevier.com/locate/apm

A research note on design of fair surfaces over irregular domains using data-dependent triangulation R. Sharma *, O.P. Sha

1

Design Laboratory, Department of Ocean Engineering and Naval Architecture, Indian Institute of Technology, Kharagpur (WB) 721 302, India Received 1 February 2006; received in revised form 1 September 2006; accepted 4 July 2007 Available online 20 July 2007

Abstract Design of fair surfaces over irregular domain is a fundamental problem in computer aided geometric design (CAGD), and has applications in engineering sciences (i.e. aircraft science, automobile science and ship science etc.). In design of fair surfaces over irregular domain defined over scattered data it was widely accepted till recently that one should use Delaunay triangulation because of its global optimum property. However, in recent times it has been shown that for continuous piecewise polynomial parametric surfaces improvements in the quality of fit can be achieved if the triangulation pattern is made dependent upon some topological property of the data set or is simply data dependent. The smoothness and fairness of surface’s planar cuts is important because not only it ensures favorable hydrodynamic drag, but also helps in reducing manhours during the production of the surface. In this paper we discuss a method for construction of C1 piecewise polynomial parametric fair surfaces which interpolate prescribed R3 scattered data using spaces of parametric splines defined on R3 triangulation. We show that our method is more specific to the cases when the projection on 2-D plane may consist of triangles of zero area. The proposed method is fast, numerically stable and robust, and computationally inexpensive. In the present work numerical examples dealing with surfaces approximated on standard curved plates, and ship hull surface have been presented.  2007 Elsevier Inc. All rights reserved. Keywords: Cubic spline; C1 continuous surface; Data-dependent triangulation; Minimum energy surface; Scattered data interpolation; Surface fitting

1. Introduction Surface interpolation over a scattered data on R3 is a problem of general interest cutting across many engineering disciplines. In recent years, this problem has been of considerable interest to a lot of researchers, e.g. [1,2], and references therein. In general, for modeling problems of surface interpolation, piecewise

*

1

Corresponding author. Tel.: +91 3222 281601; fax: +91 3222 282700. E-mail addresses: [email protected] (R. Sharma), [email protected] (O.P. Sha). Tel.: +91 3222 283788; fax: +91 3222 277190.

0307-904X/$ - see front matter  2007 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2007.07.008

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

2173

polynomial parametric surfaces are constructed which interpolate prescribed R3 scattered data using spaces of parametric splines defined on R2 triangulation. The points in R3 are projected onto a plane in R2 , and a triangulation is defined as a set of non-intersecting, and stitched triangles with the vertices at the points in the point set. The resulting surface S defined by the union of triangles in F by interpolating parametric surface of desired degree d and continuity r. A unique plane can always be defined which passes through the vertices of the triangle and the vertices of the triangle projected on the plane in R2 . However, if one wants to use the same plane for the projection of the complete point set in R3 , a situation might arise where the triangle in 3D may either be transformed into a zero area in 2D or a triangle with poor aspect ratio in 2D. Therefore, the quality of the fit, for such an application is dependent on the projection of the scattered R3 data on the plane in R2 . In general, this is undesired and may adversely affect the quality of the surface. Generally, in surface interpolation a spline of desired degree and continuity is passed over a triangulation. According to Lodha and Franke [3], for scattered data interpolation classical Delaunay triangulation can be used. However, this is suitable only when the triangulation as a parameter is used for different functions. In more specific cases (i.e. surfaces of mixed continuities) triangulation as a parameter cannot be used for different functions and the choice of triangulation depends upon the function being fitted. In this context, Dyn et al. [4,5] have observed that the quality of the fit can be improved by adjusting the triangulation to the function being fitted. The authors have also discussed that the best quality is achieved by minimizing the size of the normal derivative discontinuities across the edge of the triangulation. In the same line, Quak and Schumaker [6] have discussed the possibility of applying the construction of piecewise polynomial parametric surfaces which interpolate prescribed R3 scattered data using spaces of parametric splines defined on R2 triangulation. The authors use projection of 3-D data set on 2-D along a plane. However, they have not given any guidelines for the selection of the plane. Since, the quality of fit in case of surface design not only depends upon the continuities (i.e. geometric G0, or G1, or G2; or derivational C0, or C1, or C2), but to have practical applications in the engineering industries also upon the desired fairness properties of the surface. To have a fair surface, it has been suggested by many authors to incorporate principles of physics or optimization techniques in the design process to automatically achieve a desired smooth surface, e.g. [7–11], and references therein. In order to achieve a minimum value for an objective function, one has to solve a simultaneous system of linear, or non-linear equations. The solution of non-linear equations in general is too expensive to be used in an interactive user environment in engineering sciences, and sometimes no acceptable solution can be computed. In the present work we explore the possibility of applying the construction of piecewise polynomial parametric surfaces which interpolate prescribed R3 scattered data using spaces of parametric splines defined on R3 triangulation. Since the triangulation is defined in R3 no projection is required. We use piecewise cubic C1 surfaces and the selection of the triangulation is based on the idea of minimizing the energy of the resulting surface. We use this objective function (i.e. elastic energy of a thin plate) because it can be approximated numerically, resulting in stable and robust computations, and better computational efficiency. The remaining of this paper is organized as follows: Section 2 presents briefly about scattered data interpolation, and present work; Section 3 presents surface representation; Section 4 presents the definition of the energy of the surface; Section 5 presents the numerical examples, and discussion; the Section 6 concludes the paper by identifying some future applications, and further scope of research. Additionally, the proof of the theorem is presented in Appendix A, the computation of the energy of a piecewise triangular polynomial patch is presented in Appendix B, and some basic transformations that are used in our formulation are presented in Appendix C. 2. Brief description about scattered data interpolation, and present work The basic problem in scattered data interpolation problem is: Given a finite set of N points in R3 , find surface that interpolates the given set of points. This type of problem arises in many areas such as CAD/CAM, geometric modeling, graphics, and scientific visualization. The problem can be viewed in trivariate (i.e. for volumetric or solid geometric applications, for a comprehensive survey see [12–17]) or in bivariate (i.e. for area centric or surface applications). In general the solutions to the problem of scattered data interpolation in

2174

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

bivariate forms which are of interest in the present work, are classified into piecewise polynomial or rational parametric solutions (including interpolation at both local, and global level with a single polynomial or group of polynomials; interpolants based on data-dependent triangulations; piecewise linear solutions such as alphashapes; and interpolants on irregular meshes, for a comprehensive survey see [17–24]), algebraic solutions (for a comprehensive survey see [25]), radial basis function methods (including methods based upon Hardy’s multiquadrics; inverse multi-quadrics; and thin plate splines; for a comprehensive survey see [26–28]), Shepard’s methods (including methods based upon local solution blending; and natural neighbor interpolants, for a comprehensive survey see [29,30]), and subdivision surfaces (including methods based upon Catmull–Clark subdivisions, for a comprehensive survey see [31]). Normally, in academics and industry, the piecewise continuous polynomials or rational parametric spline solutions are used for the problem of scattered data interpolation. Theoretically, the construction of a ‘unique global polynomial interpolant’ to a given scattered data is possible, but the multivariate interpolants will be having a high degree for large data sets which is computationally expensive. Therefore, for better computational efficiency piecewise polynomial or rational interpolants of low degrees are required. Topologically, for a rectangular data set Tensor product bicubic B-splines and Non-Uniform Rational B-spline Surfaces (NURBS) give good results but since they are inherently rectangular they suffer from ‘degeneracy’ in the cases of general scattered data points on arbitrary domains. In general, for general scattered data points given on arbitrary domains a typical method is to triangulate the domain data points in the space and then construct piecewise continuous interpolants on each triangle. This method is classified as a ‘triangulation-based method’. In this class of solutions a simple solution is to construct a Delaunay triangulation of the scattered data points in R2 and construct the piecewise linear interpolant on this triangulation. Though, Delaunay triangulation has many nice properties (i.e. for example it can be computed in O(N log N) time, etc.; for details see [32]), but depending upon the nature of the data its interpolant may not necessarily have the same nice properties. To rectify this problem the piecewise linear interpolants that are built on triangulations that depend upon geometrical or topological properties of data [6,33,34] have been explored. In the present work we explore a data-dependent triangulation whose interpolant is smooth and it result into a tangent plane (i.e. first order continuous derivatives or C1 continuous) continuous surface. We focus on to construct piecewise C1 continuous interpolants over ‘triangulated’ scattered data using only minimum number patches of low degree (i.e. d = 3) as are possible depending upon the data set. We explore fitting piecewise C1 continuous cubic polynomial patches by splitting each triangle into three triangles, known as Clough–Tocher [35] split. The resulting interpolants are local and any change in a data value affects the interpolant only in a neighborhood of that data point. Additionally, we measure the quality of interpolants by energy and a better quality in the fitting is achieved via minimization of a variational function. In engineering design, a surface’s planar cuts (i.e. intersection curves of the surface with X = constant, Y = constant, and Z = constant planes) play an important role in the shape interrogation process. For example from manufacturing point of view it is highly desired that the surface has linear planar cuts as they have low cost of manufacturing, and from hydrodynamic/aerodynamic points of view it is highly desired that the surface has non-linear planar cuts as they will allow smoother flow resulting in better hydrodynamic/aerodynamic performance. In the present work, we explore the non-linearity in the surface’s planar cuts and specifically we investigate the relationships between non-linearity in the surface’s planar cuts with triangulation patterns, and surface’s energy measures. 3. Surface representation The problem addressed in this work is formally stated as: Given a set of points P ¼ fðxi ; y i ; zi ÞgNi¼1 in R3 and N a corresponding set of measurements fGi ¼ f ðxi ; y i ; zi Þgi¼1 on a function f, construct a surface S approximating f. The function f is used for generating additional points for triangulation, and as will be discussed in Section 4 is the strain energy function of thin plates under small deflections in our formulation. Following [6], let D is a triangulation associated with a set P of points as defined previously. In addition suppose the gradients (fx, fy, fz) have been estimated at each of the data points of the data set (i.e. (xi, yi, zi)), producing the values ðGxi ; Gyi ; Gzi Þ, i = 1, . . . , N. Let T be a typical triangle in the triangulation D, and its vertices

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

2175

be P1, P2, P3 and its edges are e1, e2, e3. By connecting the centroid point P 4 ¼ ðP 1 þP32 þP 3 Þ of the triangle T to each of its vertices, we split T into three sub-triangles as shown in Fig. 1. We denote by T1, T2, T3 the three sub-triangles obtained by connecting the centroid (i.e. P4) and the vertices P1, P2, P3. This split is known as Clough–Tocher split. As per Renka [36], and Lawson [37], there is a unique piecewise cubic function s1 in C1(T) which satisfies the interpolation condition as: sT ðP i Þ ¼ Gi ;

Dx sT ðP i Þ ¼ Gxi ;

Dy sT ðP i Þ ¼ Gyi ;

Dz sT ðP i Þ ¼ Gzi

ð1Þ

for i = 1, 2, 3 and which has the further property that the normal derivatives along each of the edges e1, e2, e3 are linear polynomials rather than quadratics. Since this property has been combined in the interpolation conditions it implies that the surface S obtained by stitching together the patches sT is globally C1 and satisfies the interpolation conditions as imposed in Eq. (1) for all i = 1, . . . , N. Theorem 1. On the triangle T, the spline has the following Be´rnstein–Be´zier coefficients (Fig. 1): c1 ¼ G 1 ; c2 ¼ G 2 ; c3 ¼ G 3 ;   ðx2  x1 Þ  Gx1 þ ðy 2  y 1 Þ  Gy1 þ ðz2  z1 Þ  Gz1 c4 ¼ ; 3 þ G1   ðx4  x1 Þ  Gx1 þ ðy 4  y 1 Þ  Gy1 þ ðz4  z1 Þ  Gz1 c5 ¼ ; 3 þ G1   ðx3  x1 Þ  Gx1 þ ðy 3  y 1 Þ  Gy1 þ ðz3  z1 Þ  Gz1 c6 ¼ ; 3 þ G1   ðx3  x2 Þ  Gx2 þ ðy 3  y 2 Þ  Gy2 þ ðz3  z2 Þ  Gz2 c7 ¼ ; 3 þ G2   ðx4  x2 Þ  Gx2 þ ðy 4  y 2 Þ  Gy2 þ ðz4  z2 Þ  Gz2 c8 ¼ ; 3 þ G2   ðx1  x2 Þ  Gx2 þ ðy 1  y 2 Þ  Gy2 þ ðz1  z2 Þ  Gz2 c9 ¼ ; 3 þ G2   ðx1  x3 Þ  Gx3 þ ðy 1  y 3 Þ  Gy3 þ ðz1  z3 Þ  Gz3 c10 ¼ ; 3 þ G3   ðx4  x3 Þ  Gx3 þ ðy 4  y 3 Þ  Gy3 þ ðz4  z3 Þ  Gz3 c11 ¼ ; 3 þ G3   ðx2  x3 Þ  Gx3 þ ðy 2  y 3 Þ  Gy3 þ ðz2  z3 Þ  Gz3 c12 ¼ ; 3 þ G3 P1 c1 c5 c6

c4 c16 c13 c9 c17

P2 c 2

c8

c15

P4

c19

c10

c18

c14

c7

c11

c12

Fig. 1. The Clough–Tocher split.

c3 P 3

2176

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

c13 c14 c15 c16 c17 c18 c19

  ðc5 þ c8 þ ðh1  1Þ  c1 Þ þ ð2  3h1 Þ  c4 þ ð3h1  1Þ  c9  h1  c2 ; ¼ 2   ðc8 þ c11 þ ðh2  1Þ  c2 Þ þ ð2  3h2 Þ  c7 þ ð3h2  1Þ  c12  h2  c3 ¼ ; 2   ðc11 þ c5 þ ðh3  1Þ  c3 Þ þ ð2  3h3 Þ  c10 þ ð3h3  1Þ  c6  h3  c1 ¼ ; 2 hc þ c þ c i 15 5 13 ¼ ; 3 hc þ c þ c i 13 8 14 ¼ ; 3 hc þ c þ c i 14 11 15 ¼ ; 3 hc þ c þ c i 18 16 17 ¼ ; 3

where, " h1 ¼

2

" h2 ¼ " h3 ¼

# ðx4  x1 Þ  ðx2  x1 Þ þ ðy 4  y 1 Þ  ðy 2  y 1 Þ þ ðz4  z1 Þ  ðz2  z1 Þ 2

2

ððx2  x1 Þ þ ðy 2  y 1 Þ þ ðz2  z1 Þ Þ

# ðx4  x2 Þ  ðx3  x2 Þ þ ðy 4  y 2 Þ  ðy 3  y 2 Þ þ ðz4  z2 Þ  ðz3  z2 Þ ððx3  x2 Þ2 þ ðy 3  y 2 Þ2 þ ðz3  z2 Þ2 Þ

# ðx4  x3 Þ  ðx1  x3 Þ þ ðy 4  y 3 Þ  ðy 1  y 3 Þ þ ðz4  z3 Þ  ðz1  z3 Þ 2

2

2

ððx1  x3 Þ þ ðy 1  y 3 Þ þ ðz1  z3 Þ Þ

; ; ;

hx þ x þ x i 1 2 3 ; 3 hy þ y þ y i 2 3 ; y4 ¼ 1 3 hz þ z þ z i 1 2 3 : and z4 ¼ 3 x4 ¼

The proof of Theorem 1 is given in Appendix A.

4. Definition of the energy of the surface As per Reddy [38], and Timoshenko and Woinowsky-Krieger [39]; the elastic energy of a thin plate of arbitrary shape under small deflections is given as: Z Z Z Z "Z Z 2 Z 2 # þ  2ð1  mÞ  du dv; ð2Þ uu

vv

uu

vv

uv

where the parameter m is the Poisson’s ratio depending on the material at hand. The value of m for steel or aluminum is 0.3. Though, an obvious setting of m is at 0.3, but we use m = 0 to simplify (i.e. as will be shown later it result into a symmetric formulation which is computationally efficient) the formulation. In our work we use Eq. (2) for two purposes. First we use f to be equivalent to elastic energy, and for the given point set N P ¼ fðxi ; y i ; zi Þgi¼1 we generate additional points along a rectangular or triangular grid (i.e. if the number of input points is even for example 4, 6, 8, . . . , 4n, and forms a rectangular pattern then rectangular grid; if the number of input points is odd for example 3, 5, 7, . . . , 2n + 1, and forms a triangular pattern then triangular grid; where n is an positive non-zero integer) for a cubic triangular patch by fitting f over P. Later, with these additional points we generate the surface by fitting cubic spline over the Delaunay triangulation of the

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

2177

rectangular or triangular grid point set. Second, we generate the C1 surface by further refining the triangulation using minimization of f. N Now, following [6], and [8]; suppose that D ¼ fT i gi is the first refined triangulation of D obtained by performing Clough–Tocher split on each of the triangles of D. Let s is a piecewise polynomial function over D. If we think of a surface S as a thin plate (assuming Poisson’s ratio is equal to zero), the natural elastic energy measure of the surface is given by, f ðT i Þ ¼ oðsÞ ¼

N X

oi ðsÞ;

i¼1

where for each triangle Ti of D, Z Z 2 2 2 oT i ðsÞ ¼ ½ðuxx Þ þ ðuyy Þ þ 2  ðuxy Þ  dx dy;

i ¼ 1; . . . ; N :

ð3Þ

Ti

In Eq. (3) for numerical computation s is represented in classical Be´rnstein–Be´zier coefficients. Each cubic piece si ¼ sI T i is described by a vector ci ¼ ðci;1 ; ci;2 ; . . . ; ci;10 Þ. Normally, the authors work with Eq. (3) either via non-linear numerical energy expressions or they use numerical integration techniques to evaluate the quadratic energy integral, and both of these are computationally quite expensive. Additionally, because of high computational time it restricts the use of Eq. (3) to only few cases where the point set is comparatively small. Quak and Schumaker [40] have shown that the energy contribution of this piece can be approximated by the quadratic form: T

ð3Þ

dT i ðsÞ ¼ ffcgi  fEgi  fcgi g;

ð4Þ

fEgið3Þ

where is a symmetric 10 · 10 energy matrix. However, the formulations given in [40] are for 2-D canonical triangle but these can be further extended to 3-D canonical triangle (for proof, and details see [41]). Furthermore, it has been shown that Eq. (4) is not only numerical stable, but also has low computational time. On a comparative scale, it reduces the computational time by around 40% [41]. Moreover, because of low computational time it allows the use of Eq. (3) to general cases where the point set is comparatively bigger. The numerical computation of the energy of a piecewise polynomial surface patch is given in Appendix B. 5. Numerical examples We present two algorithms for producing data-dependent triangulation based upon idea of minimizing the energy of the corresponding piecewise C1 continuous surface. We consider triangles in our implementations. The modification in the triangulation is done by swap test and is completely local in nature. Hence to compute the energies in the swap test only the three vertices of the triangle are used. Algorithm 1 N

1. Define an initial triangulation F using the given point set P ¼ fðxi ; y i ; zi Þgi¼1 . 2. Generate the additional points by fitting f over P in a rectangular or triangular grid. 3. Compute Delaunay triangulation over the domain of additional points in a rectangular or triangular grid as defined in Step (2). 4. Generate the surface by fitting a spline over Delaunay triangulation as computed in Step (3).

Algorithm 2 N

1. Define an initial triangulation F using the given point set P ¼ fðxi ; y i ; zi Þgi¼1 . 2. Generate the additional points by fitting f over P in a rectangular or triangular grid. 3. Compute Delaunay triangulation over the domain of additional points in a rectangular or triangular grid as defined in Step (2). 4. Split triangles in Delaunay triangulation using Clough–Tocher split.

2178

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

Fig. 2. Basic triangulation (F) over the given point set for non-flat rectangular plate with one corner raised (top view).

Fig. 3a. Additional point set and surface grid for a non-flat rectangular plate with one corner raised (top view).

Fig. 3b. Additional point set and surface grid for a non-flat rectangular plate with one corner raised (isometric view).

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

2179

5. Compute the energies of the each ‘triangle’ in the triangulation as defined in Step (3). 6. Arrange the energies into a vector o1 P    onp . 7. For i = 1, . . . , np, If a split (Clough–Tocher) in ‘triangle i’ is possible with a reduction in energy then split the triangle and return to Step 3. Otherwise stop. 8. Generate the surface by fitting a spline over triangulation as computed in Step (6). In implementations of Algorithms 1 and 2 as described above any standard Delaunay triangulation algorithm [2,36] can be used. In our implementations we have used subroutines of Delaunay triangulation algorithm [36] for computing the gradients at each point of the given point set, and for Algorithm 2 we have used insertion algorithm [36], and integrated that with other steps of Algorithm 2. Both the Algorithms 1 and 2 have been implemented in C++ on a on a Silicon GraphicsTM* OriginTM* 200 workstation. The surface views are presented in AutoCADTM** R2006. To illustrate the two algorithms we consider two problems. For the sake of reproducibility we give the explicit co-ordinates (i.e. Pi = (Xi, Yi, Zi)), and geometric definition with the continuities imposed.

(a) top view

(b) isometric view Fig. 4. Surface for a non-flat rectangular plate with one corner raised on Delaunay triangulation, C1 continuous surface.

2180

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

Problem 1. This is a simple case where we have used two patches to model a rectangular plate in which one corner has been raised. Let P1 = (0, 0, 1), P2 = (1, 1, 1), P3 = (0, 1, 0) and P4 = (1, 0, 0). In the initial triangulation to start with, we have divided the rectangle along the line about which the curvature is changing. Suppose the indices of the vertices of the initial triangulation Fare given by (1, 2, 3) and (3, 2, 4). Discussion: Case 1. In this case, we start with the given point set, and define the initial triangulation as shown in Fig. 2. The orientation of the Figs. (i.e. 2, 3a, 3b, 4a, 4b, 5a, 5b, 6a, 6b, 7a, 7b, 8a, 8b, 9a, 9b, 10a, 10b, 11a, and 11b) has been defined at (0, 1, 0), for X, Y, and Z axis. As defined in Algorithm 1, since the number of points in the point set is 4 and they make a rectangular pattern the additional points are generated by fitting f over P in a rectangular grid. The rectangular grid is shown in Figs. 3a and 3b. The rectangular grid of Fig. 3 is then triangulated with Delaunay triangulation, and we have fitted C1 continuous surface. The spline surface is shown in Fig. 4a and 4b. This surface can be described as the basic surface and it is used to study the behavior in other cases. This surface has energy of 53.7 and an error (in discrete maximum norm) of 0.067. In Fig. 4a, the curvature along edges is clear. Edges with greater curvature implies higher inherent energy in the surface, and hence less fairness. Since here the modeling is done without any minimization of energy, this is obvious. Fig. 4b shows an isometric view in which the curved sides, the twist in surface and the higher curvatures are visible.

(a) top view

(b) isometric view Fig. 5. Intersection curves of constant X planes with surface for a non-flat rectangular plate with one corner raised on Delaunay triangulation, C1 continuous surface.

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

2181

To further analyze the surface and explore what does high energy means in designing engineering bodies, we compute the intersection curves of the surface with constant X, Y, and Z planes at the interval of 0.1, from 0.1 to 0.9. To compute surface/plane intersection curve, we have used the method of Sharma and Sha [42]. And for the sake of clarity in this example under both cases 1 and 2, we plot the intersection curves onto the surface. Figs. 5a and 5b, 6a and 6b, and 7a and 7b, show the two views of the intersection curves at constant X, Y, and Z planes respectively. The surface is non-flat but is C1 continuous, so all the intersection curves are continuous without any break in continuity. Since the surface is a C1 continuous surface (i.e. intersection curve with continuous first order derivative at X, Y, and Z = constant plane), the intersection curves are continuous across the curve of C1 continuity (i.e. C1 continuity is imposed across the interior edge and that edge is a curve of C1 continuity). And, any section taken across this curve of continuity will show that it is across this continuous curve that the surface is C1 continuous. Though the intersection curves are continuous but because the edges are curved, and the inherent energy in the surface is high, the intersection curves are very much

Fig. 6a. Intersection curves of constant Y planes with surface for a non-flat rectangular plate with one corner raised on Delaunay triangulation, C1 continuous surface (top view).

Fig. 6b. Intersection curves of constant Y planes with surface for a non-flat rectangular plate with one corner raised on Delaunay triangulation, C1 continuous surface (isometric view).

2182

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

(a) top view

(b) isometric view Fig. 7. Intersection curves of constant Z planes with surface for a non-flat rectangular plate with one corner raised on Delaunay triangulation, C1 continuous surface.

non-linear. The non-linearity in the intersection curves mean high Gaussian curvature and that implies high cost of production (i.e. straight sections can be manufactured by rolling, pressing, and bending, and non-linear sections can only be manufactured with application of stress via line heating) in metal-based industries like ships, automobiles, and aircrafts. The non-linearity in the intersection curves is clearly visible in Figs. 5a and 5b, 6a and 6b, and 7a and 7b. Furthermore, since the surface is free from undesired oscillations so it is clear from Fig. 4a and 4b that the surface is a fair surface. Case 2. In this case, we start with the same given point set, and initial triangulation as shown in Fig. 2. The rectangular grid as shown in Figs. 3a, and 3b, is then triangulated with data-dependent triangulation via minimization of energy as defined in Algorithm 2. And we have fitted C1 continuous surface. The spline surface is shown in Fig. 8a and 8b. This surface has energy of 23.7 and an error (in discrete maximum norm) of 0.079. From Fig. 8a and 8b it is clear that the edges are less curved as compared to Figs. 4a and 4b. This surface has low inherent energy as can be seen from low curvatures along the edges which is evident from the fact that the modeling has been done by minimization of energy. However, the pattern of triangulation is similar to Fig. 4, which might be because of the reason that initially the triangulation has been defined over too many points and that may lead to a global optimal solution. If the initial triangulation is near to global optimal then there

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

2183

(a) top view

(b) isometric view Fig. 8. Surface for a non-flat rectangular plate with one corner raised on data-dependent triangulation (minimal energy), C1 continuous surface.

is less scope of refinement. Again this requires further investigation. An edge with low curvature implies low inherent energy in the surface, and hence high fairness. Fig. 8b shows the isometric view in which the less curved sides, twist in surface and the lower curvatures are visible. As previously, we compute the intersection curves of the surface with constant X, Y, and Z planes at the interval of 0.1, from 0.1 to 0.9. Figs. 9a and 9b, 10a and 10b, and 11a and 11b, show the two views of the intersection curves at constant X, Y, and Z planes respectively. Again, as earlier, the surface is non-flat but is C1 continuous, so all the intersection curves are continuous without any break in continuity. Compared to case 1, in Fig. 8a and 8b, the edges are less curved, and the inherent energy in the surface is low so there is a clear reduction in the non-linearity of the intersection curves. As a cautious observation it can be stated that release of bending energy in a surface basically implies that the edges will be straightened, and that straightening of edges will open the surface reducing its inherent energy. The low non-linearity in the intersection curves mean low Gaussian curvature and that implies low cost of production in metal-based industries. The low non-linearity in the intersection curves is clearly visible in Figs. 9a and 9b. Furthermore, the surface is free from undesired oscillations, and it is clear from Fig. 8a and 8b that the surface is fair surface.

2184

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

Fig. 9a. Intersection curves of constant X planes with surface for a non-flat rectangular plate with one corner raised on data-dependent triangulation (minimal energy), C1 continuous surface (top view).

Fig. 9b. Intersection curves of constant X planes with surface for a non-flat rectangular plate with one corner raised on data-dependent triangulation (minimal energy), C1 continuous surface (isometric view).

Problem 2. This example is from shipbuilding industries, and it is a case of a patch taken on the forward part of a hard chine fishing vessel. It is well known [43], that hard chine vessels and boats are complex objects, and their computer aided design and interrogation is challenging. Let P1 = (13.0, 3.275, 3.53), P2 = (13.0, 2.990, 0.775), P3 = (18.0, 2.672, 3.855), P4 = (18.0, 2.226, 1.512), P5 = (21.5, 1.0904, 4.202), P6 = (21.5, 0.5322, 2.58), P7 = (22.306, 0.4655, 4.2797), P8 = (22.306, 0.0, 2.90) and P9 = (22.8, 0.0, 4.35). In the initial triangulation to start with we have divided the rectangle along the line about which the curvature is changing. Suppose the indices of the vertices of the initial triangulation F are given by (1, 2, 3), (2, 3, 4), (3, 4, 5), (4, 5, 6), (5, 6, 7), (6, 7, 8), and (7, 8, 9). Discussion: As previously, in this case also we start with the given point set, and define the initial triangulation as shown in Fig. 12. The orientation of the Figs. (i.e. 12, 13a, 13b, 14, 15, and 16) has been defined at (13.0000, 2.99, 0.7750), for X, Y, and Z axis.

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

2185

(a) top view

(b) isometric view Fig. 10. Intersection curves of constant Y planes with surface for a non-flat rectangular plate with one corner raised on data-dependent triangulation (minimal energy), C1 continuous surface.

As defined in Algorithm 1, since the number of points in the point set is 9 and they make a triangular pattern the additional points are generated by fitting f over P in a triangular grid. The triangular grid is then triangulated with data-dependent triangulation via minimization of energy as defined in Algorithm 2. The spline surface is shown in Fig. 13a and b. This surface has energy of 33.7 and an error (in discrete maximum norm) of 0.076. In Fig. 13a, no higher curvature along edges is visible. An edge with no higher visible curvature implies very low inherent energy in the surface, and hence good fairness. Since here the modeling is done with minimization of energy, this implies that the designed fair surface has good quality and it can be used in industry. Fig. 13b shows the isometric view in which also there is no visible twist, or high curvature along the edges. The surface is non-flat but is C1 continuous, so the intersection curves will be non-linear but continuous. As previously, to further analyze the surface and explore does the surface have industrially relevant features to be used in shipbuilding industry for designing and analyzing ships, we compute the intersection curves of the surface with constant X, Y, and Z planes. The intersection curves of the surface with constant X planes are computed at X = 13.0000, 14.0000, 15.0000, 16.0000, 17.0000, 18.0000, 19.0000, 20, 21.0000, 21.5000, 21.7015, 21.9030, 22.1045, 22.3060, 22.4295, 22.5530, and 22.6765, and these are shown in Fig. 14. The intersection

2186

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

(a) top view

(b) isometric view Fig. 11. Intersection curves of constant Z planes with surface for a non-flat rectangular plate with one corner raised on data-dependent triangulation (minimal energy), C1 continuous surface.

curves of the surface with constant Y planes are computed at Y = 2.7292, 2.1833, 1.6375, 1.0917, 0.5458, and 0.0000, and these are shown in Fig. 15. The intersection curves of the surface with constant Z planes are computed at Z = 1.3708, 1.9667, 2.5625, 3.1583, 3.7542, and 3.8250, and these are shown in Fig. 16. And, since this example is from an industry, for the sake of better acceptance we plot the intersection curves (i.e. intersection curves of the surface with constant X planes: Body plan, constant Y planes: Profile plan, and constant Z planes: Waterline curves) as per conventions of naval architecture. As mentioned before, the surface is non-flat but is C1 continuous, so all the intersection curves are continuous without any break in continuity. Since the surface is a C1 continuous surface the intersection curves are continuous across the curves of C1 continuity. Furthermore, the surface is free from undesired oscillations and it is clear from Fig. 13a and 13b that the surface is fair surface.

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

2187

Fig. 12. Basic triangulation (F) over the given point set for the forward part of a hard-chine fishing vessel.

(a) top view

(b) isometric view Fig. 13. Surface for the forward part of a hard-chine fishing vessel on data-dependent triangulation (minimal energy), C1 continuous surface.

2188

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

Fig. 14. Intersection curves of constant X planes with surface for the forward part of a hard-chine fishing vessel on data-dependent triangulation (minimal energy), C1 continuous surface.

Fig. 15. Intersection curves of constant Y planes with surface for the forward part of a hard-chine fishing vessel on data-dependent triangulation (minimal energy), C1 continuous surface.

Fig. 16. Intersection curves of constant Z planes with surface for the forward part of a hard-chine fishing vessel on data-dependent triangulation (minimal energy), C1 continuous surface.

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

2189

6. Conclusions The present work has discussed a general problem, and explored the possibility of applying the construction of piecewise polynomial parametric surfaces which interpolate prescribed R3 scattered data using spaces of parametric splines defined on R3 triangulation to design fair surfaces over scattered data in R3 by extending, and generalizing the works of Quak and Schumaker [6], Fasshauer and Schumaker [8], and Quak, and Schumaker [40]. Since in this work the triangulation is defined on R3 no projection is required, hence the formulation is able to efficiently include the cases where the data points have a tendency to converge/collapse. We have explored two algorithms for interpolating the surface over scattered data, and the two algorithms are compared via shape interrogation of the resulting surface. Specifically, we have investigated how bending energy is related to the non-linearity in the surface’s intersection curves with constant X, Y, and Z planes. Since, surface’s planar cuts play a important role in designing engineering bodies that have favorable manufacturing or hydrodynamic related properties, the presented algorithms can be integrated with CAD/CAM design tools for engineering design. Additionally, since the resulting surface is a piecewise cubic C1 surface and the selection of the triangulation is based on the idea of minimizing the energy of the resulting surface, the designed fair surface has a physically meaningful property. Though, we have shown via theoretical analysis and implementation with two problems that the developed algorithms are stable and computationally efficient, but the complete complexity analysis of the algorithms with detailed comparative study with other similar techniques has not been attempted. In the presented algorithms the initial geometry definition as a triangulation can be used to have a local control but geometric features need to be investigated in detail to have wider applications in the areas of CAGD. In our formulation we have used canonical triangle which is simple, unique, and natural, and result into better computational efficiency. However, the optimal shape of canonical triangle which would certainly allow us to have better topological features (i.e. homomorphism, isomorphism, and equivalence, etc.) of the interpolated surface has not been investigated in the present work. Also, further improvements of the present surface interpolation techniques in the application areas of automation, anisotropic generation, adaptation and integration with numerical solution processes may be needed for development of efficient software tools. Our future work shall go in this direction, and currently this is under investigation. Trademarks and copyrights *Trademark

and copyright with Silicon Graphics Corporation, USA. AutoDesk Corporation, USA.

**Trademark

and copyright with

Acknowledgements This research was partially supported by the Defense Research and Development Organization (DRDO), MoD, GOI, India; via the grant-in-aid scheme number: ERIP/ER/0400237/M/01. The authors thank Prof. (Retd.) R.P. Gokarn for his assistance and comments. Appendix A Proof of Theorem. We consider the polynomial piece p ¼ sI T i . Then as per the coefficients shown in Fig. 1, the Be´rnstein–Be´zier representation of p is given by, pðu; v; wÞ ¼ c19 u3 þ 3  c16 u2 v þ 3  c17 u2 w þ 3  c5 uv2 þ 6  c13 uvw þ 3  c8 uw2 þ c1 v3 þ 3  c4 v2 w þ 3  c9 vw2 þ c2 w3 ; where (u, v, w) are the barycentric coordinates of a point (x, y, z) in triangle T1 with respect to its vertices P4, P1, P2. Now p takes on the value G1 at P1 if and only if c1 = G1. Considering the interpolation conditions involv-

2190

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

ing the gradients at P1, let Dr denote the directional derivative in the direction r = (x2  x1, y2  y1, z2  z1). Then using the relation between the directional derivative of a polynomial and its Be´rnstein–Be´zier coefficients [44], it follows that 2 3 ðc4 ð3  ðc1 ÞÞÞ 6 7 Dr p ¼ 4qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi5: ðx2  x1 Þ2 þ ðy 2  y 1 Þ2 þ ðz2  z1 Þ2

ð5Þ

In terms of the gradient, the directional derivative can be written as, 2 3 y x z 6ðx2  x1 Þ  G1 þ ðy 2  y 1 Þ  G1 þ ðz2  z1 Þ  G1 7 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Dr p ¼ 4 5: 2 2 2 ðx2  x1 Þ þ ðy 2  y 1 Þ þ ðz2  z1 Þ

ð6Þ

At P1, equating Eqs. (5) and (6) we get,   ðx2  x1 Þ  Gx1 þ ðy 2  y 1 Þ  Gy1 þ ðz2  z1 Þ  Gz1 c4 ¼ ; 3 þ G1 which is as given in Theorem 1. Next we consider the derivative at normal to the edge e1 and pointing into the triangle, evaluating at a point (0, 1  v, v), it can be given by, 2 0 0 13 1 2 2 c5  ð1  vÞ þ 2  c13  vð1  vÞ h1 ð1  vÞ þ 2  ðh1  1Þ 63 @ A7 A þ 3  c4  @ 7 6 2 2 7 6 v  ð1  vÞ þc8  v þ c1  ð1  vÞ  ðh1  1Þ 7 6 7 6 ! 7 6 2:h1  v  ð1  vÞ þ ðh1  1Þ 7 6 5 4 þ3  c  9 2 2 v  h1  v  c2 Dr pð0; 1  v; vÞ ¼ ; d P 4 e1 ð7Þ where d P 4 e1 is the distance of the point P4 from the edge e1 and h1 is as defined in Theorem 1. Since we are discussing C1 continuous surface the polynomial in Eq. (7) will be consisting only first degree terms. This reduces to a first degree polynomial in v if and only if the second derivative of the polynomial in Eq. (7) with respect to v is equal to zero, i.e., 0 ¼ 2  c5  4  c13 þ 2  c8 þ 2  c1 ðh1  1Þ þ c4  ð2  h1  4  ðh1  1ÞÞ þ c9  ð4  h1 þ 2  ðh1  1ÞÞ  2  h1  c2 : ð8Þ

Rearranging the terms in Eq. (8), we get,   ðc5 þ c8 þ ðh1  1Þ  c1 Þ þ ð2  3h1 Þ  c4 þ ð3h1  1Þ  c9  h1  c2 c13 ¼ ; 2 which is as given in Theorem 1. Similarly, considering the triangles T2, T3 we can get the formulas for the coefficients c6, c7, c10, c11, c12, c14 and c13. Considering C1 continuity across the interior edges, since the barycentric coordinate of the point P4 in terms of P1, P2 and P3 is ð13 ; 13 ; 13Þ, being the ‘center of gravity’ of the triangle having vertices at P1, P2 and P3, let us consider the triangle whose vertices are on c5, c13 and c15, then c16 will be at its center of gravity being the Clough–Tocher split. It will lead to, c16 ¼

hc þ c þ c i 15 5 13 : 3

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

2191

Similarly, hc þ c þ c i 13 8 14 for triangle having vertices at c8 ; c13 and c14 ; and; 3 hc þ c þ c i 14 11 15 c18 ¼ for triangle having vertices at c11 ; c14 and c15 : 3

c17 ¼

At last from any one of the C1 continuity conditions across an edge we get, hc þ c þ c i 18 16 17 c19 ¼ for triangle having vertices at c16 ; c17 and c18 : 3 Appendix B. Computation of energy of a piecewise triangular polynomial patch Because of symmetries in Eq. (4) only 13 entries are sufficient. Since we are considering a rotation invariant energy expression, let the triangle Ti be in a position as shown in Fig. 17. A canonical triangle is a triangle that has been reduced to the simplest and most significant form possible without loss of generality, and it has a ‘uniqueness’ or ‘naturalness’ inherent in its formulation. Furthermore, it is independent of coordinates. And, because of these qualities canonical triangles are widely used in computer aided geometric design, computational geometry, and graphics, e.g. [45–53]. Additionally, it can be proved that the canonical triangles of optimally-efficient shape have vertices in or on the unit circle and also on some interpolation conic whose size and shape are determined by the interpolation error and the data function through which interpolation is being done. The canonical triangle shape is identified by translating the common center to the origin and rescaling the gradient error circle to the unit circle circumscribing either equilateral triangle (all three faces equal) or isosceles triangle (two faces equal). Though we have not utilized this property in our formulation but certainly its implementation will lead to better topological features of the interpolated surface. We can note here that the computation for canonical triangle is for computational efficiency, and any general triangle can be transformed into a canonical triangle with some definitive transformations. These transformations are given in Appendix C. We may suppose that the triangle is in canonical position shown in Fig. 17 with vertices P1 = (b, c, d), P2 = (0, 0, 0), and P3 = (a, 0, 0). Now the Barycentric coordinates of any point (x, y, z) are given by (u, v, w), where, sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi y 2 þ z2 ; u¼ c2 þ d 2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 ðc  z  d  yÞ þ ðdðx  aÞ þ zða  bÞÞ þ ðyðb  aÞ þ cða  xÞÞ v¼ ; a2  ðc2 þ d 2 Þ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðc  z  d  yÞ2 þ ðd  x þ b  zÞ2 þ ðb  y þ c  xÞ2 : and w ¼ a2  ðc2 þ d 2 Þ

Y -axis

P1 (b ,c , d )

P2 (0 ,0 ,0 )

P3 (a ,0 ,0 )

Fig. 17. The canonical triangle.

X -axis

2192

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

Then the gradients are,

! !! 1 1 u ¼ ðux ; uy ; uz Þ ¼ 0; pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ;  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; c2 þ d 2 c2 þ d 2 0 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 2 2 1 1 ðb  aÞ þ d 2 1 c2 þ ða  bÞ A v ¼ ðvx ; vy ; vz Þ ¼ @ ; ; ; a a a c2 þ d 2 c2 þ d 2 0 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 1 1 b2 þ d 2 1 b2 þ d 2 A ; : and w ¼ ðwx ; wy ; wz Þ ¼ @ ;  a a c2 þ d 2 a c2 þ d 2 E

ð3Þ

1 The 13 essential entries (for details see [6,8], and [40,41]) of GS ¼ 5 9ac are then given by,

2

GS 1;1 ¼ ðU  U Þ ; GS 1;2 ¼ 2  ðU  U ÞðU  V Þ þ 0:5ðU  U Þ2 ; 2

GS 1;4 ¼ ðU  V Þ þ ðU  U ÞðU  V Þ; GS 1;5 ¼ 2  ðU  V ÞðU  W Þ þ ðU  U ÞðU  W Þ þ ðU  U Þ2 ; 2

GS 1;7 ¼ 0:5ðU  V Þ ; 2

GS 1;8 ¼ ðU  V ÞðU  W Þ þ 0  5ðU  V Þ ; 2

2

GS 2;2 ¼ 2  ðU  U ÞðV  V Þ þ 2ðU  V Þ þ 2ðU  U ÞðU  V Þ þ ðU  U Þ ; 2

GS 2;3 ¼ 2  ðU  U ÞðV  W Þ þ 2ðU  V ÞðU  W Þ þ ðU  U ÞðU  W Þ þ ðU  U ÞðU  V Þ þ 0  5ðU  U Þ ; GS 2;4 ¼ 2  ðV  V ÞðV  U Þ þ ðU  U ÞðV  V Þ þ 1:5ðU  U Þ2 þ 2ðU  U ÞðU  V Þ; GS 2;5 ¼ 2  ðV  V ÞðU  W Þ þ 2ðU  V ÞðV  W Þ þ ðU  U ÞðV  W Þ þ 2ðU  V ÞðU  W Þ þ ðU  U ÞðV  V Þ 2

þ ðU  V Þ þ 2ðU  U ÞðU  W Þ þ ðU  U ÞðU  V Þ; GS 2;6 ¼ 2  ðU  W ÞðV  W Þ þ 0  5ðU  W Þ2 þ ðU  U ÞðV  W Þ þ ðU  V ÞðU  W Þ þ ðU  U ÞðU  W Þ; 2

GS 2;9 ¼ ðU  W ÞðV  W Þ þ ðV  V ÞðU  W Þ þ ðU  V ÞðV  W Þ þ ðU  W Þ þ ðU  V ÞðU  W Þ; 2

GS 5;5 ¼ 2  ½ðV  V ÞðW  W Þ þ ðV  W Þ þ ðW  W ÞðU  W Þ þ ðU  W ÞðV  W Þ þ ðV  V ÞðU  W Þ 2

þ ðU  V ÞðV  W Þ þ ðU  U ÞðW  W Þ þ ðU  W Þ þ ðU  U ÞðV  W Þ þ ðU  V ÞðU  W Þ 2

þ ðU  U ÞðV  U Þ þ ðU  V Þ ; where, U ¼ ux i þ uy j þ uz k;

V ¼ vx i þ vy j þ vz k;

and

W ¼ wx i þ wy j þ wz k;

and since the matrix is symmetrical so, GSij = GSji. In the above expressions U, V, W are vectors and the expressions U Æ U, U Æ V etc., represent their dot products. Thus, for example, ! !!  2 2 2 2 2 1 ðb  aÞ þ d c þ ða  bÞ V  V ¼ ðv2x þ v2y þ v2z Þ ¼ ; ; 2 : a a2  ðc2 þ d 2 Þ a  ðc2 þ d 2 Þ Appendix C. Basic transformations Let the coordinates of a general triangle be (x1, y1, z1), (x2, y2, z2) and (x3, y3, z3) as shown in Fig. 18. This can be transformed to rotation invariant canonical position triangle as shown in Fig. 17 by performing the following transformations:

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

2193

Yaxis

P1 (x1 , y2 , z1 )

P2 (x2 , y2 , z2 )

X-axis

P3 (x3 , y3 , z3 )

Zaxis

Fig. 18. A general triangle.

I. Transferring the origin to any of the vertices of the triangle. The new coordinates are given by, f x

y

z

1g ¼ fx y

z

1 g  fT F 1 g;

where {TF1} is the first transformation matrix. Suppose 9 we transfer the origin to (x1, y1, z1), then the 8 1 0 0 0 > > > > = < 0 1 0 0 . {TF1} can be written as fT F 1 g ¼ 0 1 0> > 0 > > ; : x1 y 1 z1 1 II. Rotating in the X  Y plane by angle h (-clock wise about Z-axis). Since the rotation is in X  Y plane z will remain constant and, ðx2  x1 Þ cos h ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; 2 2 ðx2  x1 Þ þ ðy 2  y 1 Þ

ðy 2  y 1 Þ sin h ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 2 2 ðx2  x1 Þ þ ðy 2  y 1 Þ

The second transformation matrix can be written as 9 8 ðx2 x1 Þ ðy 2 y 1 Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 0 0> > 2 2 2 2 > > ðx x Þ þðy y Þ ðx x Þ þðy y Þ > > 2 1 2 1 2 1 2 1 > > > > = < ðy 2 y 1 Þ ðx2 x1 Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 0 0 2 2 2 2 : fT F 2 g ¼ ðx2 x1 Þ þðy 2 y 1 Þ ðx2 x1 Þ þðy 2 y 1 Þ > > > > > > 0 0 1 0 > > > > ; : 0 0 0 1 III. Rotating in the X–Z plane by angle / (-clock wise about Y-axis). Since the rotation is in X–Z plane y will remain constant and qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ðx2  x1 Þ þ ðy 2  y 1 Þ cos / ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; and 2 2 2 ðx2  x1 Þ þ ðy 2  y 1 Þ þ ðz2  z1 Þ ðz2  z1 Þ sin / ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 2 ðx2  x1 Þ þ ðy 2  y 1 Þ2 þ ðz2  z1 Þ2 The third transformation matrix can be written as,

2194

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

ffi 8 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx2 x1 Þ2 þðy 2 y 1 Þ2 ðz2 z1 Þ > p pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > > ðx2 x1 Þ2 þðy 2 y 1 Þ2 þðz2 z1 Þ ðx2 x1 Þ2 þðy 2 y 1 Þ2 þðz2 z1 Þ > > < 0 0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi fT F 3 g ¼ 2 2 > ðx x Þ ðz z Þ 2 1 þðy 2 y 1 Þ 2 1 > p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > > ðx2 x1 Þ2 þðy 2 y 1 Þ2 þðz2 z1 Þ ðx2 x1 Þ2 þðy 2 y 1 Þ2 þðz2 z1 Þ > : 0 0

0 0 1 0

9 > 0> > > > > = 0 : > 0> > > > > ; 1

Finally the transformation matrix will be as, fT F g ¼ fT F 1 g  fT F 2 g  fT F 3 g

and

f x

y

z

1g ¼ fx

y

z

1 g  fT F g:

References [1] L.L. Schumaker, Fitting surfaces to scattered data, in: G.G. Lorentz, C.K. Chui, L.L. Schumaker (Eds.), Approximation Theory II, Academic Press, New York, USA, 1976, pp. 203–268. [2] R. Franke, Recent advances in the approximation of surfaces from scattered data, in: C.K. Chui, L.L. Schumaker, F. Uterras (Eds.), Topics in Multivariate Approximation, Academic Press, New York, 1987, pp. 79–98. [3] S.K. Lodha, R. Franke, Scattered data techniques for surfaces, in: H. Hagen, G. Nielson, F. Post (Eds.), Proceedings of Dagstuhl Conference on Scientific Visualization, IEEE Computer Society Press, 1999, pp. 182–222. [4] N. Dyn, D. Levin, S. Rippa, Algorithms for the construction of data dependent triangulation, in: M.G. Cox, J.C. Mason (Eds.), Algorithms for Approximation II, Clarendon Press, Oxford, UK, 1989, pp. 21–31. [5] N. Dyn, D. Levin, S. Rippa, Data dependent triangulations for piecewise linear interpolation, IMA J. Numer. Anal. 10 (1) (1990) 137–154. [6] E. Quak, L.L. Schumaker, Cubic spline fitting using data dependent triangulations, Comput. Aided Geom. Des. 7 (1990) 293–301. [7] R.F. Sarraga, Recent methods for surface shape optimization, Comput. Aided Geom. Des. 15 (1998) 417–436. [8] E.G. Fasshauer, L.L. Schumaker, Minimal energy surfaces using parametric splines, Comput. Aided Geom. Des. 13 (1996) 45–79. [9] H. Qin, D. Terzopoulos, D-NURBS: a physics-based geometric design framework, IEEE Trans. Visualiz. Comput. Graph. 2 (1) (1996) 85–96. [10] D. Terzopoulos, H. Qin, Dynamic NURBS with geometric constraints for interactive sculpting, ACM Trans. Graph. 13 (2) (1994) 103–136. [11] R. Schmidt, Eine Methode zur Konstruktion von C1-Fla¨chen zur Interpolation unregelma¨ssig vereilter Daten, in: W. Schempp, K. Zeller (Eds.), Multivariate Approximation Theory II, Birltha¨user, Basel, 1982, pp. 343–361. [12] G.M. Nielson, J. Tvedt, Comparing methods of interpolation for scattered volumetric data, in: D. Rogers, R.A. Earnshaw (Eds.), State of the Art in Computer Graphics – Aspects of Visualization, Springer-Verlag, 1994, pp. 67–86. [13] G.M. Nielson, J. Tvedt, Modeling of scattered multivariate data, in: C. Giertson, P. Fevang (Eds.), Euro-graphics’94 State of the Art Reports, Norsied Press, Oslo, Norway, 1994, pp. 38–59. [14] G.M. Nielson, Scattered data modeling, IEEE Comput. Graph. Appl. 13 (1) (1993) 60–70. [15] G.M. Nielson, T. Dierks, Modelling and visualization of scattered volumetric data, in: E.J. Farrell (Ed.), Extracting Meaning from Complex Data: Processing, Display, Interaction II, Springer, Verlag, 1991, pp. 22–33. [16] J. Tvedt, A Software System for Comparison of Scattered Data Interpolation Methods, M.S. Thesis, Arizona State University, Tempe, USA, 1991. [17] P. Alfeld, Scattered data interpolation in three or more variables, in: T. Lyche, L.L. Schumaker (Eds.), Mathematical Methods in CAGD, 1989, pp. 1–33. [18] T.A. Foley, H. Hagen, Advances in scattered data interpolation, Surveys Math. Ind. 4 (1994) 71–84. [19] G.M. Nielson, T. Foley, Modeling of scattered multivariate data, in: C. Giertsen, P. Fevang (Eds.), Euro-Graphics’94 State of the Art Reports, Norsied Press, Oslo, Norway, 1994, pp. 38–59. [20] G.M. Nielson, A characterization of an affine invariant triangulation, in: G. Farin, H. Hagen, H. Noltemeier, W. Knoedel (Eds.), Geometric Modelling: Computing Supplementum, Springer-Verlag, 1993, pp. 191–210. [21] S. Mann, C. Loop, M. Lounsbery, D. Meyers, J. Painter, T. DeRose, K. Sloan, A survey of parametric scattered data fitting using triangular interpolants, in: H. Hagen (Ed.), Curve and Surface Design, SIAM, 1992, pp. 145–172. [22] R.E. Barnhill, G.M. Nielson, Introduction to surfaces, Rocky Mountain J. Math. 14 (1984) 1–4. [23] R. Franke, Scattered data interpolation: tests of some methods, Math. Comput. 38 (1982) 181–200. [24] L.L. Schumaker, Fitting surfaces to scattered data, in: Approximation Theory II, Academic Press, New York, 1976, pp. 203–268. [25] C.L. Bajaj, Surface fitting with implicit algebraic surface patches, in: H. Hagen (Ed.), Topics in Surface Modeling, SIAM, Philadelphia, 1992, pp. 23–52. [26] M.D. Buhmann, New developments in the theory of radial basis function interpolation, in: K. Jetter, F.I. Utreras (Eds.), Multivariate Approximation: CAGD to Wavelets, World-Scientific, Singapore, 1993, pp. 35–76. [27] R. Hardy, Theory and applications of the multiquadric–biharmonic method, Comput. Math. Appl. 19 (1990) 163–208.

R. Sharma, O.P. Sha / Applied Mathematical Modelling 32 (2008) 2172–2195

2195

[28] M.J.D. Powell, The theory of radial basis function approximation, in: W. Light (Ed.), Advances in Numerical Analysis (Volume II): Wavelets, Subdivision Algorithms, and Radial Basis Functions, Clarendon Press, Oxford, 1990, pp. 105–210. [29] G.M. Nielson, Research issues in modeling for the analysis and visualization of large data sets, in: Scientific Visualization: Advances and Challenges, Academic Press, USA, 1994, pp. 143–157. [30] G.M. Nielson, CAGD’s top ten: what to watch? IEEE Comput. Graph. Appl. 13 (1) (1993) 35–37. [31] N. Dyn, Subdivision schemes in computer-aided geometric design, in: W. Light (Ed.), Advances in Numerical Analysis (Volume II): Wavelets, Subdivision Algorithms, and Radial Basis Functions, Clarendon Press, Oxford, 1992, pp. 36–104. [32] F. Preparata, M.I. Shamos, Computational Geometry: An Introduction, Springer Verlag, 1985. [33] J.L. Brown, Vertex based data dependent triangulations, Comput. Aided Geom. Des. 8 (1991) 239–251. [34] N. Dyn, D. Levin, S. Rippa, Numerical procedures for surface fitting of scattered data by radial functions, SIAM J. Scient. Statist. Comput. 7 (1986) 639–659. [35] R. Clough, J. Tocher, Finite element stiffness matrices for analysis of plates in bending, in: Proceedings of Conference on Matrix Methods in Structural Mechanics, Air Force Institute of Technology, Ohio, USA, 1965. [36] R.J. Renka, Triangulation and interpolation at arbitrarily distributed points in the plane Algorithm 624, ACM Trans. Math. Software 10 (1984) 440–442. [37] C. Lawson, Software for C1 surface interpolation, in: J.R. Rice (Ed.), Mathematical Software III, Academic Press, New York, 1977, pp. 161–194. [38] J.N. Reddy, Theory and Analysis of Elastic Plates, Taylor and Francis Publishing Limited, UK, 1999. [39] S.P. Timoshenko, S. Woinowsky-Krieger, Theory of Plates and Shells, second ed., McGraw Hill College Division, New York, 1959, pp. 62–67. [40] E. Quak, L.L. Schumaker, Calculation of the energy of a piecewise polynomial surfaces, in: M.G. Cox, J.C. Mason (Eds.), Algorithms for Approximation II, Clarendon Press, Oxford, 1989, pp. 23–31. [41] O.P. Sha, R. Sharma, Hull surface modeling with triangular patches, Design Laboratory Memorandum 17 – 1999, Department of Ocean Engineering and Naval Architecture, Indian Institute of Technology Kharagpur, India, 1999. [42] R. Sharma, O.P. Sha, A tracing method for parametric Be´zier triangular surface/plane intersection, Int. J. Comput. Appl. Technol. 28 (4) (2007) 240–253. [43] J.S. Chalfant, T. Maekawa, Design for manufacturing using B-spline developable surfaces, J. Ship Res. 42 (3) (1998) 207–215. [44] G. Farin, Triangular Bernstein–Be´zier patches, Comput. Aided Geom. Des. 3 (1987) 83–127. [45] R.L. Goldstone, A. Jones, M.E. Roberts, Group path formation, IEEE Trans. Syst. Man Cybernet. – Part A: Syst. Humans 36 (3) (2006) 611–620. [46] M. Sharir, H. Shaul, Ray shooting and stone throwing with near-linear storage, Comput. Geom. 30 (2005) 239–252. [47] W.-S. Chun, T.J. Purtell, O.S. Cossairt, Technical Introduction to the Basic Spatial GL Rendering Pipeline for a Perspecta 1.9/2.0, Public Information Disclosure, ASI-0110, Actuality Systems, USA, 2005, pp. 1–15. [48] J. Matousˇek, On constants for cuttings in the plane, J. Discrete Comput. Geom. 20 (4) (2004) 427–448. [49] N. Wu, Hugues Libergier and his instruments, Nexus, 00/01_017-102, 2001, pp. 93–102. [50] P. Agarwal, P. Desikan, An efficient algorithm for terrain simplification, in: Proceedings of 8th ACM–SIAM Symposium on Discrete Algorithms, 1997, pp. 23–32. [51] I.J. Dejter, TMC tetrahedral types MOD 2k + 1 and their structure graphs, Graphs Combin. 12 (1996) 163–178. [52] P.K. Agarwal, M.J. Katz, M. Sharir, Computing depth orders for fat objects and related problems, Comput. Geom. 5 (1995) 187–206. [53] E.F. D’Azevedo, R.B. Simpson, On optimal triangular meshes for minimizing the gradient error, Numer. Math. 59 (1991) 321–348.