Log In Sign Up

Half-Space Power Diagrams and Discrete Surface Offsets

by   Zhen Chen, et al.

We present a novel algorithm to compute offset surfaces of shapes discretized using a dexel data structure that is both fast and trivially parallel. We achieve this by exploiting properties of half-space power diagrams, where each seed is only visible by a half space. This allows us to develop a two-stage sweeping algorithm that is both simple to implement and efficient, and avoids computing a 3D volumetric distance field. The primary application of our method is interactive modeling for digital fabrication, where printed objects have a finite resolution. However, our method can also be used for other purposes where a fast but approximate offset solution is desirable. We present experimental timings, how they compare with previous approaches, and provide a reference implementation in the supplemental material.


page 1

page 3

page 9

page 10


Tropical diagrams of probability spaces

After endowing the space of diagrams of probability spaces with an entro...

Border Algorithms for Computing Hasse Diagrams of Arbitrary Lattices

The Border algorithm and the iPred algorithm find the Hasse diagrams of ...

SurfaceVoronoi: Efficiently Computing Voronoi Diagrams over Mesh Surfaces with Arbitrary Distance Solvers

In this paper, we propose to compute Voronoi diagrams over mesh surfaces...

A simple and fast algorithm for computing discrete Voronoi, Johnson-Mehl or Laguerre diagrams of points

This article presents an algorithm to compute digital images of Voronoi,...

An Efficient, Practical Algorithm and Implementation for Computing Multiplicatively Weighted Voronoi Diagrams

We present a simple wavefront-like approach for computing multiplicative...

1 Introduction

Morphological operations such as dilation and erosion have numerous applications in various areas of computer science. They can be used to regularize shapes [Williams:2005:MMS], ensure robust designs in topology optimization [Sigmund:2007:MBB], perform collision detection [Teschner:2005:CDF], or compute image skeletons [Maragos:1986:MSR]. Surface offsetting operations are also extremely useful in the context of digital fabrication [Livesu:2017:FDM]: erosion can be used to generate support structures [Hornus:2016:TPE], to hollow an object or create a mold, closing operations can remove small holes in a model, etc.

While offset surfaces can be computed exactly by way of Minkowski sums [Hachenberger:2007:EMS], these operations can be fairly slow, especially on large models. Recent approaches [Meng:2017:ECF] can provide better results, but performances are still a concern, limiting their use in interactive applications. Conversely, approximate algorithms that rely on a discrete re-sampling of the input volume can lead to more efficient solutions, while sacrificing accuracy [Wang:2013:GBO, Martinez:2015:CSO] in a controlled way. This is especially interesting for applications in digital fabrication, where the resolution is anyhow limited by the machines and it is thus unnecessary to achieve infinite accuracy.

We present a novel algorithm to compute offset surfaces on a ray-based representation (ray-rep) of a solid. A ray-rep, such as the dexel buffer, sores the intersection results of a solid object with a set of parallel rays in a uniform 2D grid – each cell of the grid holds a list of intervals bounding the solid, as illustrated Figure 8. Ray-reps are extremely appealing in the context of digital fabrication, as they allow us to perform a number of operations directly in image-space, at the resolution of the machine. The computed result is thus identical to what is obtained by solving the problem exactly and then rounding it to machine resolution, but at a fraction of the cost. This includes CSG operations, but also support other common tasks such as support detection, toolpath palling, and infill calculations. Implicit surfaces can also be discretized directly without prior explicit meshing. Examples in the literature can be found for Computer Numerical Control (CNC) milling applications [VanHook:1986:RTS], modeling for additive manufacturing [Lefebvre:2013:IAG], hollowing, or contouring [Livesu:2017:FDM]. The offset algorithm we present in this article is both fast and accurate, in that it computes the exact offset of its input discrete dexel structure, in contrast to existing methods such as [Martinez:2015:CSO]. Our algorithm is also embarrassingly parallel, and its complexity scales well with the offset radius.

(a) Input Surface
(b) Thick Shell
(c) 3D Printed
Fig. 4: 3D printed surface shell, computed with our approach. The shell is the Boolean difference of a dilation and an erosion of the input volume. The model uses a grid of dexels, and a dilation radius of 4 dexels.

Our algorithm relies on a novel way to compute half-plane Voronoi diagrams [Fan:2011:HPV]. We first discuss the 2D case and present a simple and efficient sweeping algorithm to compute offsets in a 2D dexel structure. We then extend our sweeping procedure from 2D to 3D by leveraging the separability of the Euclidean distance, similarly to existing signed-distance function algorithms [Meijster:2002:AGA].

In Section 4, we compare running times of our algorithm with existing approaches, and discuss existing methods. Finally, we show different applications of our approach, including topological simplification and modeling for additive manufacturing.

2 Related Work

Exact Mesh Offsets

If the input shape is represented as a triangle mesh, a way to compute a dilated surface is to resolve the self-intersection in the raw offset surface [Jung:2004:SIR]. In general, a dilated mesh can also be obtained by computing the Minkowski sum of the input triangle mesh and a sphere of the desired radius [Hachenberger:2007:EMS]. This approach, which is implemented robustly in CGAL, can become slow for large models with complex geometry. Other algorithms avoiding the use of arbitrary precision arithmetic have been proposed to extract the boundary polygons of a Minkowski sum [Campen:2010:PBE]. Still, those methods are sensitive to the quality and complexity of the input mesh, and do not provide interactive runtimes.

Resampled Offsets

Rather than extracting the boundary of the Minkoswki sum directly, which is a challenging and time-consuming operation, other methods chose to resample the offset surface, by computing the isosurface of the signed-distance function to the original surface [Varadhan:2006:AMS, Peternell:2007:MSB]. A more recent approach by Meng:2017:ECF distribute and optimize the position of sites at a specified distance from the original surface, and produces a triangle from a restricted Delaunay triangulation of the sites. Calderon:2014:PMO introduced a framework for performing morphology operations directly on point sets.

Voronoi Diagrams and Signed-Distance Transforms

Our algorithm rely on the implicit computation of Voronoi and power diagrams of points and segments. Centroidal tesselations of Voronoi and power diagrams have important application in geometry processing and remeshing [Liu:2009:OCV, Xin:2016:CPD].

In a method proposed by [Fortune:1987:ASA] to a 2D Voronoi diagram of a point set, a sweep line algorithm is employed to position of the Voronoi vertices (intersection of 3 bisectors) that form the Voronoi diagram. By relying on a slight modification of the Voronoi diagram definition, called half-space Voronoi diagrams [Fan:2011:HPV], we will see in Section 3 how to compute the Voronoi diagram of a set of parallel segments directly, using two sweeps instead of one, and how it lead to the efficient computation of a discrete offset surface. The extension of our method to 3D will require the computation of power diagrams [Aurenhammer:1987:PDP], where each seed is associated a different weight. Note that our algorithm implicitly relies on the Voronoi and power diagrams of line segments. And while those can be efficiently approximated by sampling each segments with multiple points [Lu:2012:CVT], we will see how to circumvent the approximation entirely and keep the implementation simple.

Finally, our 3D twos-stage sweeping algorithm bears some similarity to existing signed-distance transform approaches such as [Meijster:2002:AGA, Maurer:2003:ALT], as it leverages the separability of the Euclidean distance to compute the resulting offset by sweeping in two orthogonal directions. However, as we operate on a dexel structure, we never need to store a full 3D distance field in memory. For a more complete review of distance transform algorithms, the reader is referred to the survey [Fabbri:2008:DED].

Ray-Based Representations and Offsets

A ray-based representation of a shape is obtained by computing the intersection of a set of rays with the given shape. In most applications, the cast rays are parallel and sampled on a uniform 2D grid, then stored in a structured called a dexel buffer. A dexel buffer is simply a 2D array, where each cell contains a list of intersection events containing the segments intersecting with the solid. See Figure 8 for an illustration. To the best of our knowledge, the term dexel (for depth pixel) can be traced back to [VanHook:1986:RTS], which introduced the dexel buffer to compute the results of CSG operations to facilitate NC milling path-planning. Similar data structures have been described in different context over the year. Layered Depth Images (LDI) [Shade:1998:LDI] are used to achieve efficient image-based rendering. The A-buffer [Carpenter:1984:TAB, Maule:2011:ASO] was used to achieve order-independent transparency. While they may differ in the way they are built, the underlying data structure remains similar. Augmenting ray-reps with normal information, the G-buffer [Saito:1991:NMW] stores a normal in every pixel for further image-processing and NC milling applications. In the context of additive manufacturing, Layered Depth Normal Images – LDI with normal information – have been proposed as an alternative way to discretize 3D models [Wang:2010:SMO, Huang:2014:AFL].

Ray-based data structures offer an intermediate representation between usual boundary representations (such as triangle meshes), and volumetric representations (such as a full or sparse 3D voxel grids). While both 3D images and dexel buffers suffer from uniform discretization errors across the volume, a dexel buffer is cheap to store compared to a full volumetric representation. This allows to discretize with very large dexel-buffer resulutions (1024 to 2048 and more) at a very low cost. Additionally, in the context of digital fabrication, 3D printers and CNC milling have limited precision. Consequently, using a dexel buffer at the resolution of the printer is enough to cover the space of shapes that can be fabricated.

Ray-reps have been used to compute and store the results of Minkoswki sums, offsets and CSG operations [Hartquist:1999:ACS]. Hui:1994:SSI use ray-reps to compute the result of a solid sweeping in image-space. In [Chen:2011:UOO], Chen:2011:UOO use LDNI to offset polygonal meshes, by filtering the result of an initial overestimate of the true dilated shape. Wang:2013:GBO compute the offset mesh as the union of spheres placed on the points sampled by the LDI. As the dilation radius increases, the number of overlapping spheres quickly increases, and they propose to decompose the dilation in a sequence of dilation with smaller radii. Finally, [Martinez:2015:CSO] approximate the dilation by a spherical kernel with the dilation by a zonotope, effectively computing the Minkoswki sum of the original shape with a sequence of segments in different directions. The method presented in this paper computes the exact offset of the input ray-rep, as opposed to [Martinez:2015:CSO]. It also differs from Wang:2013:GBO in several aspects. The latter requires a LDI sampled from 3 orthogonal directions, while our method accelerates the offset operation even when a ray-rep from a single direction is available. We describe other key differences in more details in Section 4.

3 Method

We first introduce the different concepts and notations in Section 3.1, where we present a general overview of our algorithm for offset computation of ray-reps. In Section 3.2, we describe a 2D version of our method, based on efficiently computing the Voronoi diagram of a set of parallel segments. Then, in Section 3.3, we extend the algorithm to compute the power diagram of a set of parallel segments in 2D, and observe that it is the key to extend our 2D offset algorithm to 3D.

3.1 Definitions and Overview

Dexel Representation

A dexel buffer of a shape is a set of parallel segments arranged in a 2D grid, where each cell of the grid contains a list of segments sharing the same coordinate. More specifically, we discretize as where . Each segment in the same cell have the same coordinate, and represent the intersections of the input shape with a vertical ray at the same coordinate (see Figure 8).

(a) Input Surface.
(b) Dexel Approximation.
(c) Compact Storage.
Fig. 8: A dexel data structure constructed from a triangle mesh. The input surfage LABEL:sub@fig:dexels:input is intersection with evenly spaced parallel rays LABEL:sub@fig:dexels:output, and stored compactly as a 2D grid of events in the dexel buffer LABEL:sub@fig:dexels:compact. This can be used to perform efficient CSG operations in modeling software [Lefebvre:2013:IAG].


Given an input shape and a radius we define the dilated shape as the set of points at distance less than from :



The erosion of a shape is equivalent to computing the dilation on the complemented shape , and taking the complement of the result. Since the complement operation is trivial under a ray-rep representation, we will focus on describing our algorithm applied to the dilation operation, from which all other morphological operators (erosion, closing, and opening) will follow.

Voronoi Diagram

The Voronoi diagram of a set of seeds is a partition of the space into different Voronoi cells :


We also define to be the interface between each overlapping Voronoi cell:


When the seeds are points, is a set of straight line in 2D (planes in 3D), which are equidistant to their closest seeds. When the seeds are segments, is comprised of parabolic arcs in 2D [Lu:2012:CVT], and parabolic surfaces in 3D.

In a half-space Voronoi diagram [Fan:2011:HPV], each seed point is associated with a visibility direction , and a point is considered in Equation 2 if and only if . In this work, we are interested in half-space Voronoi diagrams of seed segments, where each seed is associated to the same visibility direction . Figure 9 shows the difference between a Voronoi diagram and half-space Voronoi diagrams for point seeds. More precisely, given set of parallel seed segments , and a direction orthogonal to each segment , we define the half-space Voronoi cells and as


where is the point projected on the line . Note that the segments are parallel and orthogonal to the chosen direction , so the dot product in Equation 4 has the same sign for all the points in the line . Note that we have (see Figure 9).

[width=0.33]figs/img/voronoi_fwd.png [width=0.33]figs/img/voronoi_rev.png

Fig. 9: Voronoi diagram of seed points . From left to right: Voronoi diagram formed by the full Voronoi cells ; half-space Voronoi diagram formed by the half-space Voronoi cells and respectively. Note that . In each diagram, one Voronoi vertex (intersection between Voronoi cells) is shown with a red square.

Half-Dilated Shape

We define the half-dilated shape , as the dilation restricted to the half-space Voronoi cells of segments in :


Remember that is the -th segment in the dexel buffer approximating the shape . The half-space dilated shape is defined in a similar manner using . For simplicity, we will omit the dilation radius and simply write , unless there is an ambiguity.

Power Diagram

Finally, the power diagram of a set of seeds is a weighted variant of the Voronoi diagram. Each seed is given a weight that determines the size of the power cell associated to it:


Intuitively, one could interpret a 2D power diagram as the orthographic projection of the intersection of a set of parabola centered on each seed, where the weights determine the height of each seed embedded in . This definition extends naturally to half-space power diagrams.


The input of our algorithm is a dexel buffer and a dilation radius . The output is a new dexel buffer . A key observation in our approach is that the dilation operation can be separated into two stages. The first stage computes the dilation of the shape along the Y axis, producing in intermediate shape . Then, the second stage will compute the dilation of along the X axis by a radius , producing the final shape . This separation allows us to decompose the two dilation operations into a set of separate 2D dilations along each axis. The first stage, explained in Section 3.2, covers the basic 2D dilation. The second stage (Section 3.3) involves augmenting each ray segment with a weight that determines the radius of the dilation that should be performed along the second axis in the second stage. This decomposition is illustrated in Figure 14.

Fig. 14: Dilation operation performed in two stages LABEL:sub@fig:twopass:1. A seed is dilated by a radius along a first axis (in red). The resulting seed (in green) is then dilated along a second axis by a different radius , where is the distance between the two seeds (red point and green point). The equivalent dexel data structure for each pass is shown on right: LABEL:sub@fig:twopass:2 input dexel, LABEL:sub@fig:twopass:3 result of the first pass, LABEL:sub@fig:twopass:4 result of the second pass.

3.2 Half-Space Voronoi Diagram of Segments and 2D Offsets

Given a 2D input shape , and dilation radius , we seek to compute the dilated shape , comprised of the set of points at distance from . The input shape is given as a union of disjoint parallel segments evenly spaced on a regular grid (the dexel data structure), and we seek to compute the output shape as another dexel structure (the discretized version of the continuous dilated shape).

The dilated shape can be expressed as the union of the dilation of each individual segment of . To compute this union efficiently, our key insight is to partition the dilated shape based on the Voronoi diagram of the input segments. Within each Voronoi cell , if a point is at a distance from the seed segment , then .

Fig. 15: Fortune’s sweepline algorithm [Fortune:1987:ASA] requires transforming the point coordinates to compute the correct Voronoi diagram of points in a single sweeping pass, which makes it impossible to compute the result of the dilation in the same sweeping pass.

The Voronoi diagram of point seeds can be computed in single pass with a sweepline algorithm [Fortune:1987:ASA], the construction requires lifting coordinates in the plane according to the coordinate along the sweep direction, as illustrated in Figure 15. Unfortunately, this approach cannot be used to compute the result of a dilation operation in a single pass. Indeed, seeds need to influence the rows before and after their current position along the sweep direction. Instead, we propose a simple construction for seed points and segments, based on half-space Voronoi diagrams – which we extend to power diagrams in Section 3.3. Our key idea is to compute the dilation of a segment in its Voronoi cell as the union of the two half-dilated segments, in and . Both and can be computed efficiently in two separate sweeps of opposite direction, without requiring any transformation to the coordinate system as in [Fortune:1987:ASA].

Fig. 16: A single sweep in one direction allows us to compute the half-space Voronoi diagram of parallel line segments (the dexel data structure).


Fig. 17: The bisector of two line segments can be described by a piecewise second-order polynomial curve. Voronoi vertices are located at the intersection between two such bisectors. After this point, the middle segment will always be further away from the sweepline than its two neighbors. It will be marked as inactive and be removed from .

The idea behind our sweeping algorithm is as follows. We advance a sweepline parallel to the segments in the input dexel (Figures 17 and 16). At each step, for each seed , we compute the intersection of the current line with the points in that are at a distance . By choosing the visibility direction to be the same as the sweeping direction, then only the seeds previously encountered in the sweep will contribute to this intersection (upcoming seeds will have an empty Voronoi cell ). Figure 9 shows a Voronoi diagram of points, with two half-space Voronoi diagrams of opposite directions. To make the computation efficient, we do not want to iterate through all the seeds to compute the intersection each time we advance the sweep line (which would make the algorithm in the number of seeds). Instead, we want to keep only a small list of active seeds, that will contribute to the dilated shape and intersect the current sweepline : this will decrease the complexity to , where is the number of segments generated by the dilation process.


From an algorithmic point of view, we maintain, during the sweeping algorithm, two data structures: , the list of active seed segments , whose Voronoi cell intersects the current sweepline (and are at distance from the current sweepline), and , a priority queue of upcoming events, that indicates when an active seed can safely be removed from the as we advance the sweepline, i.e. when it no longer affect the result of the dilation. The events in can be of two types: (1) a Voronoi vertex (the intersection between Voronoi cells); and (2) a seed becoming inactive due to a distance from the sweepline. Voronoi vertices can be located at the intersection of the bisector curves between 3 consecutive seed segments , as illustrated in Figure 17. Beyond this intersection, we know that either or will be closer to the sweepline, so we can remove from (its Voronoi cell will no longer intersect the sweepline).

Both and can be represented by standard STL data structures. Note that since the segments stored in are disjoints, they can be stored efficiently in a std::set<> (sorted by the coordinate of their midpoint). A detailed description of our sweep-line algorithm is given in pseudo-code in Figures 19 and 18. We also provide a full C++ implementation of our method in supplemental. In line 15, the function DilateLine computed the result of the dilation on the current sweep line, by going through the list of active seeds, and merging the resulting dilated segments. Insertion and removal of seed segments in is handled by the functions InsertSeedSegment and RemoveSeedSegment respectively. When inserting a new active segment in , to maintain the efficient storage with the std::set<>, we need to remove subsegments which are occluded by the newly inserted seed segment (and split partially occluded segments). Indeed, the contribution of such segments is superseded by the new seed segment that is inserted. Then, we need to compute the possible Voronoi vertices formed by the newly inserted seed and its neighboring segments in . The derivation for the coordinates of the Voronoi vertex of 3 parallel seed segments in given in Appendix A.

1:2D dexel structure + dilation radius .
2:Half-dilated dexel shape .
3:function VoronoiSweepLine()
4:      Set of active segments on the sweep line
5:      List of removal events marking a seed as inactive
6:      Dilated result
7:     for  do
8:         for all  do
9:              RemoveSeedSegment;
10:         end for
11:         for all  do
12:              InsertSeedSegment;
13:         end for
14:          Dilate and merge active seeds on the current sweepline
15:          DilateLine
16:     end for
17:     return
18:end function
Fig. 18: Sweepline algorithm for computing the half-dilated shape .
Set of active seeds on the sweep line
List of removal events
Dilation radius
Seed segment to insert,
2:Updated list of active seeds and events .
3:function InsertSeedSegment()
4:     RemoveOccludedSegments
5:     SplitPartiallyOccluded
6:      Overlaps are resolved, insert into
7:      After this point, will be empty
8:     while  sequence  do
9:          VoronoiVertex See fig. 17
10:         if  then
11:               Seed becomes inactive
12:         else
13:               Remove later on
14:         end if
15:     end while
16:     while  sequence  do
17:          Repeat operation on the right side of
18:     end while
19:end function
Fig. 19: Insertion of a new seed segment into .

3.3 Half-Space Power Diagram of Points and 3D Offsets

3D Dilation

As illustrated in Figure 14, the 3D dilation process is decomposed in two stages . We first perform an extrusion along the first axis (in red), followed by a dilation along the second axis (in green). Note that the first operation is not exactly the same as a dilation along the first axis (red plane): the segments are extruded, not dilated (they map to a rectangle, not a disk).

To obtain the final dilated shape in 3D, we need to perform a dilation of the intermediate shape , where each segment is dilated by a different radius along the second axis (green plane in Figure 14), depending on its distance from parent seed (red dot in Figure 14). In the case where the first extrusion of produces overlapping segments in , the overlapping subsegments would need to be dilated by different radii, depending on which segment in it originated from. In such a case, where a subsegment has multiple parents , it is enough to dilate by the radius of its closest parent in (the one for which will be the largest). In practice, we store as a set of non-overlapping segments, as computed by our algorithm Figure 18, we a slight modification to the DilateLine function (line 15) to return the extruded seeds on the current sweepline instead of the dilated ones.

2D Power Diagram

In the rest of the discussion, we now focus on computing the 2D dilation of a shape comprised of non-overlapping segments associated with a dilation radius . In this setting, the computation of the power cells of seed segments becomes much more involved, and breaks some of the assumptions of our sweepline algorithm Figure 18. More specifically, our sweepline algorithm for Voronoi diagram makes the following assumptions about : the seeds projected on the sweepline are non-overlapping segments, and the Voronoi cells induced by the active seeds are continuous regions. Once a seed is inserted in , its cell immediately becomes active, and once we reach the first removal event in , it will become inactive and stop contributing to the dilation . For the power cells of segments, the situation is a little bit different, as illustrated on Figure 20: a power cell can contain disjoints regions of the plane. It is not clear how to maintain a disjoint set of seeds in if we need to start removing and inserting a seed multiple time, and this makes the number of cases to consider grow significantly.

Fig. 20: Special case: the power cell of a seed segment (in green) can be a disconnected region of the plane.


Fig. 21: To facilitate the computation of the power diagram in the second stage, the dilation of a segment is separated into a vertical extrusion (left), and the dilation of its two endpoints (right).

Instead, we propose to circumvent the problem entirely by making the following observation. A 2D segment dilated by a radius is in fact a capsule, which can be described by two half-disks at the endpoints, and a rectangle in the middle. We decompose the half-dilated shape into the union of two different shapes: , the result of the half-dilation of each endpoint; and , where each segment is extruded along the dilation axis by its radius . This decomposition is illustrated in Figure 21.

Fig. 22: Sweepline algorithm computing . When inserting the rightmost point in the set of active seeds , it may happen that the power vertex with its two neighbors be located further away along the sweep direction. In this case, the middle point should not be removed from , as its power cell will continue to intersect the sweepline.

Now, the dilated shape is easy to compute with a forward sweep, as there no Voronoi diagram or power diagram involved: we simply remove occluded subsegments, keeping the ones with the largest radius, and removing a seed once its distance to the sweepline is . To compute , we employ a sweep similar to the one described Section 3.2, but now all the seeds are points, which greatly simplifies the calculation of power vertices. Indeed, the power bisectors are now lines, and the power cells are intersections of half-spaces, and in particular convex. The derivations for the coordinates of a power vertex is given in Appendix B. The only special case to consider here is illustrated in Figure 22: when inserting a new seed and updating , the events corresponding to a power vertex do not always correspond to a removal. For example, in the situation illustrated in Figure 22, the power cell of the middle point will continue to intersect the sweepline after it has passed the power vertex, so we should not remove the middle seed from . Fortunately, this case is easy to detect, and we simply forgo inserting the event in in the first place.

3.4 Complexity Analysis

We analyze the runtime complexity of our algorithm for our first 2D dilation operation. For the purpose of this analysis, we assume that the dexel spacing is , so that the dilation radius is given in the same unit as the dexel numbers. The complexity of combining two dexel data structures is linear in the size of the input, since the segments are already sorted. Thus, we focus on analyzing the cost of the forward dilation operation .

Let be the number of input segments, and be the number of output segments in the dilated shape . In the worst case, each input segment generates distinct output segments, and In the worst case, each seed segment is split twice by every newly inserted segment. Since each seed can split in two separate segments at most one element of , so we have at any time that . Moreover, each seed produces at most three events in

(a Voronoi vertex with its left/right neighbors, and the moment it becomes inactive because of its distance to the sweepline). It follows that

as well. Segments in are stored in a std::set<>, thus insertion and removal (lines 9 and 12) can be performed in time. While the line dilation (line 15) is linear in the size of , the total number of segments produced by this line cannot exceed , so the amortized time complexity over the whole sweep is also . This brings the final cost of the whole dilation algorithm to a time complexity that is , and it does not depend on the dilation radius (apart from the output size ). In comparison, the offseting algorithm presented in [Wang:2013:GBO] has a total complexity that grows proportionally to .

For the second dilation operation, where each input segment is associated a specific dilation radius, the result is similar. Indeed, is computed using the same algorithm as before, so the analysis still holds. Combining the dexels in with the results from can be done linearly in the size of the output as we advance the sweepline, so the total complexity of computing is still .

For the 3D case, the result of a first extrusion is used as input for the second stage dilation, the total complexity is more difficult to analyze, as it also depends on the structure of the intermediate result. In a very conservative estimate, bounding the number of intermediate segments by

, this bring the final complexity of the 3D dilation to , where is the size of the output model. In practice, a lot of segments can be merged in the final output, especially when the dilation radius is large, and may even be smaller than (e.g. when details are erased from the surface).

Finally, we note that in each stage of the 3D pipeline, the 2D dilations can be performed completely independently in every slice of the dexel structure, making the process trivial to parallelize. We discuss in Section 4 the experimental performances of our a multi-threaded implementation of our 3D dilation operation.

4 Results

We implemented our algorithm in C++ using Eigen for linear algebra routines, and Intel Threading Building Blocks for parallelization. We ran our experiments on a desktop with a 6-core Intel® Core™ i7-5930K CPU clocked at 3.5 GHz and 64 GB of memory. The reference implementation, equipped with scripts to reproduce all our results, is attached as addition material and will be released as an open-source project to facilitate the adoption of our technique.

Baseline Comparison

We implemented a simple brute-force dilation algorithm (on the dexel grid) to verify the correctness of our implementation, and to demonstrate the benefits of our technique. In the brute-force algorithm, each segment in the input dexel structure generates an explicit list of dilated segments in a disk of radius around it, and all overlapping segments are merged in the output data structure. Figure 23 compares the two methods using a different number of threads, and with respect to both grid size and dilation radius. In all cases, our algorithm is superior not only asymptotically, but also for a fixed grid size or dilation radius. Since each slice can be treated independently in our two-stage dilation process, our algorithm is embarrassingly parallel, and scales almost linearly with the number of threads used.

We note that the asymptotic time complexity observed in Figure 23 agrees with our analysis in Section 3.4. Indeed, the dexel grid has a number of dexel is proportional to the squared grid size , while the (absolute) dilation radius grows linearly with the grid size. Since the complexity analysis in Section 3.4 uses a dilation radius expressed in dexel units, the observed asymptotic rate of for indicates that our method is indeed .

Fig. 23:

Total running time averaged on 11 testing models, compared with a direct brute-force implementation. Standard deviation is shown in overlay, and convergence rate

are reported in the log-log plots. Radii are relative to the grid size. The top row uses a relative radius of , and the bottom row uses a grid size of .

Dilation vs Erosion

In Figure 24, we compare the performance of the dilation and erosion operator on a small dataset of 11 models, provided in the supplemental material. The dilation operator has a higher cost than the erosion, both asymptotically and in absolute running time. This is likely due to the fact that the erosion operator reduces the number of dexels in the data structure, leading to a clear speedup.

Fig. 24: Total running time of the dilation vs erosion operation, for two different radii. Standard deviation shown in overlay and convergence rate reported in the log-log plot.
(a) Input
(b) Dilation
(c) Erosion
(d) Opening
(e) Closing
Fig. 30: Topological cleaning of an input shape via morphological operations, using a grid of dexels.
(a) Input Surface
(b) Grid Size
(c) Grid Size
(d) Grid Size
(e) Grid Size
(f) Grid Size
Fig. 37: Result of an erosion operation using different grid resolutions.

Comparison with Wang:2013:GBO

The most closely related work on offset from ray-reps representations is Wang:2013:GBO. It proposes to perform an offset from a LDI sampled from three orthogonal directions. The offset is computed as the union of spheres sampled at the endpoints of each segment from all three directions at once. In contrast, our method relies on a single dexel structure, which has both advantages and drawbacks: it is applicable in a situation where only one view is available, or is enough to describe the model (e.g. modeling for additive manufacturing [Lefebvre:2013:IAG]), but it will be less precise on the orthogonal directions (where the 3-views LDI will have more precise samples).

To compare our method with theirs, we matched the parameters of the experiment reported in Table 2 of [Wang:2013:GBO], and compare our results with the 4-thread CPU version reported in [Wang:2013:GBO]. At a resolution of and relative radius of the diagonal of the bounding box, our running time was 3.4s, 3.7s, and 2.7s for the Vase-lion, Filigree, and Buddha models respectively. The timings reported in [Wang:2013:GBO] for their CPU implementation are over 120s, suggesting that our CPU implementation is competitive even with their GPU implementation. Extending our method to the GPU will be challenging, but could enable real-time offsetting on large and complex dexel structures, which we believe would be an interesting venue for future work.

Topological Cleaning

Our efficient dilation and erosion operators can be combined to obtain efficient opening and closing operators (Figure 30). For example, the closing operation can be used to remove topological noise, i.e. small handles, by first dilating the shape by a fixed offset, and then partially undoing it using erosion. While most parts of the shape will recover their original shape, small holes and sharp features will not, providing an effective way to simplify the shape topology.


The compactness of the dexel representation enables us to represent and process immense volumes on normal desktop computers. An example is shown in Figure 37 for the erosion operation. Note that the results on the right have a resolution sufficiently high to hide the dexels: this resolution would be prohibitive with a traditional boundary or voxel representation.

3D Printing

The Boolean difference between a dilation and erosion of a shape produces a shell (of controllable thickness). This operation is common for 3D printing applications, since the interior of an object is usually left void (or filled with support structures) to save material. Another typical use-case for creating thick shells out of a surface mesh is the creation of molds [Malomo:2016:FAD]. We show a high resolution example in Figure 4, which has been 3D printed using PLA plastic on a Ultimaker 3 printer.

5 Future Work and Concluding Remarks

Our current implementation is restricted to uniform morphological operations, and it would be interesting to extend it to single direction thickening (for example in the normal direction only) or to directly work on a LDI offset, i.e. representing the shape with 3 dexel representations, one for each axis. A GPU implementation of our technique would likely provide a sufficient speedup to enable real-time processing of ray-reps representation at the resolution typically used by modern 3D printers.

To conclude, we proposed an algorithm to efficiently compute morphological operations on ray-rep representations, targeting in particular the generation of surface offsets. Beside offering theoretical insights on power diagrams and their application to surface offsets, our algorithm is simple, robust, and efficient: it is an ideal tool in 3D printing applications, since it can directly process voxel or dexel representations to filter out topological noise or extract volumetric shells from boundary representations.


This work was supported in part by the NSF CAREER award IIS-1652515, a gift from Adobe, and a gift from nTopology.

Zhen Chen is working towards his bachelor degree in Computational & Informational Science at University of Science and Technology of China (USTC), China. He will graduate at the beginning of June 2018. From June to August 2017, he has been a visiting student at the Courant Institute of Mathematical Sciences (New York University, USA). His research interests are 3D printing and geometry processing.

Daniele Panozzo Daniele Panozzo is an assistant professor of computer science at the Courant Institute of Mathematical Sciences in New York University. Prior to joining NYU he was a postdoctoral researcher at ETH Zurich (2012-2015). He earned his PhD in Computer Science from the University of Genova (2012) and his doctoral thesis received the EUROGRAPHICS Award for Best PhD Thesis (2013). He received the EUROGRAPHICS Young Researcher Award in 2015 and the NSF CAREER Award in 2017. Daniele is leading the development of libigl (, an award-winning (EUROGRAPHICS Symposium of Geometry Processing Software Award, 2015) open-source geometry processing library that supports academic and industrial research and practice. Daniele is chairing the Graphics Replicability Stamp (, which is an initiative to promote reproducibility of research results and to allow scientists and practitioners to immediately benefit from state-of-the-art research results. Daniele’s research interests are in digital fabrication, geometry processing, architectural geometry, and discrete differential geometry.

Jérémie Dumas Jérémie Dumas is a postdoctoral fellow at the Courant Institute of Mathematical Sciences in New York University. Prior to joining NYU he completed his PhD at INRIA Nancy Grand-Est (2017), under the direction of Sylvain Lefebvre. His doctoral thesis received the EUROGRAPHICS Award for Best PhD Thesis (2018). His work focuses on shape synthesis for digital fabrication, including shape optimization, simulation, microstructures, and procedural synthesis.

Appendix A Voronoi Vertex between Three Parallel Segments in 2D

The bisector of two parallel segment seeds is 2D is a piecewise-quadratic curve, as illustrated in Figure 17. A Voronoi vertex is a point at the intersection of the bisector curves between three segment seeds. Because the segment seeds are non-overlapping, the Voronoi vertex can be either between three points (Section A.1) or between one segment and two points (Section A.2).

a.1 Voronoi Vertex between Three Points

Let be three points, with coordinates . A point lying on the bisector of satisfies


After simplification, we get


Similarly, for a point lying on the bisector of ,


We can get the coordinate of Voronoi vertex between three points by solving the system formed by Equations 9 and 8:


a.2 Voronoi Vertex between a Segment and Two Points

Let be three seeds. We only need to consider the case where , other cases are similar. For the Voronoi vertex to intersect the bisectors where it is closest to the interior of , and not one its endpoints, we need to have . It follows the distance from to the segment is:


By definition of the Voronoi vertex,


Developing Equation 12, we get


Similarly, developing Equation 13 leads to


Now, let

We can rewrite the system of equations as


Solving Equation 17, if the roots exist, we will have

Choosing the solution that belongs to , and substituting into , we will get the -coordinate of our Voronoi vertex.

Appendix B Power Vertex between Three Points in 2D

Let be three seeds. A power vertex can be computed from the intersection of two bisector lines.

A point lying on the bisector of satisfies:


A similar equation holds for the bisector of . This translates into the following system of equations: