SignedLp-distance Fields - CiteSeerX
Signed Lp-distance Fields Alexander Belyaev a a
Pierre-Alain Fayolle b
Alexander Pasko c
School of Engineering & Physical Sciences, Heriot-Watt University, Edinburgh, UK b Computer Graphics Laboratory, University of Aizu, Japan c The National Centre for Computer Animation, Bournemouth University, UK
Abstract We introduce and study a family of generalized double-layer potentials which are used to build smooth and accurate approximants for the signed distance function. Given a surface, the value of an approximant at a given point is a power mean of distances from the point to the surface points parameterized by the angle they are viewed from the given point. We analyze mathematical properties of the potentials and corresponding approximants. In particular, approximation accuracy estimates are derived. Our theoretical results are supported by numerical experiments which reveal high practical potential of our approach.
1. Introduction Singular integrals are central to many mathematical and physical theories and constructions. In this paper, we use singular integrals to construct smooth approximations of distance functions. The paper is inspired by recent works [15,7,4] where remarkable properties of the mean value normalization function (the sum of the mean value weights in the case of a polygon/polyhedron) were discovered and studied. In this paper, we introduce signed Lp -distance fields, simple but useful generalizations of the mean value normalization function, and exploit their properties to build smooth approximations to the signed distance function. Distance and distance-based scalar fields are widely used in modeling, computer graphics, scientific visualization, robotics, medical imaging, and many other disciplines [16,8,28,23]. In particular, fast and reliable estimation of the distance function is important for efficient implementations of level-set methods [20, Chapter 7]. For many applications, the exact distance field is not needed and smooth approximations of the distance field are employed. For example, in the area of solid modeling, Shapiro and co-workers [2,27] discussed Ricci operations [24] and used R-functions [26] for building approximate distance fields from local approximations. Fayolle et al. [9] developed smooth alternatives to the classical min/max operations for constructing signed approximate distance functions. Very recently, Freytag et al. [11] introduced sampled smooth approximations of a distance function and used them to enhance the Kantorovich
method [18,19] for computational mechanics purposes. In image processing and computer graphics, Ahuja, Chuang, and co-workers [1,6] and Peng et al. [21,22] introduced generalized potential fields and applied them for shape skeletonisation and 3D texture modeling, respectively. It is also worth to mention works [14,13] where a diffusion-based distance function was introduced and used for static and dynamic shape analysis. Our approach is conceptually close to that of Ahuja and Chuang (and Peng et al. who re-introduced the AhujaChuang ideas for computer graphics applications) but instead of using generalized Newtonian potentials we introduce and study a family of generalized double-layer potentials and corresponding signed Lp -distance fields. There are several advantages of using generalized doublelayer potentials to compare with their single-layer counterparts. In particular, (i) our signed Lp -distance fields deliver smooth and accurate approximations of the signed distance function; (ii) exact formulas for the generalized double-layer potentials generated by triangle meshes (polylines in 2D) can be obtained; (iii) mathematical properties of the generalized doublelayer potentials and their corresponding signed Lp distance fields are easy to analyze; (iv) our generalized double-layer potentials can be considered as a generalization of the mean value normalization function [7,4] and can be used for high-order transfinite interpolation purposes. It is interesting that singular integrals similar to the double-layer potentials we deal with serve as a main build-
ing block for the so-called surface generalized Born model in biomolecular modeling [12,25]. 2. Generalized double-layer potentials Consider a smooth oriented surface S ⊂ Rn and a singular integral Z ny · (y − x) φk (x) = σ(y) dSy , k ≥ 0, (1) n+k S |x − y|
Fig. 1. Left and middle: illustrations for the definition of ϕk (x). Right: potential ϕk (x) generated by a segment ab and notations used.
p is, the better signed Lp -distance field ψp (x) approximates dist(x).
where x is a point outside S, y ∈ S, ny is the orientation normal at y, dSy is the surface element at y, and σ(y) is a density function defined on S. The classical double-layer potential corresponds to (1) with k = 0. If the electric dipoles are distributed over S with density σ, then φ0 (x) is proportional to the electric field generated by the dipoles at x. Let us assume that σ ≡ 1 (i.e., the dipoles are uniformly distributed over S) and introduce a generalized doublelayer potential by Z Z dΩy ny · (y − x) dS = , (2) ϕk (x) = y n+k k Ω |x − y| S |x − y|
3. Mathematical properties of approximate distance functions In this section we formulate our main mathematical results which describe approximation properties of (4). See Appendix A for proofs. Denote by vn and an the volume of the unit ball in Rn and area of the unit sphere, respectively. We have an = nvn , v0 = 1, v1 = 2, v2 = π, v3 = 4π/3, and an = 2π n/2 Γ(n/2), where Γ(·) is the Gamma function. Consider a domain D ⊂ Rn bounded by surface S.
where k > 0, Ω is the unit sphere centered at x, and dΩy is the solid angle at which surface element dSy is seen from x. We have used a simple relation dSy = ρn dΩy / h, where h = ny · (y − x) is the distance from x to the plane tangent to S at y. The two-dimensional analog of (2) can be written in polar coordinates as Z 2π dθ ϕk (x) = , ρ(θ) = |x − y| , (3) k ρ(θ) 0
Proposition 1 Assume that D is star-shaped w.r.t. x ∈ D. Let 1 < p < q < ∞. Then 1/p
an ψ1 (x) > (an )
1/p
,
p ≥ 1,
ψq (x) > dist(x).
Further ψk (x) → dist(x), as k → ∞. In other words, ψ∞ (x) = dist(x). In addition, if x ∈ D does not belong to the medial axis of D, then n − 1 ln k 1 ψk (x) = dist(x) 1 + +O , as k → ∞. 2 k k
where θ is the angle between vector y − x and a fixed direction. Note that the surface integral in (2) is correctly defined for an arbitrary smooth oriented surface S, while the integral over unit sphere Ω is properly defined if domain D is star-shaped w.r.t. x. In order to drop this limitation for the integral over Ω we follow a standard approach of alternating signs (see, for example, [7,4]). Consider a ray originated at x and intersecting S in m points y1 , . . . , ym . We set εi = 1 if the ray [x, yi ) arrives at yi from positive side of S, εi = −1 if the ray approaches yi from negative side of S, and εi = 0 if the ray is tangent to S at yi . Now let k us assume that |x − y| in the denominator of the integral P k over Ω in (2) means i εi |x − yi | . See the left and middle images of Fig. 1 below for a visual explanation. If S is a closed curve and k = 1 then (3) defines the normalization function corresponding to the transfinite mean value interpolation for the domain bounded by S. It is easy to see that ψp (x) = [1 /ϕp (x) ]
1/q
ψp (x) > (an )
Although this convergence is not fast, below we will see that, after a proper normalization, ψk (x) delivers a very accurate approximation of dist(x) near S even for small k. Fig. 2 illustrates that in the simplest case of n = 1. As 1/k seen in the left image, (an ) ψk (x) converges to dist(x) from above as k → ∞. However a different normalization of ψk (x) (no normalization is needed if n = 1) allows us to obtain a much better approximation dist(x) near the boundary.
(4) Fig. 2. Approximate distance functions (an )1/k ψk (x) (left) and ψk (x) (right) defined for [0, 1] (n = 1). See the main text for details.
vanishes as x approaches S. In the next section, we show that (4) delivers a smooth approximation for the signed distance function dist(x) from S and, furthermore, the bigger
Let us introduce a sequence 2
c0 = an /2,
c1 = vn−1 ,
ck+2 =
k+1 ck k+n
For an arbitrary positive integer k, the value of this integral can be evaluated symbolically (for example, one can use Maple or Mathematica for that). If k is odd, the integral can be expressed as a polynomial of degree k of variable t, where t = tan(γ/2). In particular, 1 1 + t, ϕ1 (x) = a b 3 1 1 1 1 1 1 + 3 t+ + t 1 + t2 ϕ3 (x) = 3 a3 b 6 a b
(5)
and define 1/k
Ψk (x) = [ck ]
−1/k
ψk (x) = [ϕk (x)/ck ]
.
Proposition 2 Let n be the outer normal for S = ∂D. If x ∈ ∂D then Ψk (x) = 0,
∂Ψk (x)/ ∂n = −1.
If ∂D contains a planar part and x is an internal point of that planar part, then, in addition, ∂ s Ψk (x)/ ∂ns = 0,
Note that ϕ1 (x) is the weight corresponding to ab and induced by the mean value coordinates [10,15]. If k is even, then more complex expressions are obtained for ϕk (x), except the case when k = 0: ϕ0 (x) = γ(x) is the standard double-layer potential induced by ab. As expected, in the 3D case, the situation is more complex. Let S be a triangulated surface and abc ∈ S be a mesh triangle. While ϕ0 (x) is just a solid angle at which abc is seen from x and the formula for ϕ1 (x) generated by triangle abc is presented in [17], deriving a closed-form expression for ϕk (x), k > 1, is a difficult and tedious task. Fortunately, a clear way to obtain analytical expressions for a family of singular integrals including ϕk (x) was proposed very recently in [5]. Given a polyline in 2D (triangle mesh in 3D), potential ϕk (x) at point x is obtained by summing the contributions of the polyline segments (mesh triangles).
s = 2, 3, . . . , k.
A simple asymptotic analysis of (5) yields 1 n − 1 ln k 1/k +O , as k → ∞, [ck ] =1− 2 k k which, in view of the last statement of Proposition 1, gives Ψk (x) = dist(x) + O (1/k) , as k → ∞,
(6)
for x outside the medial axis of D. Note that, according to Proposition 2, near a flat part of ∂D the convergence of Ψk (x) to dist(x) is much faster than (6). Our last theoretical result is presented below. In particular, it implies that the shapes of the graphs of ψk (x) and Ψk (x) are similar to that of dist(x). Proposition 3 We have
5. Numerical experiments
∆ϕk (x) = k(k + n)ϕk+2 (x) In Fig. 3, we provide the reader with a visual comparison of the level sets of dist(x) with those of Ψ1 , Ψ3 and Ψ5 for sufficiently complex planar polylines. We observe that the approximation quality improves quickly as k increases. Fig. 4 visualizes the relative error between dist(x) and ψk (x) and between dist(x) and Ψk (x) (normalized Lp distance), k = 1, 3, 5, for a squared domain. As before, the error goes down quickly as k increases. One can observe that the approximation error is relatively high near the corners. In Fig. 5, we investigate approximation properties of Lp distance fields near a corner. As expected, the approximation accuracy (as before, we calculate the relative error) at a small vicinity of a corner point is lower than near an internal point of segment. Further, the relative error increases as the corner angle becomes sharper. However, even for a sharp corner, the approximation accuracy quickly improves, as k increases. Finally, Fig. 6 illustrates our approach in 3D with the relative error between dist(x) and Ψ2 (x) plotted on a crosssection of the object.
and, therefore, ψk (x) has no local minima inside D. Note that the gradient and the other derivatives of ϕk (x) can be easily obtained by differentiating (2) by x as many times as needed. 4. Singular potentials for polylines and meshes Let us start from the 2D case and consider the generalized double-layer potentials generated by oriented segment ab. Given point x, γ = γ(x) is the angle between xa and xb, a and b are the distances from x to a and b, respectively, and ρ(θ) is the distance from x to point y ∈ ab. See the right image of Fig. 1 for a visual illustration of the notations used. The area of 4xab is equal to the sum of areas of 4xay and 4xyb. We have ab sin γ = aρ(θ) sin θ + bρ(θ) sin(γ − θ), which gives an explicit expression of ρ(θ) ab sin γ . ρ(θ) = a sin θ + b sin(γ − θ) Then generalized double-layer potential ϕk (x) generated by the segment is given by Z γ Z γ dθ 1 ϕk (x) = = (a sin θ+b sin(γ−θ))k dθ. k k ρ(θ) (ab sin γ) 0 0
6. Discussion and conclusion In this paper, we have introduced a family of generalized double-layer potentials and their corresponding 3
Fig. 3. Polygon approximation of two curves used for our experiments (left-most images). The polygon meshes consist in 201 and 29 segments respectively. Contour plots for: the signed distance, Ψ1 , Ψ3 and Ψ5 for a 2D polygon mesh (columns 2 to 5). Ψ5 (right-most images) provides already a good approximation of the distance function (second column). Note that the holes in the curve for the letter A are properly handled.
Fig. 6. Relative errors between dist(x) and Ψ2 (x) for a cube (left) and the rocker-arm model (right).
computation of Lp -distance fields can be significantly accelerated via a hierarchical approach: an accurate evaluation of a Lp -distance field at a given point does not require careful calculations of contributions of surface parts distant from the point. For polylines in 2D and triangle meshes in 3D, special correctors (boundary layers) should be added to the Lp -distance fields in order to achieve a better approximation accuracy at polyline/mesh vertices. We hope that our approach can lead to improving and enhancing existing applications of approximate distance functions in image processing, computer vision, computer graphics, geometric modeling, and several other branches of scientific computing [1,14,22]. In particular, we foresee potential applications of Lp -distance fields in computational mechanics where some new promising computational techniques rely on constructing accurate and efficient approximate smoothed distance functions and their derivatives [11, Section 3.2].
Fig. 4. Relative errors for a square. The first row corresponds to the relative error between the exact signed distance and ψ1 , ψ3 and ψ5 respectively. The second row uses the normalized functions: Ψ1 , Ψ3 and Ψ5 instead.
Appendix A: Proofs Fig. 5. Relative errors between the exact signed distance and Ψ5 for various opening angles between two segments. The angles vary between 0◦ (top row, left column) and 45◦ (bottom row, right column).
Proof of Proposition 1. An integral version of Jensen’s inequality (see, for example, [29, Problem 7.5]) states that Z Z F f (x)w(x) dx ≤ F (f (x)) w(x) dx,
Lp -distance fields. They possess interesting mathematical properties, and, in particular, the Lp -distance fields deliver smooth and accurate approximations of the distance function. This work is just a beginning of our study of generalized double-layer potentials and their applications. A practical
R
R
where F (t) is convex and w(x) is a positive weight function satisfying Z w(x) dx = 1. R
4
1 (k + 2)h2 h ∂ = k+2 − k+2 ∂h ρ ρ ρk+4 (see Fig. 7 for the notations used). Thus
Let 1 < p < q < ∞, F (t) = tq/p , R = [0, 2π], w ≡ 1/2π, p and f (θ) = (1/ρ(θ)) . Since F (t) is strictly convex, we have 1 1 Z 2π Z 2π dθ p dθ q 1 1 < , ψˆp (x) > ψˆq (x), 2π 0 ρ(θ)p 2π 0 ρ(θ)q
∂ϕk 1 = ϕk − (k + 2)hϕk+2 ∂h h or, equivalently,
1/k where ψˆk (x) = (2π) ψk (x). The multidimensional case n > 2 is considered similarly. Let x ∈ D is fixed and k → ∞. We have according to Laplace’s method (see, for example, [3]) Z Z 1 dΩ = exp k ln dΩ ϕk (x) = k ρ Ω Ω ρ
=
C1 + o(1) 1 , k (n−1)/2 ρkmin
(k + 2)ϕk+2 =
where ρmin = dist(x).
1/k k (n−1)/2 = dist(x) C1 + o(1) (n−1)/2 1/k = dist(x) [C2 + o(1)] k 1/k , as k → ∞, 1 ψk (x) = ϕk (x)
1/k
where C1 and C2 are positive constants. It remains to note that n−1 (n − 1) ln k 1 n − 1 ln k 2 1/k = exp = 1+ +O , k 2k 2 k k
with ck defined by (5), where singular potential ϕk (x) is generated by segment ab and x approaches an internal point of ab. For the potential generated by a polygon containing segment ab, we have a slightly bigger remainder: ϕk (x) = ck hk + O(1), as h → 0,
as k → ∞. Proof of Proposition 2. First let us consider the 2D case and study an asymptotic behavior of potential ϕ1 (x) generated by a single segment ab. See the left image of Fig. 7 for notations used. Assume that h, the distance from x to ab, is small: h 1. We have γ = π − α − β, a(h) = a + O(h2 ),
α = arctan (h/a) = h/a + O(h2 ),
b(h) = b + O(h2 ),
β = arctan (h/b) = h/b + O(h2 ),
γ 2ab 1 tan = + O(h), 2 a+bh
as the contribution of the other segments of the polygon is O(1). Thus 1/k ck Ψk (x) = = h + O hk+1 as h → 0. ϕk (x) In particular, if x is an internal point of ab, we have ∂ s Ψk (x)/ ∂ns = 0,
2 ϕ1 (x) = + O(h). h
In the 2D case, for a single straight segment, we have Z Z dθ h dsy ϕk (x) = = , ρ = |x − y|. ρ(θ)k ρ(θ)k+2 ∂ρ h = , ∂h ρ
∂ ∂h
1 ρk+2
=−
s = 2, 3, . . . , k.
Now let us consider a domain D with a smooth boundary S. Without loss of generality we can assume that D is convex (the case of non-convex domains can be treated similarly to how it is done in [7]). Consider a point x ∈ D on a distance h 1 from S. Let ab be a circular arc osculating S at the point closest to x such that the segment ab is perpendicular to the direction from x to its closest neighbor on S, as seen in the right image of Fig. 7. Let us also consider another circular arc cd such that it is situated inside D, ab and cd are collinear, and the direction from x to its closest neighbor on S is the bisector for cd. Now we approximate ϕk (x) given by (3) by the sum of integrals of 1/ρ(θ)k over the circular arcs ab and cd. Note that, as x approaches S, the osculating arc ab delivers a better and better local approximation of S and the leading term in asymptotics ϕk (x), as h → 0, is the same as that in the integral of 1/ρ(θ)k over the circular arc ab. The law of cosines yields q ρ(θ) = −(R − h) cos θ + R2 − (R − h)2 sin2 θ.
Fig. 7. Notations used for asymptotic analysis of potentials ϕk (x) generated by a segment (left) and a circular arc (right).
Note that p ρ = h2 + z 2 ,
(7)
Note that ϕk (x) is odd w.r.t. h. Thus the expansion of ϕk (x) w.r.t. h 1 contains only odd degrees of h. For example, 2 ϕ1 (x) = + a1 h + a3 h3 + a5 h5 + O h7 , h → 0. h Then (7) implies that 4 1 − 2a3 h − 4a5 h3 + O h5 ϕ3 (x) = 3 3 h and we can observe that (7) kills the linear w.r.t. h term in the expansion of ϕk (x). Now, applying the mathematical induction, we arrive at ϕk (x) = ck hk + O(h)
Thus
1 1 ∂ϕk ϕk − . h2 h ∂h
(k + 2)h , ρk+4 5
where R is the radius of the osculating circle. We can safely assume that R = 1. Our task now is to study the asymptotic behavior of
[9] P.-A. Fayolle, A. Pasko, B. Schmitt, and N. Mirenkov. Constructive heterogeneous object modeling using signed approximate real distance functions. Journal of Computing and Information Science in Engineering, ASME Transactions, 6(3):221–229, 2006. [10] M. S. Floater. Mean value coordinates. Computer Aided Geometric Design, 20(1):19–27, 2003. [11] M. Freytag, V. Shapiro, and I. Tsukanov. Finite element analysis in situ. Finite Elements in Analysis and Design, 47(9):957–972, 2011. [12] A. Ghosh, C. S. Rapp, and R. A. Friesner. Generalized Born model based on a surface integral formulation. J. Phys. Chem. B, 102:10983–10990, 1998. [13] L. Gorelick, M. Blank, E. Shechtman, M. Irani, and R. Basri. Actions as space-time shapes. IEEE Transactions on Pattern Analysis and Machine Intelligence, 29(12):2247–2253, 2007. [14] L. Gorelick, M. Galun, E. Sharon, R. Basri, and A. Brandt. Shape representation and classification using the Poisson equation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 28(12):1991–2005, 2006. [15] K. Hormann and M. S. Floater. Mean value coordinates for arbitrary planar polygons. ACM Transactions on Graphics, 25(4):1424–1441, 2006. [16] M. W. Jones, J. A. Baerentzen, and M. Sramek. 3D distance fields: a survey of techniques and applications. IEEE Transactions on Visualization and Computer Graphics, 12(4):581–599, 2006. [17] T. Ju, S. Schaefer, and J. Warren. Mean value coordinates for closed triangular meshes. ACM Transactions on Graphics, 24(3):561–566, 2005. ACM SIGGRAPH 2005 issue. [18] L. V. Kantorovich. Some remarks on Ritz’s method. Trydy vysshego voenno-morskogo inzhenerno-stroitel’nogo uchilishcha, (3), Leningrad, 1941. (Russian). [19] L. V. Kantorovich and V. I. Krylov. Approximate Methods of Higher Analysis. Interscience Publishers, 1958. Chapter 4, Paragraph 2. [20] S. J. Osher and R. P. Fedkiw. Level Set Methods and Dynamic Implicit Surfaces. Springer, 2002. [21] J. Peng. Thick Surfaces: Interactive Modeling of Topologically Complex Geometric Details. PhD, Department of Computer Science, New York University, 2004. [22] J. Peng, D. Kristjansson, and D. Zorin. Interactive modeling of topologically complex geometric detail. ACM Trans. Graph., 23:635–643, 2004. ACM SIGGRAPH 2004 issue. [23] T. Reiner, G. M¨ uckl, and C. Dachsbacher. Interactive modeling of implicit surfaces using a direct visualization approach with signed distance functions. Computers & Graphics, 35:596–603, 2011. SMI 2011 issue. [24] A. Ricci. A constructive geometry for computer graphics. The Computer Journal, 16(2), 1973. [25] A. N. Romanov, S. N. Jabin, Y. B. Martynov, A. V. Sulimov, F. V. Grigoriev, and V. B. Sulimov. Surface generalized Born method: A simple, fast, and precise implicit solvent model beyond the Coulomb approximation. J. Phys. Chem. A, 108(43):9323–9327, 2004. [26] V. L. Rvachev. Theory of R-functions and Some Applications. Naukova Dumka, 1982. Russian. [27] V. Shapiro. Semi-analytic geometry with R-functions. Acta Numerica, 16:239–303, 2007. [28] K. Siddiqi and S. Pizer. Medial Representations: Mathematics, Algorithms and Applications. Springer, 2008. [29] J. M. Steele. The Cauchy-Schwarz Master Class: An Introduction to the Art of Mathematical Inequalities. Cambridge University Press, 2004.
−k Zπ/2 q 2 2 (h − 1) cos θ + 1 − (1 − h) sin θ Ik (h) = dθ. −π/2
It is easy to verify that Ik (h) =
ck +O hk
1 hk−1
, as h → 0.
In the nD case, for a single (n − 1)-dimensional mesh tetrahedron (a mesh triangle in 3D), we have Z Z h dSy dΩy = , ρ = |x − y|. ϕk (x) = k ρ ρk+n So we arrive at a generalization of (7) for the multidimensional case 1 1 ∂ϕk (x) (k + n)ϕk+2 (x) = 2 ϕk (x) − . h h ∂h It is shown in [4] that vn−1 + O(1) ϕ1 (x) = h which justifies that c1 = vn−1 in (5). Proof of Proposition 3. Differentiating of (1) w.r.t x yields ∆φk (x) = k(k + n)φk+2 (x). For k = 1 and σ(y) ≡ 1 this formula coincides with that obtained in [4]. Obviously φk (x) is positive if surface ”dipole” density σ(y) is positive. Assume that φk (x) has a local maximum at some point inside the domain bounded by S. Then ∆φk (x) ≤ 0 at that point. We arrive at a contradiction which completes the proof. References [1] N. Ahuja and J.-H. Chuang. Shape representation using a generalized potential field model. IEEE Transactions on Pattern Analysis and Machine Intelligence, 19(2):169–176, 1997. [2] A. Biswas and V. Shapiro. Approximate distance fields with non-vanishing gradients. Graphical Models, 66(3):133–159, 2004. [3] N. Bleistein and R. A. Handelsman. Asymptotic Expansions of Integrals. Dover, 1986. [4] S. Bruvoll and M. S. Floater. Transfinite mean value interpolation in general dimension. J. Comp. Appl. Math., 233:1631–1639, 2010. [5] M. Carley. Potential integrals on triangles. Technical Report arXiv1201.4938, January 2012. [6] J.-H. Chuang, C.-H. Tsai, and M.-C. Ko. Skeletonization of three-dimensional object using generalized potential field. IEEE Transactions on Pattern Analysis and Machine Intelligence, 22(11):1241–1251, 2000. [7] C. Dyken and M. S. Floater. Transfinite mean value interpolation. Comp. Aided Geom. Design., 26:117–134, 2009. [8] R. Fabbri, L. da F. Costa, J. C. Torelli, and O. M. Bruno. 2D Euclidean distance transform algorithms: A comparative survey. ACM Computing Surveys, 40(1):1–44, 2008.
6
Pierre-Alain Fayolle b
Alexander Pasko c
School of Engineering & Physical Sciences, Heriot-Watt University, Edinburgh, UK b Computer Graphics Laboratory, University of Aizu, Japan c The National Centre for Computer Animation, Bournemouth University, UK
Abstract We introduce and study a family of generalized double-layer potentials which are used to build smooth and accurate approximants for the signed distance function. Given a surface, the value of an approximant at a given point is a power mean of distances from the point to the surface points parameterized by the angle they are viewed from the given point. We analyze mathematical properties of the potentials and corresponding approximants. In particular, approximation accuracy estimates are derived. Our theoretical results are supported by numerical experiments which reveal high practical potential of our approach.
1. Introduction Singular integrals are central to many mathematical and physical theories and constructions. In this paper, we use singular integrals to construct smooth approximations of distance functions. The paper is inspired by recent works [15,7,4] where remarkable properties of the mean value normalization function (the sum of the mean value weights in the case of a polygon/polyhedron) were discovered and studied. In this paper, we introduce signed Lp -distance fields, simple but useful generalizations of the mean value normalization function, and exploit their properties to build smooth approximations to the signed distance function. Distance and distance-based scalar fields are widely used in modeling, computer graphics, scientific visualization, robotics, medical imaging, and many other disciplines [16,8,28,23]. In particular, fast and reliable estimation of the distance function is important for efficient implementations of level-set methods [20, Chapter 7]. For many applications, the exact distance field is not needed and smooth approximations of the distance field are employed. For example, in the area of solid modeling, Shapiro and co-workers [2,27] discussed Ricci operations [24] and used R-functions [26] for building approximate distance fields from local approximations. Fayolle et al. [9] developed smooth alternatives to the classical min/max operations for constructing signed approximate distance functions. Very recently, Freytag et al. [11] introduced sampled smooth approximations of a distance function and used them to enhance the Kantorovich
method [18,19] for computational mechanics purposes. In image processing and computer graphics, Ahuja, Chuang, and co-workers [1,6] and Peng et al. [21,22] introduced generalized potential fields and applied them for shape skeletonisation and 3D texture modeling, respectively. It is also worth to mention works [14,13] where a diffusion-based distance function was introduced and used for static and dynamic shape analysis. Our approach is conceptually close to that of Ahuja and Chuang (and Peng et al. who re-introduced the AhujaChuang ideas for computer graphics applications) but instead of using generalized Newtonian potentials we introduce and study a family of generalized double-layer potentials and corresponding signed Lp -distance fields. There are several advantages of using generalized doublelayer potentials to compare with their single-layer counterparts. In particular, (i) our signed Lp -distance fields deliver smooth and accurate approximations of the signed distance function; (ii) exact formulas for the generalized double-layer potentials generated by triangle meshes (polylines in 2D) can be obtained; (iii) mathematical properties of the generalized doublelayer potentials and their corresponding signed Lp distance fields are easy to analyze; (iv) our generalized double-layer potentials can be considered as a generalization of the mean value normalization function [7,4] and can be used for high-order transfinite interpolation purposes. It is interesting that singular integrals similar to the double-layer potentials we deal with serve as a main build-
ing block for the so-called surface generalized Born model in biomolecular modeling [12,25]. 2. Generalized double-layer potentials Consider a smooth oriented surface S ⊂ Rn and a singular integral Z ny · (y − x) φk (x) = σ(y) dSy , k ≥ 0, (1) n+k S |x − y|
Fig. 1. Left and middle: illustrations for the definition of ϕk (x). Right: potential ϕk (x) generated by a segment ab and notations used.
p is, the better signed Lp -distance field ψp (x) approximates dist(x).
where x is a point outside S, y ∈ S, ny is the orientation normal at y, dSy is the surface element at y, and σ(y) is a density function defined on S. The classical double-layer potential corresponds to (1) with k = 0. If the electric dipoles are distributed over S with density σ, then φ0 (x) is proportional to the electric field generated by the dipoles at x. Let us assume that σ ≡ 1 (i.e., the dipoles are uniformly distributed over S) and introduce a generalized doublelayer potential by Z Z dΩy ny · (y − x) dS = , (2) ϕk (x) = y n+k k Ω |x − y| S |x − y|
3. Mathematical properties of approximate distance functions In this section we formulate our main mathematical results which describe approximation properties of (4). See Appendix A for proofs. Denote by vn and an the volume of the unit ball in Rn and area of the unit sphere, respectively. We have an = nvn , v0 = 1, v1 = 2, v2 = π, v3 = 4π/3, and an = 2π n/2 Γ(n/2), where Γ(·) is the Gamma function. Consider a domain D ⊂ Rn bounded by surface S.
where k > 0, Ω is the unit sphere centered at x, and dΩy is the solid angle at which surface element dSy is seen from x. We have used a simple relation dSy = ρn dΩy / h, where h = ny · (y − x) is the distance from x to the plane tangent to S at y. The two-dimensional analog of (2) can be written in polar coordinates as Z 2π dθ ϕk (x) = , ρ(θ) = |x − y| , (3) k ρ(θ) 0
Proposition 1 Assume that D is star-shaped w.r.t. x ∈ D. Let 1 < p < q < ∞. Then 1/p
an ψ1 (x) > (an )
1/p
,
p ≥ 1,
ψq (x) > dist(x).
Further ψk (x) → dist(x), as k → ∞. In other words, ψ∞ (x) = dist(x). In addition, if x ∈ D does not belong to the medial axis of D, then n − 1 ln k 1 ψk (x) = dist(x) 1 + +O , as k → ∞. 2 k k
where θ is the angle between vector y − x and a fixed direction. Note that the surface integral in (2) is correctly defined for an arbitrary smooth oriented surface S, while the integral over unit sphere Ω is properly defined if domain D is star-shaped w.r.t. x. In order to drop this limitation for the integral over Ω we follow a standard approach of alternating signs (see, for example, [7,4]). Consider a ray originated at x and intersecting S in m points y1 , . . . , ym . We set εi = 1 if the ray [x, yi ) arrives at yi from positive side of S, εi = −1 if the ray approaches yi from negative side of S, and εi = 0 if the ray is tangent to S at yi . Now let k us assume that |x − y| in the denominator of the integral P k over Ω in (2) means i εi |x − yi | . See the left and middle images of Fig. 1 below for a visual explanation. If S is a closed curve and k = 1 then (3) defines the normalization function corresponding to the transfinite mean value interpolation for the domain bounded by S. It is easy to see that ψp (x) = [1 /ϕp (x) ]
1/q
ψp (x) > (an )
Although this convergence is not fast, below we will see that, after a proper normalization, ψk (x) delivers a very accurate approximation of dist(x) near S even for small k. Fig. 2 illustrates that in the simplest case of n = 1. As 1/k seen in the left image, (an ) ψk (x) converges to dist(x) from above as k → ∞. However a different normalization of ψk (x) (no normalization is needed if n = 1) allows us to obtain a much better approximation dist(x) near the boundary.
(4) Fig. 2. Approximate distance functions (an )1/k ψk (x) (left) and ψk (x) (right) defined for [0, 1] (n = 1). See the main text for details.
vanishes as x approaches S. In the next section, we show that (4) delivers a smooth approximation for the signed distance function dist(x) from S and, furthermore, the bigger
Let us introduce a sequence 2
c0 = an /2,
c1 = vn−1 ,
ck+2 =
k+1 ck k+n
For an arbitrary positive integer k, the value of this integral can be evaluated symbolically (for example, one can use Maple or Mathematica for that). If k is odd, the integral can be expressed as a polynomial of degree k of variable t, where t = tan(γ/2). In particular, 1 1 + t, ϕ1 (x) = a b 3 1 1 1 1 1 1 + 3 t+ + t 1 + t2 ϕ3 (x) = 3 a3 b 6 a b
(5)
and define 1/k
Ψk (x) = [ck ]
−1/k
ψk (x) = [ϕk (x)/ck ]
.
Proposition 2 Let n be the outer normal for S = ∂D. If x ∈ ∂D then Ψk (x) = 0,
∂Ψk (x)/ ∂n = −1.
If ∂D contains a planar part and x is an internal point of that planar part, then, in addition, ∂ s Ψk (x)/ ∂ns = 0,
Note that ϕ1 (x) is the weight corresponding to ab and induced by the mean value coordinates [10,15]. If k is even, then more complex expressions are obtained for ϕk (x), except the case when k = 0: ϕ0 (x) = γ(x) is the standard double-layer potential induced by ab. As expected, in the 3D case, the situation is more complex. Let S be a triangulated surface and abc ∈ S be a mesh triangle. While ϕ0 (x) is just a solid angle at which abc is seen from x and the formula for ϕ1 (x) generated by triangle abc is presented in [17], deriving a closed-form expression for ϕk (x), k > 1, is a difficult and tedious task. Fortunately, a clear way to obtain analytical expressions for a family of singular integrals including ϕk (x) was proposed very recently in [5]. Given a polyline in 2D (triangle mesh in 3D), potential ϕk (x) at point x is obtained by summing the contributions of the polyline segments (mesh triangles).
s = 2, 3, . . . , k.
A simple asymptotic analysis of (5) yields 1 n − 1 ln k 1/k +O , as k → ∞, [ck ] =1− 2 k k which, in view of the last statement of Proposition 1, gives Ψk (x) = dist(x) + O (1/k) , as k → ∞,
(6)
for x outside the medial axis of D. Note that, according to Proposition 2, near a flat part of ∂D the convergence of Ψk (x) to dist(x) is much faster than (6). Our last theoretical result is presented below. In particular, it implies that the shapes of the graphs of ψk (x) and Ψk (x) are similar to that of dist(x). Proposition 3 We have
5. Numerical experiments
∆ϕk (x) = k(k + n)ϕk+2 (x) In Fig. 3, we provide the reader with a visual comparison of the level sets of dist(x) with those of Ψ1 , Ψ3 and Ψ5 for sufficiently complex planar polylines. We observe that the approximation quality improves quickly as k increases. Fig. 4 visualizes the relative error between dist(x) and ψk (x) and between dist(x) and Ψk (x) (normalized Lp distance), k = 1, 3, 5, for a squared domain. As before, the error goes down quickly as k increases. One can observe that the approximation error is relatively high near the corners. In Fig. 5, we investigate approximation properties of Lp distance fields near a corner. As expected, the approximation accuracy (as before, we calculate the relative error) at a small vicinity of a corner point is lower than near an internal point of segment. Further, the relative error increases as the corner angle becomes sharper. However, even for a sharp corner, the approximation accuracy quickly improves, as k increases. Finally, Fig. 6 illustrates our approach in 3D with the relative error between dist(x) and Ψ2 (x) plotted on a crosssection of the object.
and, therefore, ψk (x) has no local minima inside D. Note that the gradient and the other derivatives of ϕk (x) can be easily obtained by differentiating (2) by x as many times as needed. 4. Singular potentials for polylines and meshes Let us start from the 2D case and consider the generalized double-layer potentials generated by oriented segment ab. Given point x, γ = γ(x) is the angle between xa and xb, a and b are the distances from x to a and b, respectively, and ρ(θ) is the distance from x to point y ∈ ab. See the right image of Fig. 1 for a visual illustration of the notations used. The area of 4xab is equal to the sum of areas of 4xay and 4xyb. We have ab sin γ = aρ(θ) sin θ + bρ(θ) sin(γ − θ), which gives an explicit expression of ρ(θ) ab sin γ . ρ(θ) = a sin θ + b sin(γ − θ) Then generalized double-layer potential ϕk (x) generated by the segment is given by Z γ Z γ dθ 1 ϕk (x) = = (a sin θ+b sin(γ−θ))k dθ. k k ρ(θ) (ab sin γ) 0 0
6. Discussion and conclusion In this paper, we have introduced a family of generalized double-layer potentials and their corresponding 3
Fig. 3. Polygon approximation of two curves used for our experiments (left-most images). The polygon meshes consist in 201 and 29 segments respectively. Contour plots for: the signed distance, Ψ1 , Ψ3 and Ψ5 for a 2D polygon mesh (columns 2 to 5). Ψ5 (right-most images) provides already a good approximation of the distance function (second column). Note that the holes in the curve for the letter A are properly handled.
Fig. 6. Relative errors between dist(x) and Ψ2 (x) for a cube (left) and the rocker-arm model (right).
computation of Lp -distance fields can be significantly accelerated via a hierarchical approach: an accurate evaluation of a Lp -distance field at a given point does not require careful calculations of contributions of surface parts distant from the point. For polylines in 2D and triangle meshes in 3D, special correctors (boundary layers) should be added to the Lp -distance fields in order to achieve a better approximation accuracy at polyline/mesh vertices. We hope that our approach can lead to improving and enhancing existing applications of approximate distance functions in image processing, computer vision, computer graphics, geometric modeling, and several other branches of scientific computing [1,14,22]. In particular, we foresee potential applications of Lp -distance fields in computational mechanics where some new promising computational techniques rely on constructing accurate and efficient approximate smoothed distance functions and their derivatives [11, Section 3.2].
Fig. 4. Relative errors for a square. The first row corresponds to the relative error between the exact signed distance and ψ1 , ψ3 and ψ5 respectively. The second row uses the normalized functions: Ψ1 , Ψ3 and Ψ5 instead.
Appendix A: Proofs Fig. 5. Relative errors between the exact signed distance and Ψ5 for various opening angles between two segments. The angles vary between 0◦ (top row, left column) and 45◦ (bottom row, right column).
Proof of Proposition 1. An integral version of Jensen’s inequality (see, for example, [29, Problem 7.5]) states that Z Z F f (x)w(x) dx ≤ F (f (x)) w(x) dx,
Lp -distance fields. They possess interesting mathematical properties, and, in particular, the Lp -distance fields deliver smooth and accurate approximations of the distance function. This work is just a beginning of our study of generalized double-layer potentials and their applications. A practical
R
R
where F (t) is convex and w(x) is a positive weight function satisfying Z w(x) dx = 1. R
4
1 (k + 2)h2 h ∂ = k+2 − k+2 ∂h ρ ρ ρk+4 (see Fig. 7 for the notations used). Thus
Let 1 < p < q < ∞, F (t) = tq/p , R = [0, 2π], w ≡ 1/2π, p and f (θ) = (1/ρ(θ)) . Since F (t) is strictly convex, we have 1 1 Z 2π Z 2π dθ p dθ q 1 1 < , ψˆp (x) > ψˆq (x), 2π 0 ρ(θ)p 2π 0 ρ(θ)q
∂ϕk 1 = ϕk − (k + 2)hϕk+2 ∂h h or, equivalently,
1/k where ψˆk (x) = (2π) ψk (x). The multidimensional case n > 2 is considered similarly. Let x ∈ D is fixed and k → ∞. We have according to Laplace’s method (see, for example, [3]) Z Z 1 dΩ = exp k ln dΩ ϕk (x) = k ρ Ω Ω ρ
=
C1 + o(1) 1 , k (n−1)/2 ρkmin
(k + 2)ϕk+2 =
where ρmin = dist(x).
1/k k (n−1)/2 = dist(x) C1 + o(1) (n−1)/2 1/k = dist(x) [C2 + o(1)] k 1/k , as k → ∞, 1 ψk (x) = ϕk (x)
1/k
where C1 and C2 are positive constants. It remains to note that n−1 (n − 1) ln k 1 n − 1 ln k 2 1/k = exp = 1+ +O , k 2k 2 k k
with ck defined by (5), where singular potential ϕk (x) is generated by segment ab and x approaches an internal point of ab. For the potential generated by a polygon containing segment ab, we have a slightly bigger remainder: ϕk (x) = ck hk + O(1), as h → 0,
as k → ∞. Proof of Proposition 2. First let us consider the 2D case and study an asymptotic behavior of potential ϕ1 (x) generated by a single segment ab. See the left image of Fig. 7 for notations used. Assume that h, the distance from x to ab, is small: h 1. We have γ = π − α − β, a(h) = a + O(h2 ),
α = arctan (h/a) = h/a + O(h2 ),
b(h) = b + O(h2 ),
β = arctan (h/b) = h/b + O(h2 ),
γ 2ab 1 tan = + O(h), 2 a+bh
as the contribution of the other segments of the polygon is O(1). Thus 1/k ck Ψk (x) = = h + O hk+1 as h → 0. ϕk (x) In particular, if x is an internal point of ab, we have ∂ s Ψk (x)/ ∂ns = 0,
2 ϕ1 (x) = + O(h). h
In the 2D case, for a single straight segment, we have Z Z dθ h dsy ϕk (x) = = , ρ = |x − y|. ρ(θ)k ρ(θ)k+2 ∂ρ h = , ∂h ρ
∂ ∂h
1 ρk+2
=−
s = 2, 3, . . . , k.
Now let us consider a domain D with a smooth boundary S. Without loss of generality we can assume that D is convex (the case of non-convex domains can be treated similarly to how it is done in [7]). Consider a point x ∈ D on a distance h 1 from S. Let ab be a circular arc osculating S at the point closest to x such that the segment ab is perpendicular to the direction from x to its closest neighbor on S, as seen in the right image of Fig. 7. Let us also consider another circular arc cd such that it is situated inside D, ab and cd are collinear, and the direction from x to its closest neighbor on S is the bisector for cd. Now we approximate ϕk (x) given by (3) by the sum of integrals of 1/ρ(θ)k over the circular arcs ab and cd. Note that, as x approaches S, the osculating arc ab delivers a better and better local approximation of S and the leading term in asymptotics ϕk (x), as h → 0, is the same as that in the integral of 1/ρ(θ)k over the circular arc ab. The law of cosines yields q ρ(θ) = −(R − h) cos θ + R2 − (R − h)2 sin2 θ.
Fig. 7. Notations used for asymptotic analysis of potentials ϕk (x) generated by a segment (left) and a circular arc (right).
Note that p ρ = h2 + z 2 ,
(7)
Note that ϕk (x) is odd w.r.t. h. Thus the expansion of ϕk (x) w.r.t. h 1 contains only odd degrees of h. For example, 2 ϕ1 (x) = + a1 h + a3 h3 + a5 h5 + O h7 , h → 0. h Then (7) implies that 4 1 − 2a3 h − 4a5 h3 + O h5 ϕ3 (x) = 3 3 h and we can observe that (7) kills the linear w.r.t. h term in the expansion of ϕk (x). Now, applying the mathematical induction, we arrive at ϕk (x) = ck hk + O(h)
Thus
1 1 ∂ϕk ϕk − . h2 h ∂h
(k + 2)h , ρk+4 5
where R is the radius of the osculating circle. We can safely assume that R = 1. Our task now is to study the asymptotic behavior of
[9] P.-A. Fayolle, A. Pasko, B. Schmitt, and N. Mirenkov. Constructive heterogeneous object modeling using signed approximate real distance functions. Journal of Computing and Information Science in Engineering, ASME Transactions, 6(3):221–229, 2006. [10] M. S. Floater. Mean value coordinates. Computer Aided Geometric Design, 20(1):19–27, 2003. [11] M. Freytag, V. Shapiro, and I. Tsukanov. Finite element analysis in situ. Finite Elements in Analysis and Design, 47(9):957–972, 2011. [12] A. Ghosh, C. S. Rapp, and R. A. Friesner. Generalized Born model based on a surface integral formulation. J. Phys. Chem. B, 102:10983–10990, 1998. [13] L. Gorelick, M. Blank, E. Shechtman, M. Irani, and R. Basri. Actions as space-time shapes. IEEE Transactions on Pattern Analysis and Machine Intelligence, 29(12):2247–2253, 2007. [14] L. Gorelick, M. Galun, E. Sharon, R. Basri, and A. Brandt. Shape representation and classification using the Poisson equation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 28(12):1991–2005, 2006. [15] K. Hormann and M. S. Floater. Mean value coordinates for arbitrary planar polygons. ACM Transactions on Graphics, 25(4):1424–1441, 2006. [16] M. W. Jones, J. A. Baerentzen, and M. Sramek. 3D distance fields: a survey of techniques and applications. IEEE Transactions on Visualization and Computer Graphics, 12(4):581–599, 2006. [17] T. Ju, S. Schaefer, and J. Warren. Mean value coordinates for closed triangular meshes. ACM Transactions on Graphics, 24(3):561–566, 2005. ACM SIGGRAPH 2005 issue. [18] L. V. Kantorovich. Some remarks on Ritz’s method. Trydy vysshego voenno-morskogo inzhenerno-stroitel’nogo uchilishcha, (3), Leningrad, 1941. (Russian). [19] L. V. Kantorovich and V. I. Krylov. Approximate Methods of Higher Analysis. Interscience Publishers, 1958. Chapter 4, Paragraph 2. [20] S. J. Osher and R. P. Fedkiw. Level Set Methods and Dynamic Implicit Surfaces. Springer, 2002. [21] J. Peng. Thick Surfaces: Interactive Modeling of Topologically Complex Geometric Details. PhD, Department of Computer Science, New York University, 2004. [22] J. Peng, D. Kristjansson, and D. Zorin. Interactive modeling of topologically complex geometric detail. ACM Trans. Graph., 23:635–643, 2004. ACM SIGGRAPH 2004 issue. [23] T. Reiner, G. M¨ uckl, and C. Dachsbacher. Interactive modeling of implicit surfaces using a direct visualization approach with signed distance functions. Computers & Graphics, 35:596–603, 2011. SMI 2011 issue. [24] A. Ricci. A constructive geometry for computer graphics. The Computer Journal, 16(2), 1973. [25] A. N. Romanov, S. N. Jabin, Y. B. Martynov, A. V. Sulimov, F. V. Grigoriev, and V. B. Sulimov. Surface generalized Born method: A simple, fast, and precise implicit solvent model beyond the Coulomb approximation. J. Phys. Chem. A, 108(43):9323–9327, 2004. [26] V. L. Rvachev. Theory of R-functions and Some Applications. Naukova Dumka, 1982. Russian. [27] V. Shapiro. Semi-analytic geometry with R-functions. Acta Numerica, 16:239–303, 2007. [28] K. Siddiqi and S. Pizer. Medial Representations: Mathematics, Algorithms and Applications. Springer, 2008. [29] J. M. Steele. The Cauchy-Schwarz Master Class: An Introduction to the Art of Mathematical Inequalities. Cambridge University Press, 2004.
−k Zπ/2 q 2 2 (h − 1) cos θ + 1 − (1 − h) sin θ Ik (h) = dθ. −π/2
It is easy to verify that Ik (h) =
ck +O hk
1 hk−1
, as h → 0.
In the nD case, for a single (n − 1)-dimensional mesh tetrahedron (a mesh triangle in 3D), we have Z Z h dSy dΩy = , ρ = |x − y|. ϕk (x) = k ρ ρk+n So we arrive at a generalization of (7) for the multidimensional case 1 1 ∂ϕk (x) (k + n)ϕk+2 (x) = 2 ϕk (x) − . h h ∂h It is shown in [4] that vn−1 + O(1) ϕ1 (x) = h which justifies that c1 = vn−1 in (5). Proof of Proposition 3. Differentiating of (1) w.r.t x yields ∆φk (x) = k(k + n)φk+2 (x). For k = 1 and σ(y) ≡ 1 this formula coincides with that obtained in [4]. Obviously φk (x) is positive if surface ”dipole” density σ(y) is positive. Assume that φk (x) has a local maximum at some point inside the domain bounded by S. Then ∆φk (x) ≤ 0 at that point. We arrive at a contradiction which completes the proof. References [1] N. Ahuja and J.-H. Chuang. Shape representation using a generalized potential field model. IEEE Transactions on Pattern Analysis and Machine Intelligence, 19(2):169–176, 1997. [2] A. Biswas and V. Shapiro. Approximate distance fields with non-vanishing gradients. Graphical Models, 66(3):133–159, 2004. [3] N. Bleistein and R. A. Handelsman. Asymptotic Expansions of Integrals. Dover, 1986. [4] S. Bruvoll and M. S. Floater. Transfinite mean value interpolation in general dimension. J. Comp. Appl. Math., 233:1631–1639, 2010. [5] M. Carley. Potential integrals on triangles. Technical Report arXiv1201.4938, January 2012. [6] J.-H. Chuang, C.-H. Tsai, and M.-C. Ko. Skeletonization of three-dimensional object using generalized potential field. IEEE Transactions on Pattern Analysis and Machine Intelligence, 22(11):1241–1251, 2000. [7] C. Dyken and M. S. Floater. Transfinite mean value interpolation. Comp. Aided Geom. Design., 26:117–134, 2009. [8] R. Fabbri, L. da F. Costa, J. C. Torelli, and O. M. Bruno. 2D Euclidean distance transform algorithms: A comparative survey. ACM Computing Surveys, 40(1):1–44, 2008.
6
