The idea behind the morphological technique is to preserve the morphological features discussed in Chapter 4. Assuming the input is two solids that arer–
open and r–closed, the output should also have this property. The following procedure would yield the desired result:
• Reconstruct the original solids from their volumetric representation,
6.3 The Morphological Approach 113
• Perform the CSG operation on the reconstructed solids.
• Modify the result to ensure that the result fulfills the openness-closedness criterion.
• Voxelize once more to obtain a volumetric representation.
This scheme cannot be implemented directly, though, and, consequently, the technique operates quite differently, but produces the same result as the above.
Central to the approach is the link between morphology and propagating fronts.
For instance, dilation with a spherical structuring element can be implemented by pushing the boundary in the normal direction at unit velocity. This was exploited by Sapiro et al. who implemented Euclidean morphology in a discrete setting using Level–Set Methods [144, 143]. Another example is [164]. In some cases the front evolving in the normal direction might collide with itself causing self intersections, but if we view the propagating surface as front of burning material, self–intersections should be removed because the self–intersecting parts would propagate through material already burnt. This is known as the H¨uygens principle [156], and the Level–Set Method which will be discussed in Chapter 7 handles the problem according to that principle.
Observe that the boundary of a solid dilated with a sphere of radiusrconsists of points that are at a distance of exactlyr from the boundary of the original solid. This indicates that it is possible to implement the dilation of a solid (by a spherical structuring element) simply by finding the offset surface at a distance equal to the radius of the structuring element (See Figure 6.6). Of course, it is possible that the offset surface is self–intersecting and in that case it would not correspond exactly to the boundary of the dilated solid. Now, if the solid is r–closed, its r–level offset surface can at most touch itself as the surface is the locus of a ball rolling on the exterior side of the solid. At a point of self–
intersection the ball would be stuck (see Figure 6.5). This indicates that for r–closed solids, dilation can be implemented simply by finding ther–level offset surface. This is central to the implementation of the morphological approach to volumetric CSG. Until section 6.3.2 the approach will be discussed theoretically and only in terms of the signed distance functions.
We assume we are given two permissible,r–open andr–closed solidsS1 andS2
The union of these two solidsS1S
S2might not fulfill the criterion. However, it is easy to show that union preservesr–openness. This means that we only need to perform the close operation to ensure that the resulting solid is bothr–open andr–closed.
S=C(S1∪S2, br) (6.5)
r
r
r
Figure 6.5: A ball of radiusrrolling on the exterior side of the solid is stuck if ther–level offset surface is self–intersecting.
r S
S
d
r-d ⊕sr
p p’
Figure 6.6: Closest points on surfaces of dilated solids.
Of course, performing the close operation may ruin openness. This means that the input solids are assumed to be in a configuration so that when closed they remain open. Now, using the facts that close is dilation followed by erosion and that dilation distributes over union [151], we obtain
S= ((S1⊕br)∪(S2⊕br))⊖br (6.6) To simplify notation we will use these definitions in the following
ζ1=S1⊕br (6.7)
ζ2=S2⊕br (6.8)
ζ=ζ1∪ζ2 (6.9)
See Figure 6.7 for an illustration of theses definitions. The plan is now to (a) find the distance functionsdζ1 anddζ2 and then (b)dζ whence (c)dS is finally computed.
It is assumed that the solids S1 and S2 are r–open and r–closed. Hence, the r–level offset surface corresponds to the dilation with a ball of radius r, and
6.3 The Morphological Approach 115
I
I
ζ ζ
ζ 1
2
Figure 6.7: ζ, ζ1, ζ2, I
the shortest signed distance to the r–level offset surface is the distance to the surface minusr
dζ1(p) =dS1(p)−r (6.10)
dζ2(p) =dS2(p)−r (6.11)
Thus step (a) was trivial. It is, unfortunately, more complicated to perform step (b), but Proposition 6.1 is helpful. According to this proposition, if there is a pointq∈∂ζso thatkp−qk= min(dζ1(p), dζ2(p)) thendζ(p) = min(dζ1(p), dζ2(p)).
Now, where would the closest point be, if the condition does not hold? In that case, the next proposition (whose proof is in appendix B) gives us the answer:
Proposition 6.2 Given two permissible solidsS1 andS2 and a pointpso that
−2r < dζ(p)<0.
Bζi(p)∈/ ∂ζ =⇒ Bζ(p)∈I⊂∂ζ (6.12) where I=∂ζ1∩∂ζ2 (see Figure 6.7) andi= argminj∈{1,2}dζj.
Consequently, we can computedζ in the following way dζ(p) =
min(dζ1(p), dζ2(p)) Bζi(p)∈∂ζ
argminq∈I||p−q|| otherwise (6.13) where
i= argminj∈1,2dζj(p) (6.14)
This step was clearly a bit more complicated. In practice, we can check whether Bζ1(p)∈∂ζ simply by checking thatdζ2(Bζ1)>0 and vice versa. It is harder to find the closest point onI and that part has to be done numerically.
Fortunately, step (c) is also trivial. To perform the erosion, we simply addrto the distance value ofdζ. Thus
dS=dζ+r (6.15)
6.3.1 Examples
Assume that the input is two V–models V(S1) and V(S2) and the output is a new V–model V(S). Using the results above and the fact that the distance values are clamped we shall see how to compute the value ofV(S) at any given pointp.
The first step is to classifypaccording to the rules in table 6.1. Any point whose signed distance to the (un–dilated) solid is greater thanris called exterior. Any point whose distance is smaller than−ris called interior. As the table indicates,
state I T E
I I I I I=interior
T I TkI T T=transition
E I T E E=exterior
Table 6.1: Transition rules for volumetric union
points that are exterior to both solids remain exterior, and points that are interior with respect to either solid become interior. For instance,p2 andp3 in Figure 6.8 are exterior and interior, respectively. The values of the V–model at these points areV(S) =randV(S) =−r.
For points in the transition region of one solid which are simultaneously in the transition region or exterior to the other solid, more work is required. If the corresponding point on the surface of the dilated solid (sayζ1) is exterior to the other dilated solid, we simply use the value ofV(S1) in accordance with (6.13).
This case is exemplified byp0in Figure 6.8 whereV(S)(p0) =V(S1)(p0).
If the corresponding point on the surface of the dilated solid is an interior point in the other dilated solid, the problem is less trivial. We will call such points (or voxels)inconsistent in the following (in Figure 6.8p1is inconsistent).
For inconsistent points, we need to find the distance to I ⊂ ∂ζ according to
6.3 The Morphological Approach 117
S1 p0
S2
r r
p2
p3
I p1
∂ζ
Figure 6.8: Point classification
Proposition 6.2. If the distance is greater than 2r, we ignore the value and the point is classified as being interior.
In summary, the algorithm should work like previous volume CSG algorithms based on (6.2), except that for some points we need to estimate the distance to I. Hence, the main difficulties lie in representingIand finding the inconsistent points.
6.3.2 Implementation
We shall look at how the algorithm actually operates and take into account that the operands are not signed distance functions or V–models (i.e. clamped, signed distance functions). The algorithm works on two operands: The left operand is the voxel grid being modified, G = V(S1), and right operand is the tool,V(S2), which is an un–sampled V–model since sampling it before the CSG operation would only complicate matters. Except, of course, if we want to do cut and paste operations. CSG operations between two voxel grids is a feature which has also been implemented, but in this case the right operand is interpolated and thus takes on the rˆole of a V–model.
The result of the CSG operation is assigned back to the grid:
G←G∪vV(S2) (6.16)
The algorithm works in two passes, and both passes traverse all voxels.
First pass The goal of the first pass is to find and tag inconsistent points and to find a set of points belonging to I. For each voxel two operations are performed:
1. A foot point is estimated on either∂ζ1 or∂ζ2 depending on which distance is smaller. If the foot point is interior with respect to the other solid, the voxel is tagged as being inconsistent. My volume representation contains gradient information, which speeds up this process, although the gradient could also have been estimated.
2. If the estimated value of Bζ1(p) is within 12 vu distance to∂ζ2, the closest point inI is estimated by assuming that∂ζ1and∂ζ2are planes and finding the point on the intersection of these planes that is closest toBζ1(p). The point is added to a set of points that make up the estimate ofI.
Second pass In the second pass, the new voxel values are computed. This pass amounts to an implementation of (6.13).
For each voxel, a case analysis is performed according to table 6.1. For voxels that are in the transition region of one solid and either exterior or in the tran-sition region of the other, it is checked whether they are tagged as inconsistent.
For un-tagged voxels the new value is simply min(G,V(S2)). For inconsistent voxels, the closest point in the set of points representingI ⊂∂ζ is found, and the distance to that point is used to calculate the new voxel value.
Estimating the closest point in I is one of the trickiest parts of the algorithm, and it cannot be done simply by finding the closest point in the set ofI–estimate points generated during the first pass, because these estimates may be as far apart as 1 vu if the surfaces of bothS1 and S2 are parallel to coordinate axes in the grid. Hence, the distance to anI–estimate would sometimes deviate too much from the true distance toI. To solve this problem, we store an estimated tangent direction vector with (nearly) all points in the point set representing I. To estimate the distance from a voxel to I, we find the closestI–estimate and the projection of the voxel onto the line defined by theI–estimate and the associated direction vector. The projected point is used as the estimate of the closest point inI, and this point is used to update both the distance and normal direction of the voxel.
In some cases, the union of two solids may contain a corner. In these cases,I bends sharply and the direction vector associated with the I–estimate at the bend would be misleading. The solution is to find the closestI–estimate before and after any given I–estimate. If the angle defined by a given estimate and