Quasi-interpolation by quadratic C1 -splines on truncated octahedral partitions

Computer Aided Geometric Design 26 (2009) 825–841

Contents lists available at ScienceDirect

Computer Aided Geometric Design www.elsevier.com/locate/cagd

Quasi-interpolation by quadratic C 1 -splines on truncated octahedral partitions M. Rhein a,∗ , T. Kalbe b a b

University of Mannheim, Institute for Mathematics, 68131 Mannheim, Germany TU Darmstadt, GRIS, 64283 Darmstadt, Germany

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 5 November 2008 Accepted 22 April 2009 Available online 3 May 2009

We describe an approximating scheme for the smooth reconstruction of discrete data on volumetric grids. A local quasi-interpolation method for quadratic C 1 -splines on uniform tetrahedral partitions is used to achieve a globally smooth function. The Bernstein–Bézier coefficients of the piecewise polynomials are thereby directly determined by appropriate combinations of the data values. We explicitly give a construction scheme for a family of quasi-interpolation operators and prove that the splines and their derivatives can provide an approximation order two for smooth functions. The optimal approximation of the derivatives and the simple averaging rules for the coefficients recommend this method for high quality visualization of volume data. Numerical tests confirm the approximation properties and show the efficient computation of the splines. © 2009 Elsevier B.V. All rights reserved.

MSC: 41A15 41A25 65D07 65D17 65D18 Keywords: Piecewise quadratic polynomials Trivariate splines Visualization Gridded volume data Approximation order

1. Introduction The construction of adequate non-discrete models from discrete data on volumetric grids is a major task in a multitude of applications, e.g. medical imaging, industrial quality control, scientific visualization or numerical simulation. Modern scanning devices, such as CT and MRI allow to measure huge volume data sets. When dealing with such data, one seeks for finding efficient (real-time) reconstruction schemes suitable for fast and artifact-free visualization of the data. In numerical simulations, it is of importance for the underlying model to have guaranteed approximation properties. The most common approach in the above applications is trilinear interpolation (see Bajaj, 1999; Meissner et al., 2000; Theußl et al., 2002), which accounts for a tensor–product extension of univariate linear splines interpolating at the grid points, leading to piecewise cubic polynomials. Besides the simplicity of this local model, one of the main reasons for its widespread use is the fact that a sufficiently smooth function is approximated with order two. On the other hand, derivatives are not directly available from this model. In addition, the reconstruction is not smooth, which leads to several problems, e.g. visualization artifacts or stability issues. Tri-quadratic and tri-cubic tensor–product splines can be used to construct smooth models, but lead to piecewise polynomials of relatively high degree, namely six and nine.

*

Corresponding author. E-mail addresses: [email protected] (M. Rhein), [email protected] (T. Kalbe).

0167-8396/$ – see front matter doi:10.1016/j.cagd.2009.04.002

© 2009

Elsevier B.V. All rights reserved.

826

M. Rhein, T. Kalbe / Computer Aided Geometric Design 26 (2009) 825–841

Trivariate C 1 -splines on tetrahedral partitions, where the piecewise polynomials are given in their Bernstein–Bézier form are well suited for the purpose of volume approximation and visualization, since the Bernstein–Bézier techniques can be fully employed for the efficient computation and evaluation of the splines (see de Boor, 1987; Farin, 1986). However, these are very complex spaces and only a few results are known up to date, see Lai and Schumaker (2007) and the references therein. Furthermore, C 1 -splines with good approximation properties, that solve Hermite- and Lagrange-Interpolation problems must have a relative high degree (at least 5), see Lai and Le Méhauté (2004), Schumaker and Sorokina (2004) and Ženišek (1973), or make use of macro-elements, see Schumaker et al. (2009), Sorokina and Worsey (2008), Worsey and Farin (1987), Worsey and Piper (1988), where the tetrahedra of the partition are further subdivided, which results in complex partitions with a huge number of tetrahedra. Compared to that, quasi-interpolation (see de Boor and Fix, 1973) is a good approach which allows the usage of low degree splines on uniform partitions, while still some approximation properties can be guaranteed. The recent work of Sorokina and Zeilfelder (2007) is most comparable to our approach. Therefore, our paper follows closely their structural outline and methodology. They avoid most of the above mentioned problems by constructing a C 1 -smooth, cubic spline model based on a type-6 tetrahedral partition of the domain, which means that each data cube is subdivided into 24 disjoint tetrahedra. The splines are built from appropriate averaging formulas in a symmetric and local neighborhood of the data values in a way such that the resulting splines approximate the values and the first derivatives of a sufficiently smooth function simultaneously with order two. Previous works (see Hangelbroek et al., 2004 and Nürnberger et al., 2005) have shown that these results cannot be significantly improved. Particularly, the degree of the splines cannot be further decreased without the loss of either the approximation properties or the C 1 -smoothness of the resulting splines. In our work, we describe the first C 1 -smooth quadratic spline model on volumetric grids. Thereby, we employ a different uniform tetrahedral partition, which is obtained by subdividing each truncated octahedron of a truncated octahedral partition into 144 disjoint tetrahedra. This split is more complex than the 24-split of the cubes in the type-6 partition. On the other hand, the C 1 -splines on truncated octahedral partitions possess significantly more degrees of freedom than C 1 -splines on cubic partitions. For the approximation of the same volumetric data set, our spline method requires only approximately 50% more tetrahedra, while allowing for the lowest possible spline degree. The low degree of the polynomial patches results in a more efficient evaluation of the spline pieces, e.g., for intersecting rays with spline patches and the calculation of derivatives. We are able to describe an entire family of quasi-interpolating C 1 -splines on truncated octahedral partitions and prove that for a specific spline the same approximation properties are achieved as for the cubic splines proposed in Sorokina and Zeilfelder (2007). This means that our quasi-interpolating splines are able to approximate the derivatives of sufficiently smooth functions with optimal approximation order. The paper is organized as follows. In Section 2, we recall important facts on trivariate splines. Then, in Section 3, we introduce truncated octahedral partitions and analyze the structure of C 1 -splines on these partitions. In Section 4, we describe our quasi-interpolation scheme for the construction of a family of C 1 -splines and prove some important results for certain quasi-interpolation operators. Section 5 contains the approximation properties of our method and in Section 6 we confirm these theoretical results with numerical tests. Finally, Section 7 concludes the paper with comments on our future research and additional remarks. 2. Preliminaries Let  be a tetrahedral partition of a set Ω in R3 . In the following, we consider the space of quadratic C 1 -splines on , defined by

  S21 () = s ∈ C 1 (Ω): s| T ∈ P2 , for all T ∈  ,

where C 1 (Ω) is the set of continuously differentiable functions on Ω , and P2 denotes the 10-dimensional space of trivariate polynomials of degree two. We use the well-known Bernstein–Bézier form (B-form) of the polynomial pieces (see de Boor, 1987; Farin, 1986; Lai and Schumaker, 2007). Given a tetrahedron T := [ v 0 , v 1 , v 2 , v 3 ] ∈  with vertices v 0 , v 1 , v 2 , and v 3 , each p ∈ P2 has a unique representation



p≡

b i jk B iTjk ,

(1)

i + j +k+=2

where B iTjk ≡

2 i ! j !k!!

j

φ0i φ1 φ2k φ3 ,

i + j + k +  = 2,

are the quadratic Bernstein polynomials on T . Here, the functions φ0 , φ1 , φ2 , and φ3 denote the barycentric coordinates with respect to T , which are linear trivariate polynomials determined by

φi ( v j ) = δi , j ,

j = 0, . . . , 3,

where δi , j is the Kronecker’s symbol. The Bernstein polynomials form the partition of unity, this means



i + j +k+=2

B iTjk ≡ 1.

(2)

M. Rhein, T. Kalbe / Computer Aided Geometric Design 26 (2009) 825–841

827

As usual, we associate the Bernstein–Bézier coefficients (B-coefficients) b i jk of p in the form (1) with the domain points

ξi jk := (i v 0 + jv 1 + kv 2 +  v 3 )/2,

i + j + k +  = 2,



in T , and we call the points (ξ, bξ ) ∈ R4 control points of p. Let D = T ∈ D T be the union of the sets of domain points associated with the tetrahedra of . It is easy to see that a continuous spline s ∈ S20 () is uniquely defined by the coefficients {bξ : ξ ∈ D } and the dimension of S20 () is equal to the cardinality of D . Therefore, it follows that dim S20 () = #D = V + E , where V and E denote the vertices and edges of . If s ∈ C 1 (Ω), then these coefficients cannot be chosen arbitrarily, but have to satisfy certain smoothness conditions. Let T = [ v 0 , v 1 , v 2 , v 3 ] and  T = [ v 0 , v 1 , v 2 , v˜ 3 ] be two adjacent tetrahedra sharing a common triangular face F := T ∩  T = [ v 0 , v 1 , v 2 ] and let s be a quadratic continuous spline on T ∪  T in its ˜ i jk , respectively. Then s is C 1 ˜ = p b and b with the corresponding coefficients piecewise B-form (1): s| T = p and s| i jk  T smooth across F if and only if b˜ i jk1 = b i +1, j ,k,0 φ0 ( v˜ 3 ) + b i , j +1,k,0 φ1 ( v˜ 3 ) + b i , j ,k+1,0 φ2 ( v˜ 3 ) + b i , j ,k,1 φ3 ( v˜ 3 ),

(3)

where i + j + k = 1. In general, there are five coefficients involved in each of these three conditions. But these trivariate smoothness conditions can degenerate to lower dimensional conditions which are similar to the bivariate and univariate setting, if one or two of the barycentric coordinates vanish at the point v˜ 3 . The appearance of these degenerated conditions is typical for splines on uniform partitions (see Hangelbroek et al., 2004; Hecklin et al., 2007). Especially spline spaces with low polynomial degree are highly complex spaces, due to the fact that the smoothness conditions (3) have to be simultaneously satisfied for all interior faces of . We conclude this section with two further definitions. Let v ∈ V be a vertex of . We define by D 1 ( v ) the ball with radius 1 around v. Note that D 1 ( v ) corresponds to the set of midpoints emanating from v and v itself. An easy geometric interpretation of (3) shows that a spline s ∈ S20 () is C 1 -smooth in a vertex v, iff all control points (ξ, bξ ), ξ ∈ D 1 ( v ) lie in the same hyperplane in R4 . Let M ⊆ D and bξ , ξ ∈ M, be the associated B-coefficients of a spline s ∈ S21 (). Then, we call M a determining set for the space S21 (), if bξ = 0, ξ ∈ M, implies s ≡ 0. In the next section we will use these results to analyze quadratic C 1 -splines on a special type of uniform tetrahedral partitions, namely quadratic C 1 -splines on truncated octahedral partitions. 3. Quadratic C 1 -splines on truncated octahedral partitions Let VD be the Voronoi diagram of the Body Centered Cubic Lattice (BCC-Lattice). The BCC-Lattice is obtained from a regular cubic grid by inserting additional grid points at the center of each of the grid cells. The Voronoi cells of VD are truncated octahedra and we note that it is possible to tessellate the R3 by a uniform partition of truncated octahedra centered at the vertices of the BCC-Lattice. A truncated octahedron is an Archimedean solid with 14 faces consisting of 6 squares and 8 regular hexagons, see Fig. 1. The 24 vertices of a truncated octahedron T with radius h ∈ N centered at the origin can be described by all permutations of (0, ± h2 , ±h). Let 3 be a truncated octahedral partition and Ω ⊂ R3 the covered region in R3 . We get a tetrahedral partition  of Ω by subdividing all T ∈ 3 in the following way. We insert a vertex at the center of all faces and edges of T , as well as at the center of T itself. Then we triangulate the faces of T by connecting all vertices on the boundary of the faces with the vertex in the center of the face. Finally we connect all vertices on the boundary of T with the vertex in the center of T and achieve a tetrahedral partition , where every truncated octahedron T ∈ 3 is split into 144 tetrahedra. Each

Fig. 1. A truncated octahedron consists of 6 square faces and 8 hexagonal faces.

828

M. Rhein, T. Kalbe / Computer Aided Geometric Design 26 (2009) 825–841

tetrahedron T = [ v 0 , v 1 , v 2 , v 3 ] ∈  has one vertex, say v 0 at the center of a truncated octahedron T ∈ 3, one vertex, say v 1 at the center of a face of T , another vertex, say v 2 at the midpoint of one of the edges of T and the remaining vertex v 3 is a vertex of T . Thereby, the triangular face [ v 1 , v 2 , v 3 ] lies on a face of T . Since the faces of T are either squares or hexagons this split results in two differently shaped tetrahedra, according to the face they are associated with. We say T ∈ T S if [ v 1 , v 2 , v 3 ] lies on a square face of T and T ∈ T H if [ v 1 , v 2 , v 3 ] lies on a hexagonal face of T . In the following we want to describe the C 1 -smoothness conditions of a spline s ∈ S21 () across the interior triangular faces of . Thereby we use the notation F ν = [ v μ1 , v μ2 , v μ3 ], ν = 0, . . . , 3, μ1 , μ2 , μ3 ∈ {0, . . . , 3}\{ν }, μ1 < μ2 < μ3 . Here, F ν , ν = 0, . . . , 3, is the triangular face of T opposing the vertex v ν . Further, we denote the triangular faces of a tetrahedron T ∈  by F νS , ν = 0, . . . , 3, if T ∈ T S and by F νH , ν = 0, . . . , 3, otherwise. Note that the triangular face F 1S is the intersection between two adjacent tetrahedra of T S and T H , while the triangular face F 1H is the intersection between two adjacent tetrahedra of T H only. Let T and  T be two tetrahedra sharing a triangular face and b i jk and b˜ i jk the corresponding coefficients in the ˜ , respectively. In the case T ∩  T = F 1S we identify the representation (1) of the piecewise polynomials s| T = p and s| T = p

T ∈ T H , while in all other cases the conditions are completely symmetric. By using (3) and some tetrahedra by T ∈ T S and  elementary computations, we obtain that s ∈ S21 () iff the following conditions for its coefficients in the representation (1) are satisfied:

• Smoothness across F 0S and F 0H b1100 = 2b0200 − b˜ 1100 ,

(4a)

b1010 = 2b0110 − b˜ 1010 ,

(4b)

b1001 = 2b0101 − b˜ 1001 .

• Smoothness across

F 3S

and

(4c) F 3H

b1001 = 2b1010 − b˜ 1001 ,

(5a)

= 2b0110 − b˜ 0101 ,

(5b)

b0011 = 2b0020 − b˜ 0011 .

(5c)

b0101

• Smoothness across

F 2S

b1010 = b1100 − b˜ 1010 + b1001 ,

(6a)

b0110 = b0200 − b˜ 0110 + b0101 ,

(6b)

b0011 = b0101 − b˜ 0011 + b0002 .

(6c)

• Smoothness across b1010 = b0110 = b0011 =

1 2 1 2 1 2

F 2H

b1100 − b˜ 1010 + b0200 − b˜ 0110 + b0101 − b˜ 0011 +

3 2 3

b1001 ,

(7a)

b0101 ,

(7b)

b0002 .

(7c)

2 3 2

• Smoothness across F 1S 2˜ 4 b1100 + b1010 , 3 3 2˜ 4 b0110 = b1010 − b0110 + b0020 , 3 3 3 1 2˜ 4 b0101 = b1001 − b0101 + b0011 . 3 3 3 b1100 =

1

3 1

b2000 −

(8a) (8b) (8c)

• Smoothness across F 1H b1100 =

2

b2000 − b˜ 1100 +

4

b1010 , 3 4 b0110 = b1010 − b˜ 0110 + b0020 , 3 3 2 4 b0101 = b1001 − b˜ 0101 + b0011 . 3 3 3 2

(9a) (9b) (9c)

M. Rhein, T. Kalbe / Computer Aided Geometric Design 26 (2009) 825–841

829

Fig. 2. Left: coefficients on the outer layer of the half of an unfolded truncated octahedron T . Right: inner layer of T . ‘Red’ squares denote the determining set M on T . The other coefficients involved in the smoothness conditions (4)–(9) are shown as ‘blue’ dots, ‘green’ diamonds and ‘magenta’ stars. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Note that (4) describes the smoothness conditions across the triangular faces between neighboring truncated octahedra and (5)–(9) characterizes the smoothness conditions across triangular faces inside a single truncated octahedron. Since the smoothness conditions across the faces F 0 and F 3 degenerate to univariate conditions, we can describe these conditions simultaneously for the cases F 0S and F 0H , and F 3S and F 3H , respectively. The following lemma gives a determining set for the space S21 (). Here, T B denotes tetrahedra in  with the vertex v 3 lying on the boundary of Ω . T

T B Lemma 1. The set M := {ξ0011 , ξ0002 : T , T B ∈ } is a determining set of the space S21 ().

Proof. We have to show that for any spline s ∈ S21 (), with bξ = 0, ξ ∈ M, it follows that s ≡ 0. Let T ∈ , v 3 ∈ T and D 1 ( v 3 ) be the ball with radius 1 around v 3 . The spline s ∈ S20 () is C 1 -differentiable at v 3 iff all control points (ξ, bξ ), ξ ∈ D 1 ( v 3 ), lie in the same hyperplane in R4 . Since M ∩ D 1 ( v 3 ) contains exactly four points which do not lie on a plane in R3 , the four control points (ξ, bξ ), ξ ∈ M ∩ D 1 ( v 3 ) uniquely determine a hyperplane in R4 and it follows that bξ = 0, for all ξ ∈ D 1 ( v 3 ). This also follows by combining the smoothness conditions (4c), (6c), (7c), (8c) and (9c), and some elementary computations. Therefore bξ = 0, for all ξ = ξiT, j ,k, ∈ D ,   1. Now the smoothness conditions (5a)–(5c) imply that bξ = 0, for all ξ = ξiT, j ,k, ∈ D , k  1. Furthermore, the conditions (6a)–(6b) and (7a)–(7b) then imply bξ = 0, for all ξ = ξiT, j ,k, ∈ D , j  1. Finally, conditions (8a) and (9a) imply bξ = 0, for all ξ = ξiT, j ,k, ∈ D , (i , j , k, ) = (2, 0, 0, 0). Therefore it follows that bξ = 0, for all ξ ∈ D and s ≡ 0. 2 Fig. 2 shows the order in which the coefficients in the proof of Lemma 1 are determined. Due to symmetry, only a half unfolded truncated octahedron is shown. On the left side are the points on the boundary of the truncated octahedron and on the right side the points at distance one to the center vertex v 0 . We assume, that the truncated octahedron does not have any boundary vertices v 3 ∈ T B . The points of the determining set M are shown as ‘red’ squares ( ), while the points associated with the coefficients directly determined through the hyperplane-conditions and the coefficients bξ , ξ ∈ M, are shown as ‘blue’ dots ( ). Finally, ‘green’ diamonds ( ) and ‘magenta’ stars ( ) show the points which associated coefficients are determined through the univariate conditions (5a)–(5c) and bivariate conditions (6a)–(6b) and (7a)–(7b), respectively. 4. Quasi-interpolation scheme In this section we describe the construction of a family of quasi-interpolation operators on truncated octahedral partitions. We give a detailed description of the quasi-interpolation scheme and prove our main result, as well as some useful properties of certain operators. For n ∈ N, n  3, let

  V := v i jk = (ih, jh, kh): i , j , k = −1, . . . , (n + 1)

830

M. Rhein, T. Kalbe / Computer Aided Geometric Design 26 (2009) 825–841

Fig. 3. A truncated octahedron embedded in the cubic grid. The data points are located on the center of each hexagonal face.

be the cubic grid of (n + 3)3 points with grid size h ∈ R and 3 the truncated octahedral partition of truncated octahedra with radius h centered at the vertices of the BCC-grid



   := v i jk + h , h , h : i , j , k = 1, . . . , (n − 2), i , j , k all odd or i , j , k all even . V 2 2 2

We assume that the values of a function f ∈ C (Ω ∗ ), Ω ∗ = [−h, (n + 1)h]3 ⊆ R3 are known at the grid points of V . In each truncated octahedron T ∈ 3 lie exactly 8 points z1T , . . . , z8T of the set V located on the hexagonal faces of T (see Fig. 3),

which we denote by VT := { z1T , . . . , z8T }. Although, V contains data points that do not lie in 3, we assume in the following that these data points are also associated with truncated octahedra. Let  be the tetrahedral partition of 3 described in the previous section and Ω ⊂ Ω ∗ the covered region in R3 . In the following we describe our method to determine an approximating quadratic spline s f on  by setting all B-coefficients bξ , ξ ∈ D , to appropriate averages of the data. Thereby, for the calculation of each coefficient bξ we use only local data portions in the following sense. If the domain point ξ lies in the interior of a truncated octahedron A ∈ 3 then the associated coefficient bξ is deter), z A ∈ VA : mined by weighted averages of the 8 data values A i := f ( zA i i bξ =



w iA A i ,

w iA ∈ R.

i =1,...,8

Here, as well as in the following, the data points zT , i = 1, . . . , 8, of a truncated octahedron T are ordered by increasing i

(Euclidean) distances to the point ξ . This means z1T = minz∈VT ξ − z 2 and z8T = maxz∈VT ξ − z 2 , respectively. In the case of equal distance the order of the points can be chosen arbitrarily. Figs. 4 and 5 show a typical order of the data points for a selection of domain points on the outer layer of an unfolded truncated octahedron T . If the domain point ξ lies in the interior of a face of a truncated octahedron then the associated coefficient bξ is determined by weighted averages of the data values regarding the two truncated octahedra sharing this face. Let A, B ∈ 3 be two truncated octahedra sharing a face, then bξ , ξ ∈ int(A ∩ B ), is determined by bξ =



i =1,...,8

w iA A i +



w iB B i ,

w iA , w iB ∈ R,

i =1,...,8

where B i := f ( zB ), z B ∈ VB . Here, A denotes the truncated octahedron with ξ − z2A 2  ξ − z2B 2 . If A and B share a i i square face we have equal distances and the notation can be chosen arbitrarily. In the case that A and B share a hexagonal face, they have the common data point z1A = z1B and therefore the data value A 1 = B 1 is weighted by w 1A + w 1B .

M. Rhein, T. Kalbe / Computer Aided Geometric Design 26 (2009) 825–841

831

Fig. 4. The typical order of the data points z1T , . . . , z8T for the coefficients bξ , ξ = ξ0200 and ξ = ξ0110 (‘blue’ squares) on the outer layer of an unfolded truncated octahedron T . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

If the domain point ξ lies in the interior of an edge of a truncated octahedron then the associated coefficient bξ is determined by weighted averages of the data values regarding the three truncated octahedra sharing this edge. Let A, B , C ∈ 3 be three truncated octahedra sharing an edge, then bξ , ξ ∈ int(A ∩ B ∩ C ), is determined by bξ =



w iA A i +

i =1,...,8



w iB B i +

i =1,...,8



w iC C i ,

i =1,...,8

where w iA , w iB , w iC ∈ R, and C i := f ( zC ), zCi ∈ VC . Here, A denotes again the truncated octahedron with ξ − z2A 2 = i minT ∈3 ξ − z2T 2 . Since A and B , as well as A and C , share a hexagonal face, they have the common data points z1A = z1B and z2A = z1C . The data values A 1 = B 1 and A 2 = C 1 are then weighted by w 1A + w 1B and w 2A + w 1C , respectively. Finally, if the domain point ξ lies on a vertex of a truncated octahedron then the associated coefficient bξ is determined by weighted averages of the data values regarding the four truncated octahedra sharing this vertex. Let A, B , C , D ∈ 3 be the four truncated octahedra sharing the vertex v, then bξ , ξ = v, is determined by bξ =



i =1,...,8

w iA A i +



i =1,...,8

w iB B i +



i =1,...,8

w iC C i +



w iD D i ,

i =1,...,8

where w iA , w iB , w iC , w iD ∈ R, and D i := f ( zD ), z D ∈ VD . Due to symmetry the notation can be chosen arbitrarily here. Since i i exactly four hexagonal faces share this vertex, we also have common data points in this case. Without loss of generality we say z1A = z1B , z2A = z1C , z1D = z2B and z2D = z2C . In the following we describe a family of quasi-interpolation operators Q k : C (Ω ∗ ) → S 20 (), k  1 by giving the calculation rules for the B-coefficients of s f := Q k ( f ) in its piecewise representation (1) explicitly. Here, we assume that the data values come from a continuous function f ∈ C (Ω ∗ ). In some cases we have to distinguish which type of tetrahedra the coefficients are associated with. We denote by bξS B-coefficients associated with domain-points ξ ∈ T , where T is in T S and by bξH all other B-coefficients. Since often the same weights appear we use the compact notation A B 1,2 := A 1 + B 1 + A 2 + B 2 (analogous for A 1,2 , A BC D 1,2 , etc.). The determination of the B-coefficients is as follows: b0011 :=

1 −6   2 (k + 3) A 1,2 + (k + 1) A 3,4 + (k − 1) A 5,6 + (k − 3) A 7,8 2 k

 + 3 (k + 3) BC 1,2 + (k + 1) BC 3,4 + (k − 1) BC 5,6 + (k − 3) BC 7,8 ,    1 b0020 := 2−6 2 (k + 3) A 1,2 + k A 3,4,5,6 + (k − 3) A 7,8 + 3 (k + 3) BC 1 + (k + 2) BC 2,3 + (k + 1) BC 4 k + (k − 1) BC 5 + (k − 2) BC 6,7 + (k − 3) BC 8 ,  1 b0002 := 2−5 (k + 3) A BC D 1,2 + (k + 1) A BC D 3,4 + (k − 1) A BC D 5,6 + (k − 3) A BC D 7,8 ,

k 1 −4  := 2 (k + 3) A B 1,2 + (k + 1) A B 3,4 + (k − 1) A B 5,6 + (k − 3) A B 7,8 , k  1 S b0110 := 2−4 (k + 3) A B 1 + (k + 2) A B 2,3 + (k + 1) A B 4 + (k − 1) A B 5 + (k − 2) A B 6,7 + (k − 3) A B 8 , k  1 H b0110 := 2−4 (k + 3) A 1,2 + k A 3,4,5,6 + (k − 3) A 7,8 k S,H b0101

+ (k + 3) B 1 + (k + 2) B 2,3 + (k + 1) B 4 + (k − 1) B 5 + (k − 2) B 6,7 + (k − 3) B 8 ,

832

M. Rhein, T. Kalbe / Computer Aided Geometric Design 26 (2009) 825–841

Fig. 5. The typical order of the data points z1T , . . . , z8T for the remaining coefficients bξ (‘blue’ squares) on the outer layer of an unfolded truncated octahedron T . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

1 −4  (k + 2) A B 1,2,3,4 + (k − 2) A B 5,6,7,8 , 2 k  1 H b0200 := 2−4 (k + 3) A B 1 + (k + 1) A B 2,3,4 + (k − 1) A B 5,6,7 + (k − 3) A B 8 , k 1 −3  b1001 := 2 (k + 3) A 1,2 + (k + 1) A 3,4 + (k − 1) A 5,6 + (k − 3) A 7,8 , k 1 −3  S b1010 := 2 (k + 3) A 1 + (k + 2) A 2,3 + (k + 1) A 4 + (k − 1) A 5 + (k − 2) A 6,7 + (k − 3) A 8 , k  1 H b1010 := 2−3 (k + 3) A 1,2 + k A 3,4,5,6 + (k − 3) A 7,8 , k  1 S b1100 := 2−3 (k + 2) A 1,2,3,4 + (k − 2) A 5,6,7,8 , k  1 H b1100 := 2−3 (k + 3) A 1 + (k + 1) A 2,3,4 + (k − 1) A 5,6,7 + (k − 3) A 8 , k S b0200 :=

b2000 := 2−3 A 1,2,3,4,5,6,7,8 .

(10)

We note that all coefficients bξ , ξ ∈ D , are uniquely determined and therefore the spline s f is in will show that the weights are chosen carefully, so that s f is indeed globally C 1 on Ω .

S20 ().

The next result

Theorem 2. The quasi-interpolation operator Q k : C (Ω ∗ ) → S20 (), k  1 is a linear mapping into the space S21 (). Proof. We have to show that all conditions (4)–(9) are simultaneously satisfied. Since the complete proof is large and involves only elementary computations, we show the procedure exemplary on two typical conditions. At first we consider condition (4a) of univariate type in the case that F 0S lies on a common square face of two adjacent truncated octahedra: S S S b1100 + b˜ 1100 = 2b0200 .

Now we apply (10) for the involved coefficients and obtain on the left-hand side of the equation

1 −3  A 1,2,3,4 + (k − 2) A A 5, 6, 7, 8 . (k + 2) A 2 k  = B this is exactly two times the coefficient b S Since A 0200 and the condition is satisfied. Next, we consider condition (7a) T ∈ T H , inside a of bivariate type in the case that F 2H is the common triangular face of two adjacent tetrahedra T ∈ T H and  truncated octahedron: S H b1010 + b˜ 1010 =

1 2

H b1100 +

3 2

b1001 .

We note that one of the coefficients on the left-hand side is of type bξS and the other one of type bξH . Without loss of

S H and b˜ 1010 = b˜ 1010 . Now we apply (10) for the involved coefficients on the left-hand side of generality we say b1010 = b1010  = A we obtain the equation and considering that A

1 −3  2 (2k + 6) A 1 + (2k + 5) A 2 + (2k + 2) A 3 + (2k + 1) A 4 + (2k − 1) A 5 + (2k − 2) A 6 + (2k − 5) A 7 + (2k − 6) A 8 . k

M. Rhein, T. Kalbe / Computer Aided Geometric Design 26 (2009) 825–841

833

If we apply (10) to the coefficients on the right-hand side of the equation, we have to pay attention that only b1001 has the S H and b˜ 1010 . Therefore, for now we denote the data values regarding same ordering of the data points as the coefficients b1010

H ˆ 1 , . . . , Aˆ 8 and obtain b1100 by A

1 −4  (k + 3) Aˆ 1 + (k + 1) Aˆ 2,3,4 + (k − 1) Aˆ 5,6,7 + (k − 3) Aˆ 8 2 k

 + 3 (k + 3) A 1,2 + (k + 1) A 3,4 + (k − 1) A 5,6 + (k − 3) A 7,8 .

Considering the order of the data points we achieve the following correlations between the data values:

ˆ 1 = A1, A

ˆ 2 = A2, A

ˆ 3 = A3, A

ˆ 4 = A5, A

ˆ 5 = A4, A

ˆ 6 = A6, A

ˆ 7 = A7, A

ˆ 8 = A8. A

Applying this to the above term, we obtain 1 −4  (4k + 12) A 1 + (4k + 10) A 2 + (4k + 4) A 3 + (4k + 2) A 4 2 k

+ (4k − 2) A 5 + (4k − 4) A 6 + (4k − 10) A 7 + (4k − 12) A 8



and condition (7a) is satisfied. All other conditions (4)–(9) can be shown similarly. Only the order of the data points has to be taken into account for each coefficient individually. 2 We conclude this section with some important properties of the quasi-interpolation operators. Corollary 3. The quasi-interpolation operator Q k : C (Ω ∗ ) → S21 () is positive for k  3. Proof. For k  3 only positive weights of the data values appear in the calculation of the B-coefficients in (10). Therefore, assuming that the data values come from a positive function f ∈ C (Ω ∗ ), the quasi-interpolating spline s f := Q k ( f ) is also positive. This follows directly from the piecewise polynomial representation (1) and the non-negativity of the Bernstein polynomials. Finally, it is easy to see that f ≡ 0 implies s f ≡ 0. 2 The next lemma shows that the quasi-interpolation operators Q k , k  1, reproduce certain polynomials, which is essential for their approximation properties (see Section 5). In particular, the operator Q 2 reproduces quadratic polynomials up to a constant. Lemma 4. The following statements hold. (i) Q k (1) ≡ 1, ∀k  1; (ii) Q 2 ( p ) ≡ p , ∀ p ∈ P1 ; (iii) Q 2 ( p ) ≡ p , p ∈ {xy , xz, yz}; 2

(iv) Q 2 ( p ) ≡ p + h4 , p ∈ {x2 , y 2 , z2 }; (v) if p ∈ P2 , then there exists a constant C independent of h, such that Q 2 ( p ) ≡ p + Ch2 . Proof. (i) If all data values in the formulas (10) are equal to 1 it is easy to see that also all coefficients bξ are equal to 1. Therefore, it follows directly from (2), that Q k (1) ≡ 1, for all k  1. For (ii), we have to show that the quasi-interpolation operator Q 2 additionally to (i) reproduces the monomials {x, y , z}. Since Lemma 1 gives a determining set M for the space S21 () it is sufficient to prove the reproduction for the coefficients bξ , ξ ∈ M, only. Moreover, with (i) and the linearity of Q 2 , the reproduction is independent from the choice of the origin. Therefore we consider the tetrahedron T = [ v 0 , v 1 , v 2 , v 3 ] ∈ , where v 0 = (0, 0, 0), v 1 = (h, 0, 0), v 2 = (h, h4 , h4 ) and v 3 =

(h, h2 , 0). Here, T lies in a truncated octahedron with radius h centered at the origin. The coefficient bξ , ξ = ξ0002 , is only in M if ξ lies on the boundary of  and therefore we only need to consider the coefficient bξ ∈ M, where ξ = ξ0011 . Let p 1 (x, y , z) := x, p 2 (x, y , z) := y and p 3 (x, y , z) := z. Then these functions are univariate polynomials restricted to the edge i] [ v 2 , v 3 ] and the coefficients b[0011 , i = 1, 2, 3, in their representation (1) regarding T can be computed by [i ]

b0011 = 2p i (ξ0011 ) −

1 2



pi (v 2) + pi (v 3) ,

i = 1, 2, 3. [i ]

We obtain the following values for the coefficients b0011 , i = 1, 2, 3: [1]

b0011 = h,

[2]

b0011 =

3 8

h

[3]

and b0011 =

1 8

h.

834

M. Rhein, T. Kalbe / Computer Aided Geometric Design 26 (2009) 825–841

Now, we apply the order of increasing distances to the data points regarding the three truncated octahedra A, B and C with the common edge [ v 2 , v 3 ] and obtain



z1A

=

h h h



, h h 3h , z5A = , , 2 2 2  h h h z1B = , , , 2 2 2  h h h , z5B = − , , 2 2 2  3h h h z1C = , , , 2 2 2  5h h h , z5C = , , , ,



2 2 2

2

2 2



z2A

=

3h h h





z3A

, 3h h 3h , z6A = , , 2 2 2  h h h z2B = , ,− , 2 2 2  h h h , z6B = − , , − 2 2 2  3h h h z2C = , ,− , 2 2 2  5h h h C z6 = , ,− , 

2

2



, h 3h 3h , z7A = , , 2 2 2  h h h z3B = ,− , , 2 2 2  h h h , z7B = − , − , 2 2 2  3h h h z3C = ,− , , 2 2 2  5h h h , z7C = ,− ,

, ,

2 2

2

=

h 3h h



2

2

,

2

,

2

2

2 2



z4A

= 

z8A =

 z4B =



3h 3h h 2

,

2

,

2

3h 3h 3h 2 h

,

2

,

2

,

h

h

2

2

,− ,−

,

, h h h , z8B = − , − , − 2 2 2  3h h h z4C = ,− ,− , 2 2 2  5h h h C z8 = ,− ,− . 

2

2

2

2

Here, A is the truncated octahedron with radius h centered at (h, h, h), while B and C , are the truncated octahedra with radius h centered at (0, 0, 0) and (2h, 0, 0), respectively. Some elementary computations show that the calculation rule in [i ] (10) for the coefficient b0011 indeed reproduces the coefficient b0011 of p i , for i = 1, 2, 3. To prove (iii) and (iv) the procedure is similar. With (i), (ii) and the linearity of Q 2 , the reproduction is again independent from the choice of the origin. Now, a comparison of the calculated coefficient b0011 in (10) for Q 2 ( p ), [ p] p ∈ {xy , xz, yz, x2 , y 2 , z2 } and the coefficient b0011 of the polynomial p in its representation (1) regarding T shows that [ p]

[ p]

2

b0011 = b0011 for p ∈ {xy , xz, yz} and b0011 = b0011 + h4 for p ∈ {x2 , y 2 , z2 }. This proves (iii) and (iv), since Q 2 is linear and reproduces constants. Finally (v) follows directly by (i)–(iv) and the proof is complete. 2 The next corollary shows that the derivatives of certain polynomials are also reproduced by the quasi-interpolation operator Q 2 . In particular, the exact reproduction of the derivatives of quadratic polynomials is the reason that this operator β γ possesses the optimal approximation order for the derivatives of smooth functions. Here, we denote by D α x D y D z ( f ), where α , β, γ  0, the higher order partial derivatives of a sufficiently smooth function f . β

γ

γ

β

α Corollary 5. If p ∈ P2 , then D α x D y D z ( Q 2 ( p )) ≡ D x D y D z ( p ), for α + β + γ ∈ {1, 2}.

Proof. Let p ∈ P2 , then it follows from Lemma 4(v), that β

γ



β

γ



β

γ

α 2 Dα ≡ D αx D y D z ( p ), x D y D z Q 2 ( p ) ≡ D x D y D z p + Ch

where

α + β + γ ∈ {1, 2} and C is a constant independent of h. 2

Although, the quasi-interpolation operators Q k do not reproduce linear polynomials for k = 2, the next lemma shows that all operators Q k , k  1, reproduce the values of linear polynomials at the data points. Lemma 6. Let p ∈ P1 and k  1, then Q k ( p )( v ) = p ( v ), ∀ v ∈ V ∩ Ω . Proof. Let p ∈ P1 , v ∈ V ∩ Ω and T = [ v 0 , v 1 , v 2 , v 3 ] ∈  a tetrahedron with the vertex v 1 = v, then it follows that H H H Q k ( p )( v ) = b0200 , where b0200 is the B-coefficient of Q k ( p )| T in the representation (1). The coefficient b0200 in (10) is calculated by weighted averages of the data values A i and B i , i = 1, . . . , 8, at the data points ziA and ziB , i = 1, . . . , 8. We note that the order of the data points (see Section 3) can be chosen such that the data points ziA and ziB are point symmetric to v, i = 1, . . . , 8. This means ziA + ziB = 2v, i = 1, . . . , 8, and since the data values come from a linear polynomial it follows H that A i + B i = 2p ( v ), i = 1, . . . , 8. Applying this to the formula for the coefficient b0200 in (10) we obtain H = b0200

1 −3  2 (k + 3) + 3(k + 1) + 3(k − 1) + (k − 3) p ( v ) = p ( v ). k

2

The next result shows that the operators Q k , k  1, are uniformly bounded. This is an important result, since it guarantees that the computation of the splines is a stable process. Here as well as in the following sections we use the ∞-norm defined by

f Γ := max f (u ) , u ∈Γ

where Γ is a domain in R3 and f is a continuous function on Γ .

M. Rhein, T. Kalbe / Computer Aided Geometric Design 26 (2009) 825–841

835

Lemma 7. Let f ∈ C (Ω ∗ ) and T = [ v 0 , v 1 , v 2 , v 3 ] ∈ . Then we have

   Q k ( f )  k + 1 f Ω , T T k

k  1,

where Ω T is the cube with side length 32 h centered in v 3 . Proof. Since Ω T contains all data points needed for the calculation of the B-coefficients bξ , ξ ∈ D T , in (10) and k  1, we have



2k + 1 k + 1 f Ω T . |bξ |  max 1, , 2k

k

Then it follows from the non-negativity of the Bernstein polynomials and (2) that

     Q k ( f )  max |bξ |: ξ ∈ D T , T

and the proof is complete.

k1

2

5. Approximation properties In this section, we analyze the approximation properties of the quasi-interpolating splines Q k ( f ), k  1. We use the same notations as in the previous section. Theorem 8 gives general error bounds for the approximation of the values of smooth functions for k  1, and Theorem 9 gives special error bounds for the approximation of the values and derivatives of smooth functions for the case k = 2. In particular, we show that Q 2 ( f ) approximates the values of a sufficiently smooth function with order two, while simultaneously the first and second derivatives are optimally approximated with order two and one, respectively. In the following we let

  r     D ( f ) := max  D α D βy D γz ( f ) : α + β + γ = r , x Γ Γ

where Γ is a domain in R3 and f ∈ C r (Γ ). Theorem 8. Let f ∈ C (Ω ∗ ) and Q k ( f ), k  1, the quasi-interpolating spline defined in the previous section. Then the following statements hold. (i) If f ∈ C 1 (Ω ∗ ), then

     f − Q k ( f )  C k  D 1 ( f ) ∗ h, Ω Ω

where C k > 0 is an absolute constant independent of f and h. (ii) If f ∈ C 2 (Ω ∗ ), then



 

f − Q k ( f ) ( v )   C k  D 2 ( f )Ω ∗ h2 ,

∀ v ∈ V ∩ Ω,

where  C k > 0 is an absolute constant independent of f and h. Proof. Let f ∈ C 1 (Ω ∗ ), T = [ v 0 , v 1 , v 2 , v 3 ] ∈  an arbitrary tetrahedron with vertices as in Section 3, Ω T as in Lemma 7 and p f ∈ P0 the constant Taylor polynomial of f at v 3 . Then there exists an absolute constant C T > 0 independent of f and h such that

  f − p f ΩT  C T  D 1 ( f )Ω h. T

Furthermore, from the triangle inequality follows

     f − Q k ( f )  f − p f Ω +  Q k ( f ) − p f  . T T T

With Lemma 4(i), the linearity of Q k and Corollary 7 we obtain

     Q k ( f ) − p f  =  Q k ( f − p f )  k + 1 f − p f Ω . T T T k

Combining these inequalities leads to

     f − Q k ( f )  C k  D 1 ( f ) h, T Ω T

where C k = (1 +

k+1 )C T k 2 ∗

is an absolute constant independent of f and h. This proves statement (i).

Now, let f ∈ C (Ω ), v ∈ V ∩ Ω , T = [ v 0 , v 1 , v 2 , v 3 ] ∈  a tetrahedron with v 1 = v, Ω T as in Lemma 7 and p f ∈ P1 the linear Taylor polynomial of f at v 3 . Then there exists an absolute constant  C T > 0 independent of f and h such that





( f − p f )( v )   C T  D 2 ( f )Ω h2 . T

836

M. Rhein, T. Kalbe / Computer Aided Geometric Design 26 (2009) 825–841

Using the triangle inequality, we have





f − Q k ( f ) ( v )  ( f − p f )( v ) + Q k ( f ) − p f ( v )

and with Lemma 6, the linearity of Q k and Corollary 7 we obtain for the last term





Q k ( f ) − p f ( v ) = Q k ( f − p f )( v )  k + 1 f − p f Ω . T k

A combination of these inequalities leads to



 

f − Q k ( f ) ( v )   C k  D 2 ( f )Ω h2 , T

where  C k = (1 +

k+1  )C T k

is an absolute constant independent of f and h.

2

Theorem 9. Let f ∈ C (Ω ∗ ) and Q 2 ( f ) the quasi-interpolating spline for k = 2 defined in the previous section. Then the following statements hold. (i) If f ∈ C 2 (Ω ∗ ), then

     f − Q 2 ( f )  C  D 2 ( f ) ∗ h2 , Ω Ω

where C > 0 is an absolute constant independent of f and h. (ii) If f ∈ C 3 (Ω ∗ ), then

  r   D f − Q 2( f )    C  D 3 ( f )Ω ∗ h3−r , Ω

r = 1, 2,

where  C > 0 is an absolute constant independent of f and h. Proof. Let f ∈ C 2 (Ω ∗ ), T = [ v 0 , v 1 , v 2 , v 3 ] ∈  an arbitrary tetrahedron with vertices as in Section 3, Ω T as in Lemma 7 and p f ∈ P1 the linear Taylor polynomial of f at v 3 . Then there exist an absolute constant C T > 0 independent of f and h, such that

  f − p f ΩT  C T  D 2 ( f )Ω h2 . T

Furthermore, from the triangle inequality follows

     f − Q 2 ( f )  f − p f Ω +  Q 2 ( f ) − p f  . T T T

With Lemma 4(ii), the linearity of Q 2 and Corollary 7 for k = 2 we obtain

     Q 2 ( f ) − p f  =  Q 2 ( f − p f )  3 f − p f Ω . T T T 2

Combining these inequalities leads to

     f − Q 2 ( f )  C  D 2 ( f ) h2 , T Ω T

where C = (1 + is an absolute constant independent of f and h. This proves statement (i). Now let f ∈ C 3 (Ω ∗ ) and p f ∈ P2 the quadratic Taylor polynomial of f at v 3 . Then there exists an absolute constant  C T > 0 independent of f and h such that 3 )C T 2

 r     D ( f − p f )   C T  D 3 ( f )Ω h3−r , Ω T

T

r = 0, 1, 2.

Using the triangle inequality, we have

     r    D f − Q 2 ( f )    D r ( f − p f ) +  D r Q 2 ( f ) − p f  , T Ω T T

r = 1, 2,

and with Corollary 5 and the linearity of Q 2 we obtain for r = 1, 2,

 r     D Q 2( f ) − p f  = Dr Q 2( f − p f )  . T T

Now we apply the Markov-inequality for trivariate polynomials as well as Corollary 7 for k = 2 and obtain

    r    D Q 2 ( f − p f )   C M  Q 2 ( f − p f )  C M 3 f − p f Ω , T T T r r h

h 2

r = 1, 2,

where  C M > 0 is an absolute constant independent of f and h. A combination of these inequalities leads to

 r    D f − Q 2( f )    C  D 3 ( f )Ω h3−r , T T

r = 1, 2,

where  C = (1 + 32  C M ) C T is an absolute constant independent of f and h.

2

M. Rhein, T. Kalbe / Computer Aided Geometric Design 26 (2009) 825–841

837

Table 1 Approximation of blob by sblob for k = 3. 1/h

errblob data

errblob max

errblob mean

errblob rms

8 16 32 64 128

0.07232590 0.01885980 0.00476562 0.00119498 0.00029940

0.0809447 0.0237130 0.0083607 0.0035041 0.0016146

0.01172100 0.00398319 0.00150337 0.00065814 0.00031451

0.01859390 0.00583713 0.00207445 0.00089226 0.00042672

Table 2 Approximation of fr by sfr for k = 3. fr

fr

fr

fr

1/h

errdata

errmax

errmean

errrms

8 16 32 64 128

0.2004200 0.0565250 0.0145831 0.0036855 0.0009267

0.2005360 0.0658433 0.0196277 0.0073331 0.0031859

0.01000050 0.00339222 0.00118928 0.00048196 0.00022156

0.02202190 0.00722884 0.00238082 0.00091769 0.00041459

Table 3 Approximation of ml by sml for k = 3. 1/h

errml data

errml max

errml mean

errml rms

8 16 32 64 128

0.0935998 0.0913610 0.0453692 0.0137861 0.0033092

0.1652220 0.1020680 0.0483350 0.0151583 0.0048792

0.0607474 0.0377678 0.0126823 0.0034744 0.0010035

0.0712580 0.0450034 0.0158189 0.0043887 0.0012671

Table 4 Approximation of blob by sblob for k = 2. 1/h

errblob data

errblob max

errblob mean

errblob rms

8 16 32 64 128

0.05344710 0.01372490 0.00345431 0.00086509 0.00021647

0.05371940 0.01378350 0.00346625 0.00086701 0.00021684

0.00787331 0.00222180 0.00058516 0.00014987 0.00003791

0.01309350 0.00357384 0.00092267 0.00023371 0.00005877

6. Numerical tests In this section, the approximation properties of our splines are demonstrated with numerical examples based on smooth test functions f . To do this, we use synthetic data sampled at the grid points v ∈ V (see Section 4) to calculate the splines s f = Q k ( f ) on the partition ♦ of Ω . In our tests, we consider three functions. The blobby function 2 2 2 2 2 2 blob( v ) = e −5((x−1/2) + y +z ) + e −5(x +( y −1/2) +z )

is a smooth blending of two metaballs (Wyvill et al., 1986). As a second representative for an exponential, we have the more complex function of Franke type (Franke, 1979), 2 2 2 2 2 fr( v ) = 1/2e −10((x−1/4) +( y −1/4) ) + 3/4e −16((x−1/4) +( y −1/4) +(z−1/4) )

+ 1/2e −10((x−3/4)

2

+( y −1/8)2 +(z−1/2)2 )

− 1/4e −20((x−3/4)

2

+( y −3/4)2 )

.

Finally, we provide the results for the well-known Marschner Lobb (Marschner and Lobb, 1994) function ml( v ) =

1 − sin(π z/2) + α (1 + cos(2π f M cos(π



x2 + y 2 /2)))

2(1 + α )

,

1

α= , 4

f M = 6,

which is highly oscillating and therefore a challenging test. For each function, we proceed with a series of tests on the cubic domain Ω ∗ = [−h, 1 + h]3 ⊆ R3 translated by (− 12 , − 12 , − 12 ) with decreasing h. The first column in Tables 1 to 12 shows the increasing values of 1/h, while the remaining columns show different errors for the function values (Tables 1 to 6) and the first and second derivatives in x-direction (Tables 7 to 12), respectively. The maximal error of the function values and the derivatives at the grid points is denoted as









errdata = max D rx f ( v ) − s f ( v ) : v ∈ V ,

r = 0, 1, 2,

838

M. Rhein, T. Kalbe / Computer Aided Geometric Design 26 (2009) 825–841

Table 5 Approximation of fr by sfr for k = 2. fr

fr

fr

fr

1/h

errdata

errmax

errmean

errrms

8 16 32 64 128

0.1573410 0.0429864 0.0109869 0.0027619 0.0006914

0.15738500 0.04300920 0.01099240 0.00276327 0.00069177

0.00704709 0.00213871 0.00059565 0.00015743 0.00004047

0.01612600 0.00481450 0.00131270 0.00034245 0.00008743

Table 6 Approximation of ml by sml for k = 2. 1/h

errml data

errml max

errml mean

errml rms

8 16 32 64 128

0.0883974 0.0811994 0.0312147 0.0085849 0.0020048

0.1655600 0.0992450 0.0341531 0.0089289 0.0021327

0.06017160 0.03329990 0.00969430 0.00241578 0.00059512

0.0718670 0.0403393 0.0120682 0.0030287 0.0007481

Table 7 Approximation of D x (blob) by D x (sblob ) for k = 2. D (blob)

D (blob)

D (blob)

D (blob)

1/h

x errdata

x errmax

x errmean

x errrms

8 16 32 64 128

0.1585890 0.0411008 0.0104755 0.0026472 0.0007590

0.2417030 0.0712189 0.0188899 0.0047828 0.0012712

0.02473250 0.00708406 0.00188885 0.00048762 0.00012448

0.0385371 0.0106428 0.0027800 0.0007095 0.0001801

Table 8 Approximation of D x (fr ) by D x (sfr ) for k = 2. D (fr )

D (fr )

D (fr )

D (fr )

1/h

x errdata

x errmax

x errmean

x errrms

8 16 32 64 128

0.7384930 0.2107540 0.0574643 0.0145658 0.0036746

0.9166850 0.3006540 0.0829841 0.0212066 0.0053748

0.04234930 0.01333240 0.00370368 0.00097108 0.00024830

0.09466990 0.02900140 0.00791414 0.00205638 0.00052339

for a function f and errmax is the (approximate) maximal error of s f in the uniform norm on Ω . The latter error is computed on a fine discretization X of the domain by choosing a fixed number of 220 evenly distributed sample points in each tetrahedron:









errmax = max D rx f ( v ) − s f ( v ) : v ∈ X ,

r = 0, 1, 2.

For the sake of completeness, we also provide the mean error errmean errmean =

1 

|X |



D r f ( v ) − s f ( v ) , x

r = 0, 1, 2,

v ∈X

as well as the root mean square error



errrms =

1 

|X |



D rx f ( v ) − s f ( v ) 2 ,

r = 0, 1, 2,

v ∈X

on the same fine discretization X . As can be seen from Tables 1, 2 and 3, the approximation errors errmax , errmean and errrms for s f = Q 3 ( f ) decrease by about a factor of two for each function as the grid size goes down from h to h/2, while the error at the data points decreases by a factor of four. This confirms the results from Theorem 8, Section 5. Tables 4 to 6 show the approximation errors for s f = Q 2 ( f ), while Tables 7 to 12 show the approximation errors of the first and second derivatives in x-direction D x ( f ) and D xx ( f ), respectively. Here, the results from Theorem 9 are confirmed, since the error of the function values and the first derivatives decreases by a factor of four as the grid size goes down from h to h/2 and the error of the second derivatives decreases by a factor of two. This confirms our result that the operator Q 2 approximates the derivatives of a sufficiently smooth functions with optimal approximation order.

M. Rhein, T. Kalbe / Computer Aided Geometric Design 26 (2009) 825–841

839

Table 9 Approximation of D x (ml) by D x (sml ) for k = 2. D (ml)

D (ml)

D (ml)

D (ml)

1/h

x errdata

x errmax

x errmean

x errrms

8 16 32 64 128

5.03099 4.21588 1.52074 0.42781 0.10828

5.33248 4.58655 2.78363 0.80967 0.21132

1.575060 0.887093 0.254649 0.063386 0.015481

2.042270 1.254340 0.382565 0.097711 0.024102

Table 10 Approximation of D xx (blob) by D xx (sblob ) for k = 2. D (blob)

D (blob)

D (blob)

D (blob)

1/h

xx errdata

xx errmax

xx errmean

xx errrms

8 16 32 64 128

2.781800 1.513700 0.783794 0.398972 0.225402

7.36277 4.48919 2.34272 1.18080 0.61108

0.647459 0.366841 0.192558 0.098316 0.049920

0.944675 0.525985 0.273797 0.139153 0.070501

Table 11 Approximation of D xx (fr) by D xx (sfr ) for k = 2. D (fr )

D (fr)

D (fr )

D (fr)

1/h

xx errdata

xx errmax

xx errmean

xx errrms

8 16 32 64 128

9.48307 6.56458 3.35031 1.63083 0.81025

34.2698 24.2788 13.2468 6.82766 3.45017

1.412130 0.876069 0.480231 0.249896 0.127268

2.726820 1.685340 0.916512 0.474181 0.240706

Table 12 Approximation of D xx (ml) by D xx (sml ) for k = 2. D (ml)

D (ml)

D (ml)

D (ml)

1/h

xx errdata

xx errmax

xx errmean

xx errrms

8 16 32 64 128

223.618 147.363 63.4782 29.2791 13.8998

258.992 314.767 306.433 209.945 118.557

55.6158 43.9654 25.2862 12.7598 6.2608

77.3912 66.5142 42.1393 22.2219 11.0415

In all cases, the computation of the splines s f , even on the finest partition of Ω involving more than two million data samples, were computed in less than one minute on a standard PC, depending on the complexity of f . In the numerical tests, several techniques for fast evaluation of s f can be employed, e.g. de Casteljau’s algorithm, which provides the functional value and the derivatives simultaneously. 7. Future work and remarks Recently, a scheme for high-quality visualizations of isosurfaces from cubic trivariate C 1 -splines on type-6 partitions has been introduced (Kalbe and Zeilfelder, 2008). In addition, the authors give a visual comparison of the smooth cubic splines and the quadratic supersplines of Rössl et al. (2004), defined on that same partition. Because of the missing C 1 continuity of quadratic supersplines, visual artifacts are more apparent compared to the cubic splines. Interactive frame-rates for reasonable sized data were achieved with both methods by exploiting the high parallelism of modern programmable graphic processors. We are currently working on a similar approach for the visualization of our quasi-interpolation scheme. Due to the high symmetry of the underlying truncated octahedral partition and the fact that our splines have a lower degree, we expect better performance than for the cubic splines in Kalbe and Zeilfelder (2008), while at the same time visual quality is retained. The BCC-Lattice has gained increased popularity in numerous applications, e.g. volume rendering (Csebfalvi, 2005), skeletonization (Brunner and Brunnett, 2008) and topology preserving reconstruction (Stelldinger and Strand, 2006). It is well known that the BCC-Lattice is a superior sampling grid compared to the cubic lattice and therefore equal approximation quality can be achieved with fewer data samples. The nature of C 1 -splines on truncated octahedral partitions motivates the direct usage of data on the BCC-Lattice. Currently, we are working on an adaptation of our approximation scheme to such data sets. In this context, we are confident to achieve even better approximation properties, in particular an optimal order of approximation.

840

M. Rhein, T. Kalbe / Computer Aided Geometric Design 26 (2009) 825–841

In a next step, we analyze the structure of C 2 -splines on truncated octahedral partitions. First investigations on this topic have shown that local quasi-interpolation operators with an expected degree of four could be constructed. To our knowledge, this would be the first C 2 -smooth spline approximation scheme with a lower degree than five. Remark 10. The approximation properties in Section 5 show that only the operator Q 2 has an approximation order of two for the values of a sufficiently smooth function, while all other operators Q k , k  1, have an approximation order of one. However, in our experience data sets produced by scanning devices, such as CT or MRI, often contain large noise. Therefore, the higher approximation order of the operator Q 2 will not automatically result in a better data representation or higher visual quality. Due to their construction all quasi-interpolation operators Q k , k  1, have a different denoising effect, which is specified by the value of k. Therefore, it depends mainly on the data used, which operator should be preferred. Of special interest may be the positive operator Q 3 , since several weights in the computation of the spline coefficients in (10) become zero. The reduced number of data points needed for the calculation of the coefficients results in a noticeable faster computation of the whole spline. Remark 11. The quasi-interpolation scheme for C 1 -splines on truncated octahedral partitions described in Section 4 can be easily extended to a cubic domain by the usage of half, quarter and eighth truncated octahedra. Thereby, no additional data points are needed. However, to guarantee the approximation order of the splines at the boundary, we have to use data points outside of the domain (see Section 4). But these additional data points needed can also be constructed by local extrapolation at the boundary of the given data, such that the approximating spline covers the whole initial data volume. Remark 12. It is possible to modify our construction of an approximating spline to an interpolating spline with a simple trick. We assume that the data lie on a cubic grid with grid-size h. Then we construct a truncated octahedral partition with size h, where half of the truncated octahedra are centered at the data points. In these truncated octahedra we apply the data value in the center to the eight data points on the hexagonal faces, needed for our quasi-interpolation scheme. Therefore, we have constructed all data points needed for the calculation of the spline on the whole partition and the resulting spline interpolates the initial data values. However, we lose the approximation properties by this approach. Remark 13. We have also applied our scheme for k = 2 to the simple test functions f 1 (x, y , z) = xy + xz + yz + x + y + z + 1

and

f 2 (x, y , z) = x2

in order to verify the results of Lemma 4 and Corollary 5. To be more precise, the errors of s f 1 = Q 2 ( f 1 ) are (numerical) zero for the functional values and the derivatives, confirming Lemma 4(i), (ii) and (iii). The error of s f 2 = Q 2 ( f 2 ) is constantly h2 /4, in accordance with Lemma 4(iv) (and (v)), while the errors for the corresponding derivatives are zero, which confirms Corollary 5. For the trilinear function f 3 (x, y , z) = axyz our tests show that the errors of s f 3 = Q 2 ( f 3 ) and D x (s f 3 ) decrease by a factor of 8 for h → h/2, meaning that in this case, we have the optimal approximation order three. References Bajaj, C., 1999. Data Visualization Techniques. John Wiley and Sons. Brunner, D., Brunnett, G., 2008. Fast force field approximation and its application to skeletonization of discrete 3D objects. Computer Graphics Forum 27 (2), 261–270. Csebfalvi, B., 2005. Prefiltered Gaussian reconstruction for high-quality rendering of volumetric data sampled on a body-centered cubic grid. Visualization, VIS 05, 311–318. de Boor, C., Fix, G.J., 1973. Spline approximants by quasiinterpolants. J. Approx. Theory 8, 19–45. de Boor, C., 1987. B-form basics. In: Farin, G. (Ed.), Geometric Modelling. SIAM, pp. 131–148. Farin, G., 1986. Triangular Bernstein–Bézier patches. Comput. Aided Geom. Design 3, 83–127. Franke, R., 1979. A critical comparison of some methods for interpolation of scattered data. Naval Postgraduate School Tech. Rep. NPS-53-79-003. Hangelbroek, T., Nürnberger, G., Rössl, C., Seidel, H.-P., Zeilfelder, F., 2004. Dimension of C 1 -splines on type-6 tetrahedral partitions. J. Approx. Theory 131, 157–184. Hecklin, G., Nürnberger, G., Zeilfelder, F., 2007. The structure of C 1 spline spaces on Freudenthal partitions. SIAM J. Math. Anal. 2, 347–367. Kalbe, T., Zeilfelder, F., 2008. Hardware-accelerated, high-quality rendering based on trivariate splines approximating volume data. Computer Graphics Forum 27 (2), 331–340. Lai, M.-J., Le Méhauté, A., 2004. A new kind of trivariate C 1 spline. Advances in Comp. Math. 21, 273–292 (special issue: multivariate splines). Lai, M.-J., Schumaker, L.L., 2007. Spline Functions on Triangulations. Cambridge University Press. Marschner, S., Lobb, R., 1994. An evaluation of reconstruction filters for volume rendering. In: Proc. of IEEE Vis., pp. 100–107. Meissner, M., Huang J., Bartz D., Mueller, K., Crawfis, R., 2000. A practical comparison of popular volume rendering algorithms. In: Symposium on Volume Visualization and Graphics 2000, pp. 81–90. Nürnberger, G., Rössl, C., Seidel, H.-P., Zeilfelder, F., 2005. Quasi-interpolation by quadratic piecewise polynomials in three variables. Comput. Aided Geom. Design 22, 221–249. Rössl, C., Zeilfelder, F., Nürnberger, G., Seidel, H.-P., 2004. Reconstruction of volume data with quadratic super splines. In: van Wijk, J.J., Moorhead, R.J., Turk, G. (Eds.), Trans. Visualization and Computer Graphics 10 (4), 397–409.

M. Rhein, T. Kalbe / Computer Aided Geometric Design 26 (2009) 825–841

841

Schumaker, L.L., Sorokina, T., Worsey, A.J., 2009. A C 1 -quadratic trivariate macro-element space defined over arbitrary tetrahedral partitions. Journal of Approximation Theory 158 (1), 126–142. Schumaker, L.L., Sorokina, T., 2004. C 1 quintic splines on type-4 tetrahedral partitions. Advances in Comp. Math. 21, 421–444 (special issue: multivariate splines). Sorokina, T., Worsey, A.J., 2008. A multivariate Powell–Sabin interpolant. Adv. Comput. Math. 29 (1), 71–89. Sorokina, T., Zeilfelder, F., 2007. Local quasi-interpolation by cubic C 1 splines on type-6 tetrahedral partitions. IMJ Numerical Analysis 27, 74–101. Stelldinger, P., Strand, R., 2006. Topology preserving digitization with FCC and BCC grids. In: Combinational Image Analysis. In: Lecture Notes in Computer Science, vol. 4040. Springer, Berlin/Heidelberg, pp. 226–240. Theußl, T., Möller, T., Hladuvka, J., Gröller, M., 2002. Reconstruction issues in volume visualization. In: Data Visualization: State of the Art, Proc. of IEEE Visualization 2002, pp. 1–8. Worsey, A.J., Farin, G., 1987. An n-dimensional Clough–Tocher interpolant. Constr. Approx. 3, 99–110. Worsey, A.J., Piper, B., 1988. A trivariate Powell–Sabin interpolant. Comput. Aided Geom. Design 5, 177–186. Wyvill, G., McPhetters, C., Wyvill, B., 1986. Data structures for soft objects. The Visual Computer 2, 227–234. Ženišek, A., 1973. Polynomial approximation on tetrahedrons in the finite element method. J. Approx. Theory 7, 334–351.