fayolle@u-aizu.ac.jp, apasko@bournemouth.ac.uk

Contents

• 1. Introduction
• 2. Distance to the intersection of two orthogonal halfspace
• 3. Distance to the intersection of two arbitrary solids
• 4. Modifications for the union and the difference
• 5. Results and discussion
• 6. Conclusion
• Reference

Abstract

          We present in this paper methods to compute the signed Euclidean distance to surfaces obtained by the intersection (respectively union or difference) of two solids (in two or three dimensions). These implementations can replace min/max or R-functions traditionally used to model set operations used with implicit surfaces.

          Keywords:  shape and solid modeling, constructive modeling.

1. Introduction

          In shape modeling, a solid can be defined by the sign of a continuous function: the set {p : f(p) ≥ 0} defines the interior and the boundary and the set {p : f(p) < 0} the exterior (the so called implicit surfaces or F-Rep, see [3], [11] and the references therein). The particular case where f is the signed Euclidean distance to the boundary of the solid is of special interest. Distance based models are extremely useful in many applications such as: constant-radius offsetting and blending operations [16], surface metamorphosis and smoothing [13], object reconstruction from a set of cross-sections [9], rendering with sphere tracing [5], generation of skeletal shape representation [24], heterogeneous object modeling [2] and others. We present in this paper methods for calculating the signed Euclidean distance from a point to the surface of a solid constructed by applying set-theoretic operations (union, intersection, difference) to primitives defined by distance functions. That is: if d₁ and d₂ represent the distances to the surfaces of two primitives S₁ and S₂, we present an algorithm for calculating the distance d to the surface of the union (respectively intersection, difference) of S₁ and S₂.

1.1 Related works

1.1.1 The distance function

Let d(p), p R² or R³ be the signed distance function to an oriented closed surface M. The function d is the vanishing viscosity solution of the Eikonal equation [23, 21, 19]:

(1)

where ∥.∥₂ – the Euclidean norm.

Let c be the closest point to p in the surface M, the signed distance is then ϵ∥p - c∥₂, where ϵ = -1 if p is outside the solid delimited by the surface M, 1 otherwise. If the surface is smooth, then p - c is orthogonal to the surface. The signed Euclidean distance function is continuous but not everywhere differentiable.

Expressions for the distance function to most of the typical surfaces of a CSG system (sphere, cylinder, cone) are known analytically [5] and distance to general quadrics and ellipsoids can be computed by a numerical procedure [6].

In general, if the surface M is available as an oriented point-set or a mesh of polygons, it is possible to solve the Eikonal equation (Eq. (1)) on a finite grid. There exist various numerical algorithms such as the fast marching method [19], the fast sweeping method [21, 23], or the characteristics / scan conversion algorithm [10]. Algorithms exploiting the GPU have also been designed in order to compute efficiently the Euclidean distance function on a grid [7, 20]. Once a grid is obtained, with the signed Euclidean distance to M sampled at each grid node, it is possible to apply spline interpolation or approximation to get a function approximating the distance [15].

1.1.2 Constructive geometry with real valued functions

In constructive geometry, complex solids are built by applying successively set operations (union, intersection, difference) to primitives. For solids described using real valued functions, like in implicit surfaces or F-Rep modeling systems, expressions for these operations have been proposed by Sabin [18], Ricci [14] and Rvachev [17].

Sabin [18] and Ricci [14], independently proposed the use of the functions min and max to implement the intersection and union of implicit surfaces. If d₁ and d₂ are the functions defining two solids, then the union is defined by the function max(d₁,d₂) and the intersection by min(d₁,d₂) (the difference operation is obtained by min(d₁,-d₂), i.e. by replacing d₂ by -d₂). Rvachev proposed the R-functions [17]:

	d1 ∨αd2 = (d1 + d2 + )	(2)
	d1 ∧αd2 = (d1 + d2 -)	(3)
	d1 \αd2 = d1 ∧α(-d2)	(4)

Neither the scalar field obtained by using min/max nor the R-functions correspond to the Euclidean distance to the surface of the constructed solid. This is illustrated in Fig. 1, where the approximate distance field to the surface made by the intersection of the halfspaces: x ≥ 0 and y ≥ 0 is computed using an R-function (left) and min (right). The red curve corresponds to the solid boundary (i.e. corresponds to the distance 0) and the blue curve to the set of points at a value of -1 outside of the intersection.

Figure 1.     Approximate distance to the intersection of the halfspaces x ≥ 0 and y ≥ 0. The red curve corresponds to the boundary of the halfspaces intersection and the blue curve corresponds to the set of points at a value -1 outside of the intersection. Left: R-function is used to implement the intersection. Right: min is used.

Scalar fields constructed using the min/max functions keeps a better approximation of the distance field than when the R-functions are used. This can be illustrated by computing the union of a disk with itself and inspecting the value of the function (representing the distance to the union of the disks) at the center of the disk. The distance to a circle of radius 1 and center [0,0] is: d(p) = 1.0 - and the union of the disk with itself is defined by: max(d(p),d(p)) (when using min/max) or: d(p) ∨₀d(p) (when using the R-functions). In the former case, the distance at the center is: 1.0 while in the latter it is: 3.41421.

The functions min/max are not differentiable on the set of points corresponding to the equality of their arguments (i.e. for all points p where d1(p) = d2(p)). Because of this, R-functions are sometimes preferred in modeling systems; R-functions are not differentiable only on the set of points where their arguments are both equal to 0. Smoothness of the R-functions is used in defining some solid modeling operations such as for example when implementing blending between shapes (see for example the description proposed in [12]).

Some works tried to address this issue by modifying the contour lines of the functions min(x,y) and max(x,y): functions proposed in the works [8, 1] were designed for defining blending operations, while functions proposed in [4] were designed for keeping the distance approximation of min/max while removing the points where the function is not C¹.

1.2 Overview and main contributions

The main contributions of this paper are methods to compute (in two and three dimensions) the Euclidean distance to the surface of solids defined by the union (respectively intersection, difference) of solids defined by distance functions. That is: if S₁ and S₂ are solids defined by the distance functions d₁ and d₂, we describe methods for computing d the signed distance to the boundary of S = S₁ ∪ S₂ (respectively S₁ ∩ S₂ and S₁\S₂). These expressions for the set operations can be used instead of min/max or the R-functions in implicit surfaces or F-Rep modeling systems.

Computation for the union and intersection operations are symmetric and the difference operation can be obtained from the expression of the intersection, so we limit the discussion to the construction of an expression for the intersection. First, we investigate in section 2 the case of the intersection of two orthogonal halfspaces. Then in section 3, we describe how to compute the distance to the intersection of two general solids. In section 4, modifications for obtaining the distance for the union and difference operations are given. Finally, two and three dimensional examples are used to illustrate the behavior of these functions and to compare them with min/max and R-functions.

 2. Distance to the intersection of two orthogonal halfspaces

Before describing the computation of the distance to the intersection of two solids defined by arbitrary functions, we start by a simple example: the distance to the intersection of two orthogonal halfspaces given by: x ≥ 0 and y ≥ 0 (as illustrated in Fig. 1 and 2). In particular, we explain how the scalar field constructed by min(x,y) differs from the exact distance to the intersection.

2.1 Intersection of two orthogonal halfspaces

Let us consider the example used in Fig. 1: two halfspaces are defined by x ≥ 0 and y ≥ 0. Their intersection is the domain given by {(x,y) : x ≥ 0 & y ≥ 0}.

Figure 2 (left) illustrates some contour lines of the distance field to the boundary of the domain previously defined. The distance field is signed such that points inside the intersection of the halfspaces are positive and points outside are negative. Fig. 2 (right) shows some contour lines of the scalar field defined by min(x,y).

Figure 2.     Left: contour lines of the Euclidean distance to the boundary of the intersection of two halfspaces x ≥ 0 and y ≥ 0. Right: contour lines of the scalar field min(x,y). For both cases, contour lines vary from -1.5 to 1.5 with step of 0.5.

As illustrated by Fig. 2, the scalar field obtained by min corresponds to the Euclidean distance field everywhere except in the domain defined by: Q := {(x,y) : x < 0 & y < 0. For each point in that domain, the closest point on the boundary is the point obtained by the intersection of the two lines x = 0 and y = 0. Contour lines in that domain are circular arcs centered at this point.

The expression of Q for two orthogonal halfspaces is simple. We need now to define it for any general shapes to be intersected and then explain how to compute the closest point on the surface in that domain.

3. Distance to the intersection of two arbitrary solids

3.1 Determination of the domain where min differs from the Euclidean distance

Figure 3 illustrates the case of the local intersection of two general (non orthogonal) curves intersecting at a point O. The boundary (d = 0), one contour line (d = -1.0) and the normals to each curve at O are drawn.

If the closest point is known, then the region of interest can be identified using the normal to each curve at the intersection point as seen in Fig. 3. However, since it requires to know the closest point, this method is not feasible. Instead, we proceed in the reverse way by first determining the domain of interest and then computing the closest point in that domain.

Given a point p R² or R³ (see Fig. 4 for an illustration in two dimensions), we first compute the projection p₁ of p on the first curve c₁. Using the signed distance field d₁ to the curve and its gradient ∇d₁(p), the projection is given by:

(5)

Similarly, we compute p₂ the projection of p on the second curve c₂. p belongs to the zone of interest Q, if d₁(p₂) < 0 and d₂(p₁) < 0.

Figure 4.     Two curves c1 and c2 (solid lines) defined by the iso-lines of the distance fields d1 = 0 and d2 = 0. A given point p and the contour lines of d1 and d2 passing through p. p1 and p2 the projections of p on c1 and c2.

After determining if a point is in the domain where min applied to the distance fields differs from the Euclidean distance field, we need to determine the closest point on the boundary of the solid obtained by the intersection of the two solids.

3.2 Determination of the closest point

In all the above examples, the closest point on the boundary of the intersection from a point in the domain Q happened to be at the intersection of the two curves. This is not always the case. The closest point (x,y,z) on the boundary of the intersection of S₁ and S₂ from a point (x₀,y₀,z₀) in Q is the solution to the following constrained optimization problem:

Minimize:		(6)
such that:	d1(x,y,z) ≥ 0	(7)
and:	d2(x,y,z) ≥ 0	(8)

Where: d₁ and d₂ are the distance fields to the boundary of the solids S₁ and S₂. In Mathematica^tm [22], such a constrained minimization problem can be solved with the command NMinimize. The point (x₀,y₀,z₀) can be used as a starting point. Since this point is in general close to the solution, the convergence is relatively fast.

3.3 Distance to the intersection

Regrouping all the results obtained so far, a procedure for computing the distance to the boundary of the intersection of two solids defined by signed distance functions is:

4. Modifications for the union and the difference

Modification of the previous algorithm for obtaining the distance to the difference of two solids consists in replacing the function d₂(.) by the function -d₂(.) everywhere.

For the union, the following changes need to be made to the previous algorithm for the intersection:

• Line 7 needs to be changed to: inQ = d₁(p₂) > 0 AND d₂(p₁) > 0,
• Line 10 needs to be replaced by: vd = ∥p - p_c∥₂

Finally, the closest point computed in line 9 is obtained by solving the following constrained minimization problem:

Minimize:		(9)
such that:	d1(x,y,z) ≤ 0	(10)
and:	d2(x,y,z) ≤ 0	(11)

5. Results and discussions

5.1 Visualization of the resulting distance fields

Scalar fields corresponding to the intersection and union of two solids is illustrated in two and three dimensions in the following examples. Min/max, R-functions and the methods described in this paper are used to compute these scalar fields.

5.1.1 Contour plots in two dimensions

Figure 5 illustrates the contour plot of signed distance functions obtained by the intersection of two non-orthogonal planes and by the intersection of two disks. Intersection was implemented using: R-functions (left), min/max (middle) and the method proposed in this paper (right). The functions built using R-functions clearly differ from the Euclidean distance to the boundary of the constructed solid. The functions built using min/max fail also to do so but only in a subset of R². In comparison, the method described in this paper allows the construction of an accurate signed distance field to the boundary of the constructed solids as illustrated by the circular arcs in the contour lines outside of the constructed solids.

Figure 5.     First row: contour plots of the scalar fields for the intersection of two halfplanes built using(from left to right): R-functions, min and the method described in this work. Second row: Result for the intersection of two disks.

Figure 6 illustrates contour plots of the distance fields to the boundary of the union of two halfplanes (left) and two disks (right) using the method introduced here.

Figure 6.     Left: contour plot of the distance to the union of two halfplanes with our method. Right: contour plot of the distance to the union of two disks.

5.1.2 Contour plots in three dimensions

Visualizing scalar fields in three dimensions is more difficult. We visualize instead contour plots of the fields in a planar section. Fig. 7 (right) illustrates the proposed method applied in three dimensions to compute the distance to the intersection of two spheres. Contour plots need to be compared with the contour plots shown in fig. 7 (left) and (middle) obtained with respectively R-functions and min.

Figure 7.     Distance to the intersection of two spheres using for the intersection operation (from left to right): R-functions, min and the method proposed in this work.

Figure 8 illustrates the result of the proposed method to compute the distance to the union of two spheres by showing a contour plot of the constructed distance field on a section by the plane y = 0.

5.2 Discussion

5.2.1 Differentiability

The distance function is continuous but not everywhere differentiable; for example: the function d(p) = 1 -∥p - p_c∥₂ is not differentiable at the point p = p_c. The algorithms described above use the gradient of the functions passed as arguments. Points at which the gradients can not be computed need to be handled differently. These points are on the medial axis of the solid and have more than one closest point on the surface of the solid. A possible solution is to select one of the closest point on the surface and define the gradient as the unit vector ending at this point.

5.2.2 Successive applications of set operations

It is possible to successively apply the proposed methods for intersection, union or difference to solids built using these operations. The gradient of the function constructed by the methods described in this paper needs to be computed. In the example of the intersection, the function is not differentiable at the points: {p : pQ and d₁(p) = d₂(p)}. This point-set corresponds to the medial axis of the resulting solid. It is possible to use the same method as described above. An alternative solution is to use modified versions of min and max. In the case of min, it is possible to write: min(d₁,d₂) = (d₁ + d₂ + ) and to add a small perturbation ϵ: (d₁,d₂) = (d₁ + d₂ + ). This function is C¹ however it has the inconvenient of displacing each iso-contour; for example: assuming that d₁ = 0 and d₂ > 0, then (d₁,d₂) = (d₂ + ) < 0 instead of 0.

5.3 Complexity

The proposed expressions for the set operations are in general slower than min/max or the corresponding R-functions. This is due to the necessity to solve a constrained minimization problem in some cases. We have only done experiments with the NMinimize function of Mathematica^tm and some simple algorithms to solve this minimization problem. Since the objective function is not too complicated, there should be some more efficient methods. In addition, it is possible to provide a good approximate guess for the starting point, making the convergence relatively fast.

6. Conclusion

We have presented in this paper algorithms for computing in two and three dimensions the distance to a solid obtained by applying set operations (intersection, union or difference) to primitives. Contrary to the existing implementations of these operations (min/max and R-functions), the proposed methods allow to compute the signed Euclidean distance field to the boundary of the resulting shape.

We are considering extending this work by experimenting with the iterative application of these operations in order to construct complex solids. The use of these functions should also allow to easily implement rolling ball blend and double offsetting.

References