Multivariate spline approximation of the signed distance function

Journal of Computational and Applied Mathematics 265 (2014) 276–289

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Multivariate spline approximation of the signed distance function Ren-Hong Wang a , Qing-Jie Guo a,∗ , Chun-Gang Zhu a , Chun-Jing Li b a

School of Mathematical Sciences, Dalian University of Technology, Dalian 116024, China

b

Department of Mathematics, Tongji University, Shanghai 200092, China

article

info

Article history: Received 12 October 2012 Received in revised form 27 February 2013 Dedicated to Professor Benyu Guo on the occasion of his 70th Birthday

abstract The signed distance function can effectively support many geometry processing tasks such as decimates, smoothing and shape reconstruction since it provides efficient access to distance estimates. In this paper, we present an adaptive method to approximate the signed distance function of a smooth curve by using polynomial splines over type-2 triangulation. The trimmed offsets are also studied. © 2013 Elsevier B.V. All rights reserved.

Keywords: Multivariate spline Signed distance function Type-2 triangulation Trimmed offset Quasi-interpolation

1. Introduction The signed distance function (SDF) of a closed curve Γ in the plane assigns to any point x the shortest distance between x and Γ , with a positive sign if x is inside Γ and a negative one otherwise [1]. The SDF is widely used to deal with the problems in image processing, geometric computing, surface reconstruction, etc. [2,3]. Russo et al. [4] presented a method for the reconstruction of the SDF in the context of level set methods. It was a modification of the algorithm which made use of the PDE equation

∂φ = sgn(φ 0 )(1 − |∇φ|), ∂t φ(x, 0) = φ 0 (x) for the distance function introduced by Sussman et al. [5]. It was based mainly on the use of a truly upwind discretization near the interface. Belytschko et al. [6] developed an algorithm for smoothing the surfaces in finite element formulations of contact-impact. The smoothed signed distance functions were constructed by a moving least-squares approximation with a polynomial basis. There is a close relationship between SDF and the level set method. For example, the distance function preservation has several advantages in the geometric and numerical point of view [7–9]. Song et al. [1,10] presented a method to approximate the SDF of a smooth curve or surface and compute the medial axis by using polynomial splines over hierarchical T -meshes (PHT splines); the theory about PHT splines refers to Deng et al. [11].

∗

Corresponding author. E-mail addresses: [email protected] (R.-H. Wang), [email protected] (Q.-J. Guo), [email protected] (C.-G. Zhu), [email protected] (C.-J. Li). 0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.09.043

R.-H. Wang et al. / Journal of Computational and Applied Mathematics 265 (2014) 276–289

277

Offset curve is defined as the locus of the points which are at constant distance along the normal from the given curve. It is used as the CNC machine tool path in the CAD/CAM. There are many literatures dealing with various aspects of offset curves [12,13]. It is well known that the offsetting may produce curves with self-intersections, which have to be detected and eliminated in applications. After eliminating these unwanted segments of the offset curve, trimmed offset curve is obtained [1]. In this paper, we propose a method to compute the approximation of the SDF based on the quasi-interpolation operator (2) 2) ¯ (mn ). This method is simple since the derivatives are not required. We can on the bivariate spline space S21 (∆mn ) and S21 (∆ approximate the SDF using non-uniform spline space with the adaptive refinement. Moreover, this method can be applied to compute an implicit approximation of the trimmed offset curves of any given curve. The paper is organized as follows. In (2) Section 2, we introduce the definition and properties of SDF. We recall some facts about the bivariate spline space S21 (∆mn ) (2)

¯ mn ), and their quasi-interpolation operators in Section 3. In Section 4, we present a method to construct the and S21 (∆ approximation of the SDF by bivariate splines. Some examples are provided in Section 5. Finally, we conclude the paper in Section 6. 2. The signed distance function Consider a simple closed curve Γ in the Euclidean plane R2 . The curve Γ divides R2 into three parts: the interior Γ + , the exterior Γ − , and Γ . We define the distance function, and the signed distance function of Γ as follows:

Definition 1 ([14]). Let Γ be a simple closed curve in the Euclidean plane R2 . The distance function to Γ is the function x ∈ R2 → d(x, Γ ) defined by: d(x, Γ ) = min |x − y|.

(2.1)

y∈Γ

Definition 2 ([14]). Let Γ be a simple closed curve in the Euclidean plane R2 . The signed distance function to Γ is the function x ∈ R2 → u(x, Γ ) defined by:

 −d(x, Γ )

u(x, Γ ) = 0 d(x, Γ )

if x ∈ Γ + , if x ∈ Γ , if x ∈ Γ − .

(2.2)

Distance functions have a kink at the interface where d = 0 is a local minimum, causing problems in approximation derivatives on or near the interface. On the other hand, SDFs are monotonic across the interface and can be differentiated there with significantly higher confidence. SDFs share all the properties of implicit functions. In addition, there are a number of new properties that only SDFs possess. The main property of SDF is

|∇ u| = 1.

(2.3)

Eq. (2.3) is true only in a general sense. It is not true for points that are equidistant from at least two points on the interface. This set of points is often called a skeleton in the literature and is a zero measure set. Definition 3 ([15]). The set {x ∈ R2 | there exist at least two distinct points y and z in Γ such that ∥x − y∥ = ∥x − z ∥ = d(x, Γ )} is called the skeleton of Γ . Some interesting properties of SDF are listed below (cf. [14]): (1) The unit normal vector:

∇u ∇u − → = = ∇ u. n = |∇ u| 1 (2) Mean curvature:

κ=∇·



∇u |∇ u|



= ∇ · (∇ u) = 1u.

(3) Closest point on interface: given x ∈ R2 , then the closest point y on Γ is

− →

y = x − u(x) n . (4) Differentiability on Γ : u is differentiable on Γ almost everywhere. (5) Convex: if Γ + is a convex region, then u is a convex function.

278

R.-H. Wang et al. / Journal of Computational and Applied Mathematics 265 (2014) 276–289

Fig. 1. The uniform type-2 triangulation.

(2)

( 2)

¯ mn ) 3. An introduction to bivariate spline space S21 (∆mn ) and S21 (∆ (2)

In this section, we will introduce the bivariate spline space, the spline space over uniform type-2 triangulation S21 (∆mn ), and the spline space over non-uniform type-2 triangulation

S21

(2)

¯ mn ). (∆

Definition 4 ([16]). For a simply connected domain D ⊂ R2 , let ∆ be a partition of the domain D consisting of finite irreducible algebraic curves, and Di , i = 1, 2, . . . , T be all of the cells. For integers k and µ with k > µ ≥ 0, we define µ

Sk (∆) = {s(x, y) ∈ C µ (D)|s(x, y)|Di ∈ Pk , ∀Di , i = 1, 2, . . . , T }. It is called the bivariate spline space with degree k and smoothness µ over the domain D with respect to the partition ∆, where Pk denotes the collection of all bivariate polynomials of real coefficients with total degree k. Wang discovered the following fundamental theorem on bivariate splines (cf. [17]). µ

Theorem 1 (Cf. [17]). s(x, y) ∈ Sk (∆) if and only if the following conditions are satisfied: 1. for each interior edge of ∆, which is defined by Γi : li (x, y) = 0, there exists the so-called smoothing cofactor qi (x, y) such that µ+1

pi,1 − pi,2 = li

qi ,

where the polynomials pi,1 and pi,2 are determined by the restriction of s(x, y) on the two cells Di,1 and Di,2 with Γi as the common edge, and qi ∈ Pα−(µ+1)ni , α = max(degrees of pi,1 and pi,2 ), ni = degree of li ; 2. for any interior vertex vj of ∆, the following conformality conditions are satisfied

 (j) [li (x, y)]µ+1 qi(j) (x, y) ≡ 0, where the summation is taken over all the interior edges Γi are refixed in such a way that when a point crosses Γi

(j)

(j)

(j)

passing though vj , and the signs of the smoothing cofactors qi

from Di,2 and Di,1 , it goes around vj in a counter-clockwise manner.

In this paper, without loss of generality, we consider the rectangular region D = [0, 1] × [0, 1]. (2)

3.1. Bivariate spline space S21 (∆mn ) over the uniform type-2 triangulation (2)

The uniform type-2 triangulation ∆mn of D is yielded by the following partition lines (see Fig. 1): mx − i = 0,

ny − i = 0,

where i = · · · , −1, 0, 1, . . . .

ny − mx − i = 0,

ny − mx − i = 0,

R.-H. Wang et al. / Journal of Computational and Applied Mathematics 265 (2014) 276–289

(a) G: Support of B(x, y).

279

(b) B-spline base. (2)

Fig. 2. The support and B-spline base B(x, y) over ∆mn .

Using the dimension formulae on the cross-cut partition [16,18], we have 2) dim S21 (∆(mn ) = (m + 2)(n + 2) − 1.

(3.1) (2)

We introduce a locally supported spline in S21 (∆mn ) whose supported octagon G centers at (0, 0) as shown in Fig. 2(a). It is known that a bivariate polynomial of degree 2 on a triangle can be uniquely determined by the values of three vertices, and three midpoints of the edges. In Fig. 2(a), the values are denoted on some triangles, and other values can be obtained by symmetry of lines x = 0, y = 0, x + y = 0, x − y = 0. Let B(x, y) be a piecewise polynomial defined on G with the values shown in Fig. 2(a). Clearly, B(x, y) ∈ C 1 (R2 ), and B(x, y) > 0 inside of G (Fig. 2(b)). Hence, B(x, y) is a bivariate B-spline over (2) the triangulation ∆mn . Moreover, we can point out the support of B(x, y) is the smallest one [16,18]. Let



Bi,j (x, y) = B mx − i +

1 2

, ny − j +

1



2

,

(3.2)

then

A = {Bi,j (x, y) : i = 0, . . . , m + 1; j = 0, . . . , n + 1} (2)

(2)

is a subset of S21 (∆mn ). Note that each element of A is a nontrivial element of S21 (∆mn ), and the number of elements in A is (m + 2)(n + 2). By (3.1), A must be a linearly dependent set. There are the following results (cf. [16]). Theorem 2 ([16]). The bivariate B-spline functions in A defined by (3.2) satisfy m+1 n+1

 (−1)i+j Bi,j (x, y) = 0. i=0 j=0

For any i0 , j0 , 0 ≤ i0 ≤ m + 1, 0 ≤ j0 ≤ n + 1, the collection

Ai0 ,j0 = {Bi,j (x, y) ∈ A : (i, j) ̸= (i0 , j0 )} (2)

is a basis of S21 (∆mn ). Theorem 3 ([16]). Bivariate B-spline functions in A have the property m+1 n+1

 i=0 j=0

Bi,j (x, y) ≡ 1.

(3.3)

280

R.-H. Wang et al. / Journal of Computational and Applied Mathematics 265 (2014) 276–289

Fig. 3. The non-uniform type-2 triangulation.

(2)

By Theorem 2, for each bivariate spline s(x, y) ∈ S21 (∆mn ), there must exist ci,j ∈ R, 0 ≤ i ≤ m + 1, 0 ≤ j ≤ n + 1, such that m+1 n+1

s(x, y) =



ci,j Bi,j (x, y).

i=0 j=0

(2)

(2)

Considering two bivariate quasi-interpolation operators (cf. [16]) Vmn : C (D) → S21 (∆mn ), Wmn : C (D) → S21 (∆mn ) defined as follows: m+1 n+1

Vmn (f ) =



f (xi , yj )Bi,j (x, y),

(3.4)

i=0 j=0 m+1 n+1

Wmn (f ) =



λi,j (f )Bi,j (x, y),

(3.5)

i =0 j =0

 −1 , 2j2n is the center of the support of Bi,j (x, y), and            1 i j i−1 j i−1 j−1 i j−1 2i − 1 2j − 1 − f +f +f +f . λi,j (f ) = 2f , , , , ,

where (xi , yj ) =

 2i−1 2m

2m

2n

4

m n

m

n

m

n

m

n

Theorem 4 (Cf. [16]). Vmn (f ) ≡ f , Wmn (f ) ≡ f ,

f = 1, x, y, xy,

∀f ∈ P 2 . (2)

¯ mn ) over the non-uniform type-2 triangulation 3.2. Bivariate spline space S21 (∆ The non-uniform type-2 triangulation (see Fig. 3) of D is yielded over the original non-uniform rectangular partition 0 = x0 < x1 < · · · < xm = 1, 0 = y0 < y1 < · · · < yn = 1. Denote by hi = xi − xi−1 , kj = yj − yj−1 , where hi are unnecessary to be the same as each other, so do kj .

R.-H. Wang et al. / Journal of Computational and Applied Mathematics 265 (2014) 276–289

(a) The support of non-uniform Bi,j (x, y).

281

(b) Non-uniform B-spline Bi,j (x, y). (2)

¯ mn . Fig. 4. The support and B-spline base Bi,j (x, y) over ∆

Let Ai =

hi hi + hi+1

,

A′i =

hi +1 hi + hi+1

,

Bj =

kj kj + kj+1

,

B′j =

kj+1 kj + kj+1

.

As we know that a bivariate quadric polynomial defined on any given triangle can be determined by its values at three vertices, and midpoints of 3 edges of this triangle. So, for all −1 ≤ i ≤ m and −1 ≤ j ≤ n, every non-uniform B-spline 2) ¯ (mn ), (u, v) ∈ [xi−1 , xi+2 ] × [yj−1 , yj+2 ] can be represented as follows (see Fig. 4). Bi,j (x, y) ∈ S21 (∆ 1. 7.

1 16 1

B′j+1 ;

2.

1 16

3.

B′j+1 ;

(5 + 4A′i+1 )B′j+1 ;

1

1 4

Ai B′j+1 ;

4.

1 ′ Bj+1 ; 4

5. 1

1 ′ Ai+1 B′j+1 ;

4

6.

1 16

(5 + 4Ai )B′j+1 ;

(2 + Ai + A′i+1 )B′j+1 ; 11. A′i+1 B′j+1 ; 16 4 4 1 1 ′ 1 1 Ai ; 14. (5 + 4B′j+1 )Ai ; (Ai+1 + Bj + 5Ai + 5B′j+1 + 4Ai B′j+1 ); 15. 13. 12. A′i+1 B′j+1 ; 4 16 16 16 1 1 ′ 1 17. 16. (5A′i+1 + Bj + Ai + 5B′j+1 + 4A′i+1 B′j+1 ); (5 + 4B′j+1 )A′i+1 ; 18. Ai+1 ; 16 16 16 1 1 1 1 20. (2 + Bj + B′j+1 )Ai ; 19. Ai ; 21. (A′i+1 + Bj + B′j+1 + Ai ); 22. (2 + Bj + B′j+1 )A′i+1 ; 4 4 4 4 1 1 1 ′ 1 24. Ai ; 25. (5 + 4Bj )Ai ; 26. (Ai+1 + 5Bj + 5Ai + B′j+1 + 4Ai Bj ); 23. A′i+1 ; 4 16 16 16 1 1 1 ′ 27. (5A′i+1 + 5Bj + Ai + B′j+1 + 4A′i+1 B′j+1 ); 28. (5 + 4Bj )A′i+1 ; 29. Ai+1 ; 16 16 16 1 1 30. Ai Bj ; 31.Ai Bj ; 32. (2 + Ai + A′i+1 )Bj ; 33. A′i+1 Bj ; 4 4 1 1 34. A′i+1 Bj ; 35. (5 + 4Ai )Bj ; 4 16 1 1 1 1 1 1 36. (5 + 4A′i+1 )Bj ; 37. Ai Bj ; 38. Bj ; 39. A′i+1 Bj ; 40. Bj ; 41. Bj . 16 4 4 4 16 16 We have 8.

Ai B′j+1 ;

9. Ai B′j+1 ;

Theorem 5 ([16]). Let x−2 < x−1 < 0 = x0 < · · · < xm = 1 < xm+1 < xm+2 , and y−2 < y−1 < 0 = y0 < · · · < yn = 1 < yn+1 < yn+2 ,

10.

282

R.-H. Wang et al. / Journal of Computational and Applied Mathematics 265 (2014) 276–289

then, (1) for every −1 ≤ i ≤ m and −1 ≤ j ≤ n, Bi,j (x, y) is unique and strictly positive inside its support, and each function in (2)

¯ mn ) supported by the support of Bi,j (x, y) is a constant multiple of Bi,j (x, y). S21 (∆ (2) ∀(x, y) ∈ D = [0, 1] × [0, 1], we have m  n 

Bi,j (x, y) ≡ 1.

i=−1 j=−1

(3) ∀(x, y) ∈ D, we have m  n 

(−1)i+j hi+1 kj+1 Bi,j (x, y) = 0.

i=−1 j=−1

(4) for any i0 , j0 , − 1 ≤ i0 ≤ m, − 1 ≤ j0 ≤ n, the collection of the following B-splines

Ai0 ,j0 = {Bi,j (x, y) ∈ A : (i, j) ̸= (i0 , j0 )} (2)

¯ mn ). is a basis of S21 (∆ (2)

(2)

¯ mn ), Wmn : C (D) → S21 (∆ ¯ mn ) Considering two bivariate quasi-interpolation operators (cf. [16]) Vmn : C (D) → S21 (∆ defined as follows:   m  n  xi + xi+1 yj + yj+1 f , Bi,j (x, y),

Vmn (f ) =

2

i=−1 j=−1

Wmn (f ) =

m  n 

2

λi,j (f )Bi,j (x, y),

(3.6)

(3.7)

i=−1 j=−1

where



xi +xi+1 2

,

yj +yj+1

λi,j (f ) = 2f

2





is the center of the support of Bi,j (x, y), and

xi + xi+1 yj + yj+1 , 2 2



1

− (f (xi , yj ) + f (xi+1 , yj ) + f (xi , yj+1 ) + f (xi+1 , yj+1 )). 4

(2)

(2)

¯ mn ) becomes S21 (∆mn ), (3.6), (3.7) are changed into (3.4), (3.5) respectively. Clearly, if hi = h, kj = k, S21 (∆ Theorem 6 ([16]). Vmn (f ) ≡ f , Wmn (f ) ≡ f ,

f = 1, x, y, xy,

∀f ∈ P2 .

Denote the continuous module of f ∈ C (D) by

ω(f , ϱ) = sup{|f (x, y) − f (u, v)| : (x, y), (u, v) ∈ D, |(x, y) − (u, v)| < ϱ}, and

ϱmn = max(hi , kj ), i ,j   ′ ϱmn = max( 9h2 + k2 , h2 + 9k2 ), h = max(hi ), k = max(kj ), i j   k ∂ f ϱmn ωk f = max ω , , l=0,...,k ∂ xk−l ∂ yl 2    ∂ kf  ∥Dk f ∥ = max sup  k−l l (x, y) . l=0,...,k (x,y)∈D ∂ x ∂y For the sufficiently larger m and n (m, n ≥ N0 ), the support center of B-spline Bi,j (x, y) is located in the interior of D. Let ∥ · ∥ be the supremum over D.

R.-H. Wang et al. / Journal of Computational and Applied Mathematics 265 (2014) 276–289

283

We have the following results (cf. [19]). Theorem 7. Let f ∈ C (D), then for all m, n ≥ N0 , we have ′ ∥f − Vmn (f )∥ ≤ ω(f , ϱmn ).

If f ∈ C 1 (D), then we have

∥f − Vmn (f )∥ ≤ 2ϱmn ω1 (f ). If f ∈ C 2 (D), then we have 2 ∥f − Vmn (f )∥ ≤ ϱmn ∥D2 f ∥.

Theorem 8. Let f ∈ C (D), then for all m, n ≥ N0 , if f ∈ C 2 (D), then we have 2 ∥f − Wmn (f )∥D ≤ 2ϱmn ω2 (f ).

If f ∈ C 3 (D), then we have

∥f − Wmn (f )∥D ≤

2 3

3 ϱmn ∥D3 f ∥.

4. The computation of the approximation signed distance function In this section, we describe how we approximate the SDF of a given curve Γ by using the bivariate spline quasiinterpolant. Without loss of generality, we consider the simple closed curve Γ (t ) = (x(t ), y(t )), t ∈ [0, 2π ] lying in the unit square D = [0, 1] × [0, 1].    In order to get Vmn , Wmn by (3.6), (3.7) respectively, we need to compute u

xi +xi+1 2

,

yj +yj+1 2

,Γ

and u((xi , yj ), Γ ),

where u((x, y), Γ ), (x, y) ∈ D is the SDF of Γ (t ). Denote

 B =

xi + xi+1 yj + yj+1 , 2 2



, i = −1, . . . , m, j = −1, . . . , n



∪ {(xi , yj ), i = −1, . . . , m + 1, j = −1, . . . , n + 1}. Algorithm 1 (Closest Point Computation). Input: the number of sampling points N, the point P, the maximum iterations N1 and the threshold τ . Step 1 Sampling N points {Pi = (x( 2Nπ i ), y( 2Nπ i ))}Ni=1 on Γ (t ). Then, denote with P0 the closest point of P in {Pi }Ni=1 . Take P0 as the initial point. Step 2 Let t0 denote the parameter of P0 on the curve, and h(t ) = (Γ (t ) − P ) · Γ ′ (t ) be a scalar function. Put n = 0. h(t ) Step 3 Computing tn+1 = tn − h′ (tn ) . If |h(tn+1 )| < τ or n > N1 , stop and output P˜ = Γ (tn+1 ); Otherwise, n = n + 1. Then n return to Step [3].

˜ Output: P. After we get the closest point of every point in B , we can easily compute Vmn and Wmn . The algorithm of the computation of the approximation signed distance function consists of three steps. Algorithm 2 (The Computation of the Approximation Signed Distance Function). Input: original partition of D



x0 = −

2 m

1

, − , 0, . . . , 1, m

m+1 m+2 m

,

m



,



y0 = −

2 m

(2)

1

, − , 0, . . . , 1, m

n+1 n+2 n

,

n



,

the corresponding spline space S 0 := S21 (∆mn ), the curve Γ (t ), the threshold σ , ρ and the maximum number of cells M. For k = 0, 1, . . .. Step 1 Using the space S k and Algorithm 1, obtain Vmn (Γ ) by (3.6).

284

R.-H. Wang et al. / Journal of Computational and Applied Mathematics 265 (2014) 276–289 Table 1 The fitting errors of Γ1 (t ) using uniform splines. k

m

n

ek (V )

ek (W )

0 1 2 3

4 8 16 32

4 8 16 32

5.46 × 10−2 2.72 × 10−2 1.34 × 10−2 6.45 × 10−3

2.12 × 10−2 1.03 × 10−2 5.05 × 10−3 2.66 × 10−3

Step 2 Computing the fitting error of Vmn (Γ ), choose the sample point Pik,j as the center point in every cell [xki , xki+1 ] × [ykj , ykj+1 ]. That is,

 Pik,j

=

k k xki + xki+1 yj + yj+1 , 2 2

 ,

i = 0, . . . , m − 1, j = 0, . . . , n − 1.

The error in every cell is ek (V )i,j = |Vmn (Γ ) − u|P k , i,j

where u is the approximation function obtained from the Newton’s method. The global fitting error is denoted as ek (V ) = max ek (V )i,j . i ,j

Step 3 If e (V )i,j < σ or min{xki+1 − xki , ykj+1 − ykj } < ρ holds for every cell [xki , xki+1 ] × [ykj , ykj+1 ] or m × n > M, stop and output Vmn (Γ ); Otherwise, for the uniform partition, make m → 2m, n → 2n; For the non-uniform case, insert k

x k +x k

yk +yk

knots i 2 i+1 and j 2 j+1 into xk , yk respectively. Denote by xk+1 , yk+1 the final knot sequence respectively. Then we obtain the corresponding S k+1 , and return to Step [1]. Output: Vmn (Γ ). We implement the algorithm for the quasi-interpolant Wmn (Γ ) in the same way. As an important application, we deal with the problem of finding the offset curve of the curve Γ . For a given curve Γ , the interior δ -offset curve Γδ is the locus of point in Γ + whose distance to Γ is δ , i.e.,

Γδ = {p ∈ Γ + | dist(p, Γ ) = δ}. The exterior offset curve can be defined in the same way with a negative distance to Γ . After we get the approximation of the SDF of the given curve, we can easily compute an implicit representation of the approximate trimmed offset curve. More specifically, let s(x, y) be the resulting bivariate spline, the trimmed offset curve of Γ can be approximately represented as a level set of the scalar function s(x, y),

Γδ = {(x, y) ∈ R2 | s(x, y) = δ}. 5. Examples In this section, we demonstrate the performance of our method by some representative examples. We choose N = 40, σ = 2 × 10−3 , ρ = 2 × 10−4 , τ = 1 × 10−5 , M = 1 × 103 , N1 = 5 and δ = 0.01. The given curve and the trimmed offsets of the curve are shown in blue color and red color respectively. The fitting errors of the given curves are given in Tables 1–3 respectively. (2) Firstly, we give some examples using the uniform spline S21 (∆mn ). 5 Example 1. The closed curve is Γ1 (t ) = ( 21 + 49 cos t , 21 + 11 sin t ), t ∈ [0, 2π ]. Fig. 5 shows the given curve Γ1 (t ) (Fig. 5(a)), the SDF of Vmn (Γ1 ) (Fig. 5(b)) and the SDF of Wmn (Γ1 ) (Fig. 5(c)) which are obtained by Algorithm 2. The corresponding trimmed offsets of Γ1 (t ) are shown in Fig. 6. 1 1 (3 + 2 cos t ) cos t , 25 + 10 (3 + 2 cos t ) sin t ), t ∈ [0, 2π ]. The given curve Γ2 (t ) Example 2. The test curve is Γ2 (t ) = ( 15 + 10 (Fig. 7(a)), the SDF of Vmn (Γ2 ) (Fig. 7(b)) and the SDF of Wmn (Γ2 ) (Fig. 7(c)) are shown in Fig. 7. Fig. 8 shows the corresponding trimmed offsets of Γ2 (t ). 9 9 Example 3. The closed curve is Γ3 (t ) = ( 20 + 12 cos 3t cos t , 20 + 12 cos 3t sin t ), t ∈ [0, 2π ], which has a self-intersection Γ3 (0). Fig. 9 shows the given curve Γ3 (t ) (Fig. 9(a)), the SDF of Vmn (Γ3 ) (Fig. 9(b)) and the SDF of Wmn (Γ3 ) (Fig. 9(c)). The trimmed offsets of Γ3 (t ) are given in Fig. 10.

R.-H. Wang et al. / Journal of Computational and Applied Mathematics 265 (2014) 276–289

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

–0.1

–0.1

–0.2

–0.2 0

0 0.2

0.2 0.4

y

0.6

0.6 0.8

0

0 0.2

0.2

0.4

0.4

x

y

0.8

0.4 0.6

0.6 0.8

1 1

(a) Γ1 (t ).

285

x

0.8 1 1

(b) Vmn (Γ1 ).

(c) Wmn (Γ1 ).

Fig. 5. Γ1 (t ) and the approximate SDFs of Γ1 (t ).

(a) Wmn (Γ1 ) + 0.01 = 0.

(b) Vmn (Γ1 ) + 0.01 = 0.

(c) Wmn (Γ1 ) − 0.01 = 0.

(d) Vmn (Γ1 ) − 0.01 = 0.

Fig. 6. The trimmed offsets of Γ1 (t ).

Table 2 The fitting errors of Γ2 (t ) using uniform splines. k

m

n

ek (V )

ek (W )

0 1 2 3

4 8 16 32

4 8 16 32

1.01 × 10−1 3.47 × 10−2 2.05 × 10−2 1.01 × 10−3

2.53 × 10−2 1.21 × 10−2 4.31 × 10−3 2.84 × 10−3

(2)

The main disadvantage of approximation by S21 (∆mn ) is that it cannot refine the partition locally. So we can use the non(2)

¯ mn to approximate the SDF adaptively and refine the partition locally. We give some examples. The uniform spline over ∆ fitting errors of the given curves are given in Tables 4–7. The resulting non-uniform type-2 triangulations are shown in Figs. 11 and 12.

286

R.-H. Wang et al. / Journal of Computational and Applied Mathematics 265 (2014) 276–289

0.2 0.1 0 –0.1 –0.2 –0.3 –0.4

0.2 0.1 0 –0.1 –0.2 –0.3 –0.4

0 0.2

0

0.2

0.6

0.4

y

x

0.2 0.4

y

0.8

0.6

0 0.2

0.4

0

0.8

0.8 1 1

(a) Γ2 (t ).

0.4 0.6

0.6

x

0.8 1 1

(b) Vmn (Γ2 ).

(c) Wmn (Γ2 ).

Fig. 7. Γ2 (t ) and the approximate SDFs of Γ2 (t ).

(a) Wmn (Γ2 ) + 0.01 = 0.

(b) Vmn (Γ2 ) + 0.01 = 0.

(c) Wmn (Γ2 ) − 0.01 = 0.

(d) Vmn (Γ2 ) − 0.01 = 0.

Fig. 8. The trimmed offsets of Γ2 (t ). Table 3 The fitting errors of Γ3 (t ) by uniform splines. k 0 1 2 3

m 5 10 20 40

n

ek (V )

ek (W )

5 10 20 40

7.60 × 10 3.25 × 10−2 1.54 × 10−2 8.50 × 10−3 −2

2.64 × 10−2 1.19 × 10−2 7.33 × 10−3 3.61 × 10−3

(2)

5 ¯ mn ) are given Example 4. The closed curve is Γ1 (t ) = ( 21 + 49 cos t , 21 + 11 sin t ). The fitting errors of the Γ1 (t ) using S21 (∆ in Tables 4 and 5. Fig. 11 shows the resulting non-uniform type-2 triangulations of Γ1 (t ).

Example 5. The test curve is Γ2 (t ) = ( 15 + (2)

1 10

(3 + 2 cos t ) cos t , 25 +

1 10

(3 + 2 cos t ) sin t ). The fitting errors of the Γ2 (t )

¯ mn ) are given in Tables 6 and 7. The resulting non-uniform type-2 triangulations of Γ2 (t ) are shown in Fig. 12. using S21 (∆

R.-H. Wang et al. / Journal of Computational and Applied Mathematics 265 (2014) 276–289

0

0

–0.1

–0.1

–0.2

–0.2

–0.3

–0.3 –0.4

–0.4 0

0 0.2 0.6

y 0.6 0.8

0

0 0.2

0.2 0.4

0.2

0.4

0.4

x

y

0.4 0.6

0.6 0.8

0.8 1 1

(a) Γ3 (t ).

287

x

0.8 1 1

(b) Vmn (Γ3 ).

(c) Wmn (Γ3 ).

Fig. 9. Γ3 (t ) and the approximate SDFs of Γ3 (t ).

(a) Wmn (Γ3 ) + 0.01 = 0.

(b) Vmn (Γ3 ) + 0.01 = 0.

(c) Wmn (Γ3 ) − 0.01 = 0.

(d) Vmn (Γ3 ) − 0.01 = 0.

Fig. 10. The trimmed offsets of Γ3 (t ). Table 4 (2 )

¯ mn ). The fitting errors of Γ1 (t ) using Wmn of S21 (∆ k

m

n

e (W )

0 1 2 3

10 12 14 16

10 12 14 16

8.19 × 10−3 4.05 × 10−3 2.26 × 10−3 1.43 × 10−3

k

6. Conclusions In this paper, a method to compute the approximation of the signed distance function is presented, which is based on (2) 2) ¯ (mn ). This method is simple since we do the quasi-interpolation operator on the bivariate spline space S21 (∆mn ) and S21 (∆ not need to compute the derivatives. The approximate function obtained by our method has the advantage of low degree and high smoothness. We can approximate the signed distance function using non-uniform spline space with the adaptive

288

R.-H. Wang et al. / Journal of Computational and Applied Mathematics 265 (2014) 276–289

(2 )

(2 )

¯ mn of Γ1 (t ) using Wmn . (a) The resulting ∆

¯ mn of Γ1 (t ) using Vmn . (b) The resulting ∆

Fig. 11. The resulting non-uniform type-2 triangulations of Γ1 (t ).

Table 5 (2 )

¯ mn ). The fitting errors of Γ1 (t ) using Vmn of S21 (∆ k

m

n

e k (V )

0 1 2 3 4

10 20 26 28 30

10 20 26 30 32

2.17 × 10−2 1.07 × 10−2 5.04 × 10−3 2.24 × 10−3 1.07 × 10−3

Table 6 (2)

¯ mn ). The fitting errors of Γ2 (t ) using Wmn of S21 (∆ k

m

n

e (W )

0 1 2 3 4

10 14 17 18 19

10 14 20 22 24

8.53 × 10−3 5.05 × 10−3 2.87 × 10−3 2.04 × 10−3 1.45 × 10−3

k

Table 7 (2 )

¯ mn ). The fitting errors of Γ2 (t ) using Vmn of S21 (∆ k

m

n

e k (V )

0 1 2 3

10 20 27 32

10 20 28 38

2.12 × 10−2 1.39 × 10−2 6.03 × 10−3 3.50 × 10−3

refinement. To compute an implicit approximation of the trimmed offset curves of the given curve is also studied. The (2) (1) proposed method can be easily generalized to S42 (∆mn ) and S31 (∆mn ) [16]. As the future work, we will use the trivariate 3 spline to approximate the SDF of a closed surface in R [20]. Acknowledgments This work is partly supported by the National Natural Science Foundation of China (Nos. 11271060, U0935004, U1135003, 211071031, 11290143, 11301052), the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State

R.-H. Wang et al. / Journal of Computational and Applied Mathematics 265 (2014) 276–289

(2)

¯ mn of Γ2 (t ) using Wmn . (a) The resulting ∆

289

(2 )

¯ mn of Γ2 (t ) using Vmn . (b) The resulting ∆

Fig. 12. The resulting non-uniform type-2 triangulations of Γ2 (t ).

Education Ministry, the Fundamental Research Funds for the Central Universities (Nos. DUT13LK07, DUT13LK45), and National Engineering Research Center of Digital Life, Guangzhou 510006, China. References [1] X. Song, B. Jüttler, A. Poteaus, Hierachical spline approximation of the signed distance function, in: Proceedings of the Shape Modeling International Conference 2010, IEEE Computer Society, 2010, pp. 241–245. [2] C.L. Bajaj, F. Bernardini, G. Xu, Automatic reconstruction of surfaces and scalar fields from 3D scans, in: Proc. SIGGRAPH, ACM, New York, USA, 1995, pp. 109–118. [3] J.D. Boissonnat, F. Cazals, Smooth surface reconstruction via natural interpolation of distance functions, in: Proc. SoCG, ACM, New York, USA, 2000, pp. 223–232. [4] G. Russo, P. Smereka, A remark on computing distance functions, Journal of Computational Physics 163 (2000) 51–67. [5] M. Sussman, P. Smereka, S. Osher, A level set method for computing solutions to incompressible two phase flow, Journal of Computational Physics 119 (1994) 146–159. [6] T. Belytschko, W.J.T. Daniel, G. Ventura, A monolithic smoothing-gap algorithm for contact-impact based on the signed distance function, International Journal for Numerical Methods in Engineering 55 (2002) 101–125. [7] J. Gomes, O.D. Faugeras, Reconciling distance functions and level sets, in: Proc. SCALE-SPACE, Springer, London, 1999, pp. 70–81. [8] S. Osher, R.P. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, Springer, 2002. [9] H.K. Zhao, T. Chan, B. Merriman, S. Osher, A variational level set approach to multiphase motion, Journal of Computational Physics 127 (1) (1996) 179–195. [10] X. Song, B. Jüttler, A. Poteaus, Medial axis computation using a hierarchical spline approximation of the signed distance function, http://www. industrial-geometry.at/uploads/nrn_report_105.pdf. [11] J. Deng, F. Chen, X. Li, C. Hu, W. Tong, Z. Yang, Y. Feng, Polynomial splines over heiarchical T -meshes, Graphical Models 70 (4) (2008) 76–86. [12] H.I. Choi, C.Y. Han, H.P. Moon, K.H. Roh, N. Wee, Medial axis transform and offset curves by Minkowski Pythagorean hodograph curves, ComputerAided Design 31 (1999) 59–72. [13] T. Maekawa, An overview of offset curves and surfaces, Computer-Aided Design 31 (1999) 165–173. [14] S. Patty, Level set framework, signed distance function and various tools, 2010. http://impact.byu.edu/Image%20Processing%20Seminar/W02_ SignedDistanceFunction.pdf. [15] C. Dapogny, P. Frey, Computation of the signed distance function to a discrete contour on adapted triangulation, Calcolo 49 (3) (2012) 193–219. [16] R.H. Wang, Mutivariate Spline Functions and their Applications, Science Press, Kluwer Pub., Beijing, New York, 2001. [17] R.H. Wang, The structural characterization and interpolation for multivariate splines, Acta Mathematica Sinica 18 (1975) 91–106. [18] R.H. Wang, The dimension and basis of spaces of mutivariate splines, Journal of Computationl and Applied Mathematics 12–13 (1985) 163–177. [19] C.J. Li, Multivariate splines on special triangulations and their applications, Ph.D. Dissertation, Dalian university of technology, 2004. [20] R.H. Wang, X.Q. Shi, A kind of cubic c 1 -interpolations in the n-dimensional finite element methods, Journal of Mathematical Research and Exposition 9 (2) (1989) 173–179.