Support-Free Hollowing for 3D Printing via Voronoi Diagram of Ellipses

08/22/2017 ∙ by Mokwon Lee, et al. ∙ HanYang University USTC 0

Recent work has demonstrated that the interior material layout of a 3D model can be designed to make a fabricated replica satisfy application-specific demands on its physical properties such as resistance to external loads. A widely used practice to fabricate such models is by layer-based additive manufacturing (AM) or 3D printing, which however suffers from the problem of adding and removing interior supporting structures. In this paper, we present a novel method for generating support-free elliptic hollowing for 3D shapes which can entirely avoid additional supporting structures. To achieve this, we perform the ellipse hollowing in the polygons on parallel section planes and protrude the ellipses of one plane to its neighboring planes. To pack the ellipses in a polygon, we construct the Voronoi diagram of ellipses to efficiently reason the free-space around the ellipses and other geometric features by taking advantage of the available algorithm for the efficient and robust construction of the Voronoi diagram of circles. We demonstrate the effectiveness and feasibility of our proposed method by generating interior designs for and printing various 3D shapes.



There are no comments yet.


page 3

page 8

page 12

page 15

page 20

page 21

page 22

page 25

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Additive manufacturing (AM), also called 3D printing, has been advanced rapidly to make customized 3D models. Comparing to traditional manufacturing techniques, it offers enormous geometrical freedom for designers to create highly optimized components with various functionality.

Model hollowing is a typical practice for purposes of reducing printing material and time in 3D printing of light-weighted artifacts and various methods on generating optimized interior have been developed during the last few years Stava et al. (2012); Wang et al. (2013); Lu et al. (2014). The mainstream of 3D printing technologies, such as Fused Deposition Modeling (FDM) and Stereolithography (SLA), requires additional supporting structures to avoid the falling of relatively large overhanging parts during the printing process Strano et al. (2013); Dumas et al. (2014); Vanek et al. (2014). Generally the extra supporting structures have to be removed either manually or by dissolving dissolvable material from the printed objects.

However, there is no way at all to remove the supporting material inside interior voids of a printed object. A naïve method is to first decompose the model into a few subparts, print them and remove the supporting material individually, and then glue them together. Obviously it may largely affect the computed physical properties. Recently, there are quite a few attempts to create support-free interior voids or structures by constraining boundary slopes Langelaar (2016); Wu et al. (2016a); Reiner and Lefebvre (2016).

A noticeable work adopts rhombic cell structures, where the slope angles of all rhombic cells are smaller than a prescribed maximum overhang-angle, as a infill pattern to generate support-free interior voids inside the objects Wu et al. (2016a). However, the rhombic cells have only boundaries, which suffers serious problem of stress concentration Pilkey and Pilkey (2008). The stress around a discontinuity, e.g., a corner, will be excessively high when compared to the stresses at the other smooth areas as shown in Figure LABEL:fig:stress-concentration. For example, hatches and doorways in airplanes are oval to stay away from being broken easily Karthikeyan (2016)

. Rounded corners are structurally more beneficial than sharp corners and also reduce the probability of crack development unlike sharp corners.

stress-concentration0.8 Stress analysis on two thin plates by conducting two external loads (represented by two arrows respectively) on them. There is a discontinuity in the concave region of the upper shape while the concave region of the lower shape has a smooth boundary. From the stress maps, it is seen that the stress around the discontinuity is excessively high.

Our Work

To this end, we present a novel method for hollowing 3D shapes with support-free smooth elliptic interior voids. To make it simple, we first derive a class of support-free ellipses in 2D case, based on the observation on sticky property of printing material. Then we develop a novel approach for packing these support-free ellipses in the interior of 2D shapes, which is a very challenging problem. To achieve this goal, we develop a greedy but efficient algorithm on adding these ellipses successively in 2D shapes via the Voronoi diagram of ellipses in polygons using that of circular disks. Then the hollowed ellipses are extruded into 3D volume and thus a support-free hollowed 3D shape is generated. Various experimental results have demonstrated the feasibility and applicability of our proposed method.


Our contributions are summarized as follows.

  • We develop efficient algorithms for computing the Voronoi diagram of polygon and that of ellipses within a polygon;

  • We develop an efficient method for packing 2D ellipses with derivation of support-free constraint;

  • We propose a method for generating support-free elliptic voids for 3D shapes.

To the best of our knowledge, our work is the first to offer a framework to generate a smooth, elliptic support-free interior voids for printing 3D shapes via an efficient computational scheme. This provides a practically feasible operator for support-free hollowing of 3D shapes and exemplifies research along this direction.

2 Related Work

3D printing technology has drawn a lot of attention in geometric and physical modeling and optimization in computer graphics community. We discuss about the work closely related to our study and give specific discussion about their strength and limitations.

Supporting Structures for 3D Printing

Extra supporting structures are required to make 3D shapes printable for shapes with large overhanging parts, which leads to material waste and longer printing time. Some methods, which adopt various supporting structures such as scaffold-like structures Dumas et al. (2014) and tree-like ones Vanek et al. (2014), have been developed to generate economic usage of supporting structure for 3D models. Vanek et al. Vanek et al. (2014) searches an optimal printing direction by reducing the total area of facing down regions where require additional supporting structures while Zhang et al. Zhang et al. (2015) develops a training-and-learning model to determine the optimal printing direction considering multiple factors such as contact area, viewpoint preference, and visual saliency. The work of Hu et al. (2015) optimizes the shape of an input model to reduce the area of facing down regions for less supporting structures. However, adding supporting structure for interior voids will suffer the serious problem that it is impossible to remove these extra structures without breaking the object into pieces. Instead, we generate interior support-free voids, which completely avoids the usage of extra supports.

Interior Hollowing

Significant work has been done in generating interior structure of a model to meet various geometric and/or physical properties. Stava et al. Stava et al. (2012) hollows a 3D-printed object while maintaining its structural strength by adding some internal struts. The interior is optimized by a reduced-order parameterization of offset surfaces in Musialski et al. (2015). Various internal structures, such as the skin-frame structure Wang et al. (2013), the honeycomb-like structure Lu et al. (2014), and the medial axis tree structure Zhang et al. (2015), were developed for cost-effective purposes while preserving the structural strength of printed objects. Both static balance and dynamic balance have been studied by designing the interior infills as well as changing the model shapes Prévost et al. (2013); Christiansen et al. (2015); Bächer et al. (2014). Instead of designing hollowing structure explicitly, a lot of efforts have been put on topology optimization to obtain distributions of material according to certain performance criteria during the last three decades Deaton and Grandhi (2014); Wu et al. (2016b). However, these works have not handled the problem of avoiding large overhangs. We study this problem by developing a carving operator via support-free elliptic voids.

Support-Free Structure

Hu et al. Hu et al. (2014) proposes a method to decompose a 3D object into support-free pyramidal subparts. Reiner and Lefebvre Reiner and Lefebvre (2016) proposes an interactive sculpting system for designing support-free models. Recent attempts have been put on creating support-free interior structure for 3D printing. Wu et al. Wu et al. (2016a) develops a method to generate support-free infill structures on adaptive rhombic cells. A concurrent work of Langelaar Langelaar (2016) considers an overhang angle threshold in topology optimization and generates a support-free material distribution. However, these works only involve overhang angle and generate boundaries inevitably, resulting in large stress concentration in discontinuities. In this paper, we develop a special class of support-free ellipses and adopt them as an interior carving operator to create infill structures.

Ellipse Packing

Our method of ellipse carving is quite related to the problems of ellipse and ellipsoid packing which are NP-hard. In science, the problem was prevalently approached from the packing ratio, particle size distribution, or jamming point of view to understand material properties by enforcing contacts between ellipses and using Monte Carlo method Buchalter and Bradley (1992); Schreck et al. (2010), Molecular Dynamics Donev et al. (2004, 2007), local and greedy algorithms Hit (2013); Birgin et al. (2016) with sampling points on the ellipse boundary Delaney et al. (2005) in either regular containers or arbitrary domains Lozano et al. (2016). However, these methods do not balance the computation cost and packing density outcome well because research interests were primarily on the discovery of new phenomena. This is quite different from the circle packing problem which have been extensively studied from the view points of both solution quality (i.e. packing ratio) and computation time, from early study using an event-driven algorithm with a bucket acceleration Jodrey and Tory (1985); Lubachevsky and Stillinger (1990); Lubachevsky (1991) to recent ones using a systematic reasoning of empty space Specht (2015); Sugihara et al. (2004). However, the infilled ellipses in our work have special constraints which make current methods infeasible to use. In this work, we develop a new method for packing ellipses in arbitrary shapes with the mathematical tool of Voronoi diagram (VD). Specifically, we perform an ellipse packing in the polygons on parallel section planes of a 3D model and protrude the ellipses of one plane to its neighboring planes in the shape volume to generate the hollowed results. To better pack the ellipses in a polygon, we construct the Voronoi diagram of ellipses to efficiently reason the free-space around the ellipses and other geometric features by taking advantage of the available algorithm for the efficient and robust construction of the Voronoi diagram of circular disks which approximate the ellipses Lee et al. (2016). The algorithm inherits the disk packing algorithm using the VD of disks Sugihara et al. (2004).

3 Notations and Overview

For the sake of simplicity, we first elaborate on our method in a 2D setting. The extension to 3D is realized by extrusion in Section 5.

3.1 Support-Free Ellipses

Fabrication and Material Parameters

Denote as the printing precision, i.e., the thickness of each fabrication layer, which is for general FDM printers. Due to the sticky property of printing material, a short length of horizontal hangover can be successfully printed without extra supporting structure. Denote as the maximally allowed overhang-angle. The material-dependent parameters and can be measured by experiments, for example, and for plastic PLA material. The minimum wall thickness is set as .

Support-Free Ellipses

Denote and as the horizontal axis and vertical axis of an ellipse , respectively. A hollowed ellipse is called support-free if it can be printed without any extra supporting structure. Given an ellipse shown in Figure LABEL:fig:Ellipse-SupportFree(a), denote and as the two points whose tangent lines have an overhang angle of . By many experiments we have proved the following observation: if the distance between and is no larger than , then this ellipse is support-free, that is, the elliptic arc between and (shown in red) can be safely printed without any extra support. On the other hand, too small interior ellipses cannot be printed and thus we set a lower bound of as . Based on the observations, we have derived the conditions for a support-free ellipse as:


It is worthwhile to mentioning that a uniform shrinking of a support-free ellipse is still support-free as the support-free conditions in Equation (1) are convex.

Ellipse-SupportFree0.350.35 (a) The condition of a support-free ellipse: where and are the two points with tangent lines of overhang angle ; (b) Ellipses with and within the gray region are support-free.

VD-illustration0.470.470.470.47 Illustration of VDs in this study. (a) The VD of a polygon (); (b) The VD of a disk set within (); (c) The VD of the disk approximation of (); (d) The VD of ellipses within after the V-faces belonging to a P-edge are merged ().

idea-overview-5-bunnies0.380.380.380.380.38 The ellipse packing process for the bunny polygon using the VD of ellipses. (a) The VD of boundary disks. (b) The VD of the bunny interior. (c) The in-disks of the first ellipse within the clearance probe are incremented. (d) The fifth ellipse located within the clearance probe after the four ellipses were incremented. (e) The 100-th ellipse to be incremented.

3.2 Ellipse Packing Method


Given a polygon , our goal is to find an optimal hollowing, i.e., packing, of a set of ellipses within so as to maximize the sum of the areas all ellipses (the packing ratio) while satisfying mechanical and physicochemical constraints.

2D Polygon

Given a 3D mesh representing the boundary of an oriented solid, possibly with handles and interior voids, and a fabrication orientation (-axis), we project onto planes parallel to -axis and choose the projection direction with the largest projected silhouette area. Then is represented as a sequence of parallel cross-sections, i.e., 2D polygonal shapes , along the chosen projection direction. We choose the polygon in with the largest area and denote it as which may have internal holes (voids). Denote where and are the sets of vertices and edges, thus denoted as P-vertices and P-edges, respectively. If contains holes inside, its boundary is represented as a few closed loops: One outer loop (with counter-clockwise vertices) and a few inner ones (with clockwise vertices). There can be more than one polygon on a section plane. In the following, we present our method on hollowing with support-free ellipses, i.e., packing of support-free ellipses in .

Voronoi Diagrams (VD)

We use Voronoi diagrams (VD) to compute the hollowed ellipses because it is the most efficient and compact data structure for spatial reasoning among particles. The VD of a generator set is the tessellation of space such that each cell of the tessellation consists of the locations closer to a corresponding generator than to others. Various types of VD can be defined by generalizing the generator type, the distance definition, and the dimension. For details, see Okabe et al. (1999).

Voronoi Diagrams of Disks and Ellipses

In this study, we use the VDs of disks and ellipses within a polygon with the Euclidean -distance in 2D. See Figure LABEL:fig:VD-illustration: (a) The VD of the polygon (); (b) The VD of a disk set within (); (c) The VD of the disk approximation of (); (d) The VD of ellipses within the approximation (). Note that both and can be correctly, efficiently, and robustly computed. However, we compute instead of because of the challenges involved in the computation of V-vertices and V-edges of , which will be explained in detail in Section 5. Hence, we construct the approximation instead of where each P-vertex is associated with a disk called at-disk (the red filled-circles in Figure LABEL:fig:VD-illustration(c,d)) and each P-edge is associated with more than two disks called on-disks (the blue filled-circles). In Figure LABEL:fig:VD-illustration(d), the blue thick circle is the maximum clearance probe , centered at a V-vertex of , that guarantees its inscribing ellipse is intersection-free with any other existing ellipses and the red thick circle is 70% shrunken probe in which the actual ellipse is created. The shrinking ratio is given by users according to their intention to control the overall distribution of ellipse heights. Be aware that is associated with a disk set where and are the sets of at-disks and on-disks, respectively.

Idea of the Packing Algorithm

Let be the VD of the interior of (Figure LABEL:fig:idea-overview-5-bunnies (b)) obtained by trimming the exterior part of the entire VD in Figure LABEL:fig:idea-overview-5-bunnies (a) and that of within . Let be the biggest empty circle centered at a V-vertex of where its radius is given by the distance from to the boundary of its generators. Let be the possible biggest one, called the maximum clearance probe, among all such empty circles from the V-vertices of . Let be the V-vertex corresponding to . Starting from , we first find the V-vertex with (the blue solid circle in Figure LABEL:fig:idea-overview-5-bunnies(c)) and its shrunken probe (the red one) and place one new ellipse within (Figure LABEL:fig:VD-illustration(d) and LABEL:fig:idea-overview-5-bunnies(c)). Then, we construct by inserting into . We repeat the clearance-probe-finding, the ellipse-placement, and Voronoi diagram update processes for a sufficient number of times until a termination condition is met, as shown in Figure LABEL:fig:idea-overview-5-bunnies(d,e).

4 Ellipse Hollowing via Voronoi Diagram


There are two major computational phases for using VD in the ellipse hollowing: first, the construction of and second, the construction of .

4.1 Voronoi Diagram of a Polygon


It is well-known that the algorithm for an efficient and robust construction of is not trivial due to the influence of numerical error on maintaining correct topology of Voronoi diagram Kim et al. (1995); Okabe et al. (1999). In this study we developed and implemented a new, simple, and efficient yet robust algorithm based on the topology-oriented approach Sugihara and Iri (1992); Held (2001); Held and Huber (2009); Lee et al. (2016) for . This is because we eventually need to construct which is lacking in existing codes such as CGAL CGAL (2016); Alliez et al. (2007) and VRONI Held (2001).

TOI-D Algorithm

The proposed algorithm takes advantage of the Voronoi diagram of circular disks Kim et al. (2001a, b), particularly the recently reported topology-oriented incremental (TOI) algorithm for computing the Voronoi of circular disks, thus abbreviated as the TOI-D algorithm, which takes time in the worst case but time on average for disks Lee et al. (2016). The idea is to approximate target geometric entities using circular disks in a sufficient resolution, construct the VD of the disks using the TOI-D algorithm, and merge some V-cells. While a similar idea was used to curved objects using the ordinary Voronoi diagram of points which were sampled from curves Sugihara (1993); Held (2001); Emiris et al. (2013), the proposed algorithm using the TOI-D algorithm is much powerful as circles can significantly reduce problem size and complexity.

Approximating Polygon with Disks

We represent a polygon as follows. Let be the shortest P-edge with its length . Suppose that we cover with two open disks with the diameter by placing their centers on and the boundary of each disk coincides either one of the two extreme points of . For each of the other P-edges with the length , we place non-overlapping open disks in a sequel on the P-edge ( denotes a floor function). The uncovered remaining segment with the length on the P-edge is then covered by one, and only one, smaller open disk. We place this smaller disk in the middle of the P-edge for algorithmic simplicity. Hence, the P-edge is entirely covered by mutually exclusive open disks called an on-disk (The blue ones in Figure LABEL:fig:VD-illustration(c)). A disk is a child of an associated P-edge and the parent of . We also place a same size disk, called at-disk, at each P-vertex (The red ones in Figure LABEL:fig:VD-illustration(c)): and also have a child-parent relationship. Each disk knows its parent and each parent knows its children disks via pointers. We eventually have a set of children disks representing where no disk contains any other while two disks may intersect, .

The computation speed of the idea above, which is sensitive to the shortest P-edge, can be improved by enforcing fewer on-disks for each P-edge. We initially allocate only three on-disk: two near the extreme points of the P-edge with the radii identical to the at-disks for P-vertices and the third with the diameter covering the entire rest segment of the P-edge. If the third on-disk intersects any other disk, either an at-disk or on-disk, we subdivide it (to avoid computational complications) until the intersection is resolved by employing a bucket system to accelerate the intersection check. In this way, the number of children disks can be significantly reduced. This approach of using on-disks with non-uniform sizes works well for the polygons produced from fine mesh models such as bunny.

However, there are models such that for the the minimum wall thickness . In such a case, the proposed algorithm both leaves an undesirable bulk material in printed artifacts and may run into a computational complication and thus we subdivide each P-edge into multiple shorter P-edges with . Some engineering models with may have frequently large planar facets possibly with smaller ones. We also subdivide long P-edges into several short ones with the same rule above.

VD-illustration-further0.440.2160.440.2160.440.216 Three important steps of the TOI-P algorithm for . (a,b) after the merge process. (c,d) After some V-vertices are relocated. (e,f) After outside V-edges are trimmed, the V-vertex coordinates and V-edge equations are computed.

TOI-P Algorithm

We construct of the disk set where and using the TOI-D algorithm (Figure LABEL:fig:VD-illustration(c)) and merge the V-cells of the children on-disks of each P-edge (Figure LABEL:fig:VD-illustration-further(a,b)). As shown from the figure, the resulting VD structure has a unique V-cell for both each P-edge and each P-vertex and thus its structure is close to from both topology and geometry point of views. However, some V-vertices are off P-vertices which play the role of V-vertices (Figure LABEL:fig:VD-illustration-further(b)). Hence, they should be moved to the related P-vertices. Figure LABEL:fig:VD-illustration-further(c,d) shows the VD after relocating those V-vertices to the polygon boundary and flipping some V-edges to get the correct topology, both taking at most time for disks. Note that a V-edge flipping is required to get the correct topological structure of (Compare Figure LABEL:fig:VD-illustration-further(b) and (d)). Then, removing the exterior part of the VD and computing the V-vertex coordinates and the V-edge equations transforms the intermediate VD structure to the correct (Figure LABEL:fig:VD-illustration-further(e,f)). Note that some previously linear V-edges are curved. Removing the at-disks and on-disks results in the VD of Figure LABEL:fig:VD-illustration(a).

There are three cases of V-edges: i) If a V-edge is defined between two P-vertices, is a line segment; ii) Between two P-edges, is also a line segment; iii) Between a P-vertex and P-edge, is a parabolic arc. As the V-edges of are quadratic curve segments, they can be represented as a rational quadratic Bézier curve Kim et al. (1995). There are four cases of V-vertices: i) Among three line segments; ii) Among two line segments and one point; iii) Among one line segment and two points; iv) Among three points. Each V-vertex coordinate and V-edge equation can be correctly computed in time Kim et al. (1995).

Time Complexity

The TOI-D algorithm for the construction of the VD of disks takes time on average and time in the worst case Lee et al. (2016). With the polygon of P-vertices and P-edges, , the merge process takes time because each merge of two adjacent V-cells takes time. As V-vertex shift and coordinate correction takes time, the entire process for takes time provided that is given. As depends either on the length of the shortest P-edge or on the minimum wall thickness, the proposed TOI-P algorithm is input-sensitive. All time complexities in this paper are in the worst case sense unless otherwise stated.

4.2 Voronoi Diagram of Ellipses in a Polygon

Regarding for as , we first find the V-vertex with the maximum clearance probe and place an ellipse which inscribes so that its center coincides . We incrementally insert at of to get . Instead of directly inserting into , we insert each in-disk, one by one, into where . As the ellipse increment process goes on, contains all the in-disks of the ellipses incremented so far in addition to .


Construction of involves the computation of V-vertices and V-edges and the topological structure among them where each of these tasks is challenging. Consider a V-vertex defined by three ellipses. It is known that there can be up to 184 complex circles that are simultaneously tangent to three conics in the plane and each corresponds to a root of a polynomial of degree 184 Emiris and Tzoumas (2005). It is hard to expect to find the roots of a polynomial of such a high degree both exactly and efficiently. Given 10-bit precision to represent the coefficients of three random ellipses, each coefficient of the resultant necessary for the exact computation of a V-vertex is, on average, 4603-bit integers Emiris and Tzoumas (2005). Hence, the exact and efficient computation of the correct coordinate of itself is not an easy problem at all. We are not aware of any method to solve this resultant exactly and efficiently and thus an exact computation approach to construct the VD of ellipses seems impractical. The V-edge between two ellipses can be more complicated than one might expect. Even the bisector between a point and an ellipse can be very complicated Farouki and Johnstone (1994): It may have cusps and self-intersections and is disconnected if the point is located outside the ellipse. The bisector between two rational curves can be non-rational and even a two-dimensional object Alt et al. (2005). Therefore, the computation of V-vertices, V-edges, and the their association through the topological structure among ellipses in a free-space is a challenge, not to mention about the ellipses within a polygon. As far as we know, no study has been reported for constructing .

Approximating an Ellipse with Circles

We represent an ellipse as an approximation with a set

of an odd number of disks, called in-disks which inscribes

(The yellow ones in Figure LABEL:fig:VD-illustration(d) and Figure LABEL:fig:idea-overview-5-bunnies(c)). In-disks may intersect but none is contained by another. The in-disks are generated as follows. The first in-disk is the maximal inscribing disk which is centered at the center of . Let be an approximation error defined as the horizontal distance between and . Then, a point can be located for an a priori defined error, say . Hence, a second in-disk passes through while inscribing . We alternate this calculation up and down of to get and , respectively. Repeating this calculation produces in-disks with a strictly controlled error bound . Given , a shorter ellipse has fewer in-disks than a longer one does.

TOI-EinP Algorithm

The idea is very simple as follows. We insert each in-disk in into of existing disks and then merge the V-cells of the in-disks of . If a sufficient number of in-disks approximates each ellipse, the VD of the disks well-approximate the VD of ellipses from both topology and geometry point of views. As ellipses do not intersect, the in-disks from distinct ellipses do not intersect and the in-disks do not intersect both on-disks and at-disks on the polygon boundary, either.

The increment of an in-disk is done using the TOI-D algorithm. When we increment an in-disk, we maintain a dual representation of both and . In other words, and are carefully synchronized in the following sense. We first incrementally update until all in-disks of in (thus those of ) are exhausted to get . Then, we merge the V-cells of the in-disks of to produce the topology of . After the merge process, each of the remaining V-vertices of becomes the V-vertices of and a connected subset of some appropriate remaining V-edges becomes the V-edge of . Hence, we carefully maintain the correspondence of the V-vertices and V-edges between and . Note that the merge can also be done incrementally as soon as after an in-disk is incremented. Figure LABEL:fig:sequential-illustration-hollow-6-dogs shows the process of incrementing ellipses: (a) The clearance probe (the large blue circle), its shrunken probe (the red circle), the ellipse within , and the biggest in-disk (blue filled circle) is incremented into the VD; (b) The second in-disk (the blue filled circle) is incremented into the VD (The previously incremented in-disk is yellow now); (c) After four in-disks are incremented; (d) After all in-disks of the first ellipse are incremented; (e) After the second ellipse is incremented; (f) After the third ellipse is incremented.

The V-vertices of remaining after the V-cell merge have their coordinates inheriting from which are computed from a triplet of in-disks. Thus, they are not necessarily correct for ellipses and it is necessary to compute their correct coordinates for the successful packing of next ellipse because the maximum clearance probe needs to be found from these V-vertices. The six cases of generator combination for a V-vertex in among ellipses, line segments (i.e. P-edges), and points (i.e. P-vertices) becomes a unified case of among three ellipses in because of and . The six cases are as follows: among three ellipses, among two ellipses and one line segment, among two ellipses and one point, among one ellipse and two line segments, among one ellipse and two points, and among one ellipse, one line, and one point.

sequential-illustration-hollow-6-dogs0.3840.3840.3840.3840.3840.384 The ellipse increment process of the TOI-EinP algorithm through the increments of in-disks into the VD structure. (a) Identification of the clearance probes, the ellipse, and the increment of the biggest in-disk. (b) the increment of the second in-disk. (c) the increment of the fourth in-disk. (d) after the increment of all in-disks of the first ellipse. (e) after the increment of the second ellipse. (f) the increment of the third ellipse.

V-vertex Coordinate among Three Ellipses

Consider a V-vertex defined by three ellipse generators. We iteratively find the correct location of in starting with its initial coordinate provided by . In other words, is initially equidistant from three in-disks where each is a child of each of three ellipses. We project to each of the three ellipses to find its footprint (which is the closest location on an ellipse from ). Then, we compute the circumcircle, say , which passes through the three footprints and use the center of as the new coordinate of . Provided that each ellipse is approximated by a sufficient number of in-disks, the iteration of this footprint-projection and circumcircle-finding process quickly converges to the correct coordinate of due to the convexity of ellipse. Experiment shows that the initial coordinate of is already very close to the converged coordinate. The situation that a generator(s) of is at-disk or on-disk can be handled as an easier special case

V-edge Between Two Ellipses

Due to the current theoretical limitations, it is practically inevitable to approximate a V-edge with a sequence of passing points computed through the envelopes of families of point/curve bisectors Farouki and Johnstone (1994, 1992) or a sequence of curve segments Farouki and Ramamurthy (1998). There exist other approaches to trace bisectors for VD and medial axis transformations Omirou and Demosthenous (2006); Cao and Liu (2008); Cao et al. (2009).

We emphasize that the topological structure of is already known. In other words, for each V-edge , its starting and ending V-vertices are known with correct coordinates along with its two elliptic generators. We approximate each V-edge as a sequence of points by tracing the V-edge in a way conceptually similar to the tracing algorithm of the intersection curve between two free-form surfaces Barnhill et al. (1987). Tracing V-edges in this study is, however, much simpler than tracing general intersection curve in that i) the coordinates of two V-vertices of each V-edge are known, ii) V-edges are planar, and iii) each V-edge is -continuous between two V-vertices. The case that is defined between one ellipse and one at-disk (or on-disk) is an easier special case.

Finding Footprints

Finding footprints is a key building block. Suppose that is a point outside an ellipse and is a line passing through . It is known that there are four, three, or two locations on that perpendicularly intersects depending on whether lies inside the evolute, lies on the evolute but not at a cusp, or lies on a cusp or outside the evolute, respectively Emiris and Tzoumas (2005). Finding the perpendicular intersection between and can be formulated as a root-finding problem of a quartic polynomial thus taking time. The footprint of is obviously one of these locations which determines the minimum distance.

Time Complexity

Suppose that has elements of and ellipses and suppose that has children disks and there are in-disks for all incremented ellipses so far. Given a new ellipse , with in-disks in , the increment of all in-disks into takes time. Then, it is followed by the merges of the V-cells among the in-disks of taking time. This completes the computation of the topological structure of . Then, the computation of the geometry of each of V-vertices and V-edges takes time. Hence, given the synchronized and , the increment of an ellipse takes time. As and , the increment of one ellipse takes time in the worst case. Hence, the increment of ellipses takes time.

4.3 Implementation Issues

Criteria of Adding Ellipses

During the ellipse increment, the distance between two adjacent interior ellipses should be no less than the minimum wall thickness . An intuitive scheme is to maximally pack the polygon with ellipses and shrink each ellipse by . As an ellipse is not offset-invariant, the shrunk ellipse needs to be approximated but can be effectively computed.

We instead use a computationally easier yet equally effective scheme as follows. We have two parameters to control the height of each ellipse: the minimum wall thickness and the shrink ratio of the maximal clearance probe . Given , we multiply to to get a shrunken probe with the shrunken radius . If , we use to produce an inscribing ellipse and increment in the VD. Otherwise, we reduce the radius of by and produce the inscribing ellipse. The purpose of is to provide users a convenient handle to control the overall distribution of ellipse heights because, in addition to the NP-hardness of the optimal ellipse packing, we never know which way is best even if we only consider geometry. Moreover, experienced users may want to have a control of overall shape distribution by tuning .

Terminating Condition

An ellipse should not be too small to be printed. We regard (see Equation 1) as the terminating condition of inserting new ellipses. In other words, whenever is produced and shrunken, we check its horizontal axis and terminate if it is less than .

5 Extrusion to 3D

5.1 3D Hollowing

Extrusion of Ellipses

The hollowed polygon is lifted to 3D by extruding the ellipses orthogonal to the 2D plane in both directions. Denote as the index of , i.e., . For each ellipse in , we project it onto . If totally lies in within a distance of minimal wall thickness , we keep it in . Otherwise, we shrink it with a factor so as to the shrunk ellipse lies in within a distance of . We can also enlarge the ellipse if there is much space around it. Note that the enlargement of the ellipse should meet the support-free condition (Equation 1). This operation is successively applied for the other cross sectional polygons.

Hollowing Other Polygons

After we complete the extrusion of all ellipses for all cross sectional polygons, we check each polygon and choose the one with the largest available region which can insert more ellipses. Then we set it as input and add more support-free ellipses in it using our method, and then extrude the newly-added ellipses to its neighborhood polygons. The above process is iteratively performed until there is no any ellipse which can be added into the polygons.

Hollowed Volume

After we obtain the hollowed polygons, we connect the corresponding ellipses on successive polygons and thus generate a hollowed volume of . As each polygon is support-free, the obtained hollowed 3D volume is also support-free.

5.2 Functional Constraints

The printed objects are generally required to meet some functional constraints such as static balance and mechanical stiffness Wang et al. (2013); Wu et al. (2016a). However, integration of these constraints into the VD computation is computationally expensive. Thus we handle them as a postprocess after we generate the support-free hollowed volume.

Design Variables

We define a design variable , as a shrinking factor, for each ellipse inside . A value means that is totally filled with solid. The basic idea is that shrinking of ellipses can shift the center of gravity of the model and can improve its mechanical stiffness as more material is filled. Thus we can easily formulate the optimization according to a specific objective function. In particular, we discuss about the optimization with respect to static balance as an example. The optimization for other constraints can be similarly achieved.

Static Balance

An object is self-balanced when the vertical projection of its gravity center lies in the convex hull of its contact points with the ground. As shown in Figure LABEL:fig:balance, it is intuitive that shrinking ellipses on the left hand side of the gravity center will shift it leftwards, i.e., closer to the convex hull of its contact points. It is easy to formulate an optimization of minimizing the horizontal distance between the gravity center and the boundary of the convex hull. Our optimizer thus tries to reach a balance by shrinking some ellipses and thus shifting the gravity center into the convex hull.

balance0.360.36 Optimization of static balance (2D case). (a) The hollowed object cannot stand by itself (left). (b) It is optimized to a self-balanced object using our optimizer (right). The red dots denote the gravity center.

6 Experimental Results

Computational Platform

We have implemented our algorithm in C++ on a standard desktop PC with Intel(R) Core(Tm) i7-4790K CPU@4.0 GHz and 16 GB of RAM. Thanks to the efficient implementation of TOI-EinP, the VD generation of polygons and ellipses is fast and the ellipse hollowing is less than 30 seconds for all examples.

Printer Configuration and Parameters

We fabricate the objects using a commercial FDM 3D printer: The Ultimaker 2+ with tray size of . The printable layer thickness of the printer (printing precision) ranges from to and we use a value of . We test the plastic PLA material used in the 3D printing and set the maximally allowed overhang-angle as and the maximal length of printable horizontal hangover as .

Manufacture Setting

After we generate the hollowed models, we add supporting structures for the exterior part of the models but do not add any interior support. After the models are fabricated, we manually remove the exterior supports. All models have been successfully printed which reveals that the hollowed interior of the models is printed without any problem. We also validate this by printing and checking only half of the models which will be shown later.


Figure LABEL:fig:hollowing-volume shows an example of hollowed bunny model with different cross-sections. Figure LABEL:fig:cross-sections-hollowing shows a 3D hollowed volume of kitten model with a set of cross-sections, and then adding more in later passes. Figure LABEL:fig:balance-3D shows two hollowed hanging balls. The left one cannot stand by itself. After using our optimizer, it can be optimized to be well balanced in the right by filling material in the elliptic voids of the column. Figure LABEL:fig:gallery shows photos of the fabricated bunny model and kitten model which are hollowed by our method. The models are successfully fabricated without adding any extra support in the interior voids.

hollowing-volume1.2A hollowed bunny model. (a) the hollowed model; (b-d) different cross-sections.

cross-sections-hollowing0.450.450.450.45 A hollowed bunny model. (a) the hollowed model; (b-d) different cross-sections.

balance-3D0.480.48 (a) The hollowed hanging ball may fall down without balance optimization (left) (b) while the optimized one is standing stable.

gallery1.0 Photos of the fabricated bunny model and kitten model hollowed by our method.

Performance of the TOI-EinP Algorithm for Constructing VD

We conducted computational efficiency test using two models: the bunny and the hanging ball models. From the bunny model, we produced 121 planes resulting 219 polygons as some planes contain more than one polygon. The smallest and the largest polygons have 8 and 655 P-edges, respectively. For experimental purpose, we enforced to pack 100 ellipses into each polygon ignoring the mechanical and physicochemical constraints. Figure LABEL:fig:computation-time-profile(a) shows computation time vs. polygon size in terms of the input P-edges. The top-most black curve: the total time; The next red one: the time for incrementing all the in-disks of the ellipses into the VD structure; The next green one: that for finding the maximum clearance probe; The blue one: that for constructing the VD of the at-disks and on-disks. Note that the total time is weakly super-linear mainly due to the increment process of in-disks which is believed to be caused by the mapping mechanism of equivalent V-vertices between and . We used the map in the C++ template which is implemented by a binary search tree, taking time for each query for entities. With this model, we used non-uniform on-disk generation method using a bucket system for the acceleration.

From the hanging ball model, we produced 88 planes resulting 141 polygons: The smallest and the largest polygons have 25 and 232 P-edges, respectively. As this model consists of several large planar faces together with smaller ones, the polygons have several long P-edges together with short ones. This is common in many engineering products. Hence, we subdivided each P-edge into a set of P-edges of the length defined by the previously stated rule. Figure LABEL:fig:computation-time-profile(b) also shows computation time vs. polygon size in terms of the input P-edges. Note that the correlation is very weak compared to the bunny model as is expected because of the big variation of the P-edge lengths. Figure LABEL:fig:computation-time-profile(c) shows computation time vs. the number of subdivided P-edges: The curves are fairly well correlated in a slightly super-linear fashion. The big gab between the two clusters of data is due to many subdivided P-edges through very long input P-edges. For example, the left-most data in the right cluster in Figure LABEL:fig:computation-time-profile(c) corresponds to a polygon with 234 input P-edges which was subdivided into 513 shorter P-edges. Figure LABEL:fig:computation-time-profile(d) shows the number of subdivided P-edges vs. the number of input P-edges. The bunny model is expectedly a straight line whereas the hanging ball model shows bumpy curve which is similar to the curves in Figure LABEL:fig:computation-time-profile(b). The relationship between Figure LABEL:fig:computation-time-profile(b) and (c) can be explained by Figure LABEL:fig:computation-time-profile(d). Figure LABEL:fig:computation-time-profile(e) shows the packing ratio of the 100 ellipses in each polygon. The bunny model is expectedly smooth with decreasing pattern for bigger polygons as we incremented only 100 ellipses whereas the curve of the hanging ball model is bumpy.

computation-time-profile0.4680.4680.4680.4680.468 Computation time profile of the TOI-EinP algorithm. (a) Bunny model (Time vs. # input P-edges). (b) Hanging ball model (Time vs. # input P-edges). (c) Hanging ball model (Time vs. # subdivided P-edges). (d) # P-edges changes of both Bunny and Hanging ball models after subdivision. (e) Packing ratio of polygons in both Bunny and Hanging ball models.

Comparison to rhombic cell structure

The work of Wu et al. (2016a) adopts rhombic cell structure, which have discontinuity on boundaries, to generate support-free interior voids for 3D shapes. We compare our method with this method as shown in Figure LABEL:fig:Comparison-P-Model. We apply two methods on the same P model with similar hollowing ratios. Then we fix the bottom of the model and conduct an identical external load on it, respectively. From the stress map we can see that the result generated by Wu et al. (2016a) suffers the problem of stress concentration at the region marked in red, which generally happens in discontinuity. This does not happen in our method.

Comparison-P-Model0.85 Comparison with [Wu et al. 2016b]. Upper row: the hollowed P model using [Wu et al. 2016b] with a hollowing ratio of 24.3%; Lower row: the hollowed P model using our method with a hollowing ratio of 25.7%. The bottom of the model is fixed and one identical external load is conducted on the right, respectively, as shown by the arrow. The right figures show the color maps of stress. It is seen that the region marked in red in the upper-right figure suffers the problem of stress concentration, i.e., the stress there is very high.

7 Conclusions and Future Work

In this paper we propose a novel approach for generating support-free interior hollowing for general 3D shapes. The generated 3D shapes can be directly fabricated with FDM 3D printers without any usage of extra supports in interior voids. This is based on the observation of a family of support-free ellipses and is achieved by hollowing 2D shapes with these ellipses. Then the interior ellipses are extruded into volume for generating hollowed 3D shapes. We also develop a new, efficient and robust algorithm for the Voronoi diagram polygons and the first algorithm for the Voronoi diagram of ellipses within a polygon, both based on the topology-oriented incremental approach, which are quite useful for generating the ellipse packing in 2D shapes. With the sizes of ellipses as design variables, the optimization according to a specific objective function, e.g., static stability, can be easily formulated. Experimental results have shown the practicability and feasibility of our proposed approach.

Limitation and Future work

Our research opens many directions for future studies. First, the packing results can be further optimized by optimizing the positions and sizes of the ellipses for the purpose of increasing packing ratio. Second, it is expected to extend our approach for generating support-free ellipsoids for 3D shapes. This is feasible but needs more effort. Last but not the least, we are interested in studying general support-free shapes for additive manufacturing, which is a promising direction for geometric modeling and processing.


This work was supported by the National Research Foundation of Korea(NRF) grant funded by the Korea government(MSIP) (No. 2017R1A3B1023591).



  • Stava et al. (2012) Stava, O., Vanek, J., Benes, B., Carr, N., Měch, R.. Stress relief: Improving structural strength of 3D printable objects. ACM Trans Graph 2012;31(4):48:1–48:11.
  • Wang et al. (2013) Wang, W., Wang, T.Y., Yang, Z., Liu, L., Tong, X., Tong, W., et al. Cost-effective printing of 3D objects with skin-frame structures. ACM Trans Graph 2013;32(6):177:1–177:10.
  • Lu et al. (2014) Lu, L., Sharf, A., Zhao, H., Wei, Y., Fan, Q., Chen, X., et al. Build-to-last: Strength to weight 3D printed objects. ACM Trans Graph 2014;33(4):97:1–97:10.
  • Strano et al. (2013) Strano, G., Hao, L., Everson, R.M., Evans, K.E.. A new approach to the design and optimisation of support structures in additive manufacturing. International Journal of Advanced Manufacturing Technology 2013;66(9):1247–1254.
  • Dumas et al. (2014) Dumas, J., Hergel, J., Lefebvre, S.. Bridging the gap: Automated steady scaffoldings for 3D printing. ACM Transactions on Graphics 2014;33(4):1–10.
  • Vanek et al. (2014) Vanek, J., Galicia, J.A.G., Benes, B.. Clever support: Efficient support structure generation for digital fabrication. Comput Graph Forum 2014;33(5):117–125.
  • Langelaar (2016) Langelaar, M.. Topology optimization of 3d self-supporting structures for additive manufacturing. Additive Manufacturing 2016;12;Part A:60–70.
  • Wu et al. (2016a) Wu, J., Wang, C.C.L., Zhang, X., Westermann, R.. Self-supporting rhombic infill structures for additive manufacturing. Computer-Aided Design 2016a;80:32–42.
  • Reiner and Lefebvre (2016) Reiner, T., Lefebvre, S.. Interactive modeling of support-free shapes for fabrication. In: EUROGRAPHICS. 2016,.
  • Pilkey and Pilkey (2008) Pilkey, W.D., Pilkey, D.F.. Peterson’s stress concentration factors. 3 ed.; Wiley; 2008. ISBN 978-0-470-04824-5.
  • Karthikeyan (2016) Karthikeyan, K.. Why hatches and doorways in ships and airplanes are oval? 2016.
  • Zhang et al. (2015) Zhang, X., Xia, Y., Wang, J., Yang, Z., Tu, C., Wang, W.. Medial axis tree-an internal supporting structure for 3D printing. Computer Aided Geometric Design 2015;35-36(5):149–162.
  • Hu et al. (2015) Hu, K., Jin, S., Wang, C.C.. Support slimming for single material based additive manufacturing. Computer-Aided Design 2015;65:1–10.
  • Musialski et al. (2015) Musialski, P., Auzinger, T., Birsak, M., Wimmer, M., Kobbelt, L.. Reduced-order shape optimization using offset surfaces. ACM Transactions on Graphics 2015;34(4):102;1–9.
  • Prévost et al. (2013) Prévost, R., Whiting, E., Lefebvre, S., Sorkine-Hornung, O.. Make it stand: Balancing shapes for 3D fabrication. ACM Trans Graph 2013;32(4).
  • Christiansen et al. (2015) Christiansen, A.N., Schmidt, R., Bærentzen, J.A.. Automatic balancing of 3d models. Computer-Aided Design 2015;58:236–241.
  • Bächer et al. (2014) Bächer, M., Whiting, E., Bickel, B., Sorkine-Hornung, O..

    Spin-it: Optimizing moment of inertia for spinnable objects.

    ACM Trans Graph 2014;33(4):96:1–96:10.
  • Deaton and Grandhi (2014) Deaton, J.D., Grandhi, R.V.. A survey of structural and multidisciplinary continuum topology optimization: post 2000. Structural and Multidisciplinary Optimization 2014;49(1):1–38.
  • Wu et al. (2016b) Wu, J., Dick, C., Westermann, R.. A system for high-resolution topology optimization. IEEE Transactions on Visualization and Computer Graphics 2016b;22(3):1195–1208.
  • Hu et al. (2014) Hu, R., Li, H., Zhang, H., Cohen-Or, D.. Approximate pyramidal shape decomposition. ACM Trans Graph 2014;33(6):213:1–10.
  • Buchalter and Bradley (1992) Buchalter, B.J., Bradley, R.M.. Orientational order in amorphous packings of ellipses. Journal of Physics A: Mathematical and General 1992;26(3):L1219–L1224.
  • Schreck et al. (2010) Schreck, C.F., Xu, N., O’Hern, C.S.. A comparison of jamming behavior in systems composed of dimer- and ellipse-shaped particles. Soft Matter 2010;6(13):2960–2969.
  • Donev et al. (2004) Donev, A., Torquato, S., Stillinger, F.H., Connelly, R.. Jamming in hard sphere and disk packings. Journal of Applied Physics 2004;95(3):989–999.
  • Donev et al. (2007) Donev, A., Connelly, R., Stillinger, F.H., Torquato, S.. Underconstrained jammed packings of nonspherical hard particles: Ellipses and ellipsoids. Physical Review E 2007;75(5):051304.
  • Hit (2013) Optimized dropping and rolling (odr) method for packing of poly-disperse spheres. Applied Mathematical Modelling 2013;37(8):5715–5722.
  • Birgin et al. (2016) Birgin, E.G., Lobato, R.D., Martínez, J.M.. Packing ellipsoids by nonlinear optimization. Journal of Global Optimization 2016;65(4):709–743.
  • Delaney et al. (2005) Delaney, G., Weaire, D., Hutzler, S., Merphy, S.. Random packing of elliptical disks. Philisophical Magazine Letters 2005;85(2):89–96.
  • Lozano et al. (2016) Lozano, E., Roehl, D., Celes, W., Gattass, M.. An efficient algorithm to generate random sphere packs in arbitrary domains. Comput Math Appl 2016;71(8):1586–1601.
  • Jodrey and Tory (1985) Jodrey, W.S., Tory, E.M.. Computer simulation of close random packing of equal spheres. Phys Rev A 1985;32:2347–2351.
  • Lubachevsky and Stillinger (1990) Lubachevsky, B.D., Stillinger, F.H.. Geometric properties of random disk packings. Journal of Statistical Physics 1990;60(5):561–583.
  • Lubachevsky (1991) Lubachevsky, B.D.. How to simulate billiards and similar systems. Journal of Computational Physics 1991;94(2):255–283.
  • Specht (2015) Specht, E.. A precise algorithm to detect voids in polydisperse circle packings. In: Proc. R. Soc. A. 2015, p. 20150421.
  • Sugihara et al. (2004) Sugihara, K., Sawai, M., Sano, H., Kim, D.S., Kim, D..

    Disk packing for the estimation of the size of a wire bundle.

    Japan Journal of Industrial and Applied Mathematics 2004;21(3):259–278.
  • Lee et al. (2016) Lee, M., Sugihara, K., Kim, D.S.. Topology-oriented incremental algorithm for the robust construction of the voronoi diagrams of disks. ACM Transactions on Mathematical Software 2016;43(2):14:1–14:23.
  • Okabe et al. (1999) Okabe, A., Boots, B., Sugihara, K., Chiu, S.N.. Spatial Tessellations: Concepts and Applications of Voronoi Diagrams. 2nd ed.; Chichester: John Wiley & Sons; 1999.
  • Kim et al. (1995) Kim, D.S., Hwang, I.K., Park, B.J.. Representing the Voronoi diagram of a simple polygon using rational quadratic Bzier curves. Computer-Aided Design 1995;27(8):605–614.
  • Sugihara and Iri (1992) Sugihara, K., Iri, M.. Construction of the Voronoi diagram for “one million” generators in single-precision arithmetic. Proceedings of the IEEE 1992;80(9):1471–1484.
  • Held (2001) Held, M.. VRONI: An engineering approach to the reliable and efficient computation of Voronoi diagrams of points and line segments. Computational Geometry - Theory and Applications 2001;18(2):95–123.
  • Held and Huber (2009) Held, M., Huber, S.. Topology-oriented incremental computation of Voronoi diagrams of circular arcs and straight-line segments. Computer-Aided Design 2009;41(5):327–338.
  • CGAL (2016) CGAL, . CGAL Library Homepage. 2016.
  • Alliez et al. (2007) Alliez, P., Delage, C., Karavelas, M.I., Pion, S., Teillaud, M., Yvinec, M.. Delaunay tessellations and Voronoi diagrams in CGAL (draft). In: Rien van de Weijgaert Gert Vegter, J.R., Icke, V., editors. Tessellations in the Sciences. Virtues, Techniques and Applications of Geometric Tilings. 2007,.
  • Kim et al. (2001a) Kim, D.S., Kim, D., Sugihara, K.. Voronoi diagram of a circle set from Voronoi diagram of a point set: I. topology. Computer Aided Geometric Design 2001a;18:541–562.
  • Kim et al. (2001b) Kim, D.S., Kim, D., Sugihara, K.. Voronoi diagram of a circle set from Voronoi diagram of a point set: II. geometry. Computer Aided Geometric Design 2001b;18:563–585.
  • Sugihara (1993) Sugihara, K.. Approximation of generalized voronoi diagrams by ordinary voronoi diagrams. Graphical Models and Image Processing 1993;55(6):522–531.
  • Emiris et al. (2013) Emiris, I.Z., Tsigaridas, E.P., Tzoumas, G.M.. Exact voronoi diagram of smooth convex pseudo-circles: General predicates, and implementation for ellipses. Computer Aided Geometric Design 2013;30(8):760–777.
  • Emiris and Tzoumas (2005) Emiris, I.Z., Tzoumas, G.M.. Algebraic study of the apollonius circle of three ellipses. In: EuroCG. 2005, p. 147–150.
  • Farouki and Johnstone (1994) Farouki, R.T., Johnstone, J.K.. The bisector of a point and a plane parametric curve. Computer Aided Geometric Design 1994;11:117–151.
  • Alt et al. (2005) Alt, H., Cheong, O., Vigneron, A.. The voronoi diagram of curved objects. Discrete Comput Geom 2005;34:439–453.
  • Farouki and Johnstone (1992) Farouki, R.T., Johnstone, J.K.. Computing point/curve and curve/curve bisectors. In: IMA Conference on the Mathematics of Surfaces. 1992,.
  • Farouki and Ramamurthy (1998) Farouki, R.T., Ramamurthy, R.. Degenerate point/curve and curve/curve bisectors arising in medial axis computations for planar domains with curved boundaries. Computer Aided Geometric Design 1998;15(6):615–635.
  • Omirou and Demosthenous (2006) Omirou, S.L., Demosthenous, G.A..

    A new interpolation algorithm for tracing planar equidistant curves.

    Robotics and Computer-Integrated Manufacturing 2006;22:306–314.
  • Cao and Liu (2008) Cao, L., Liu, J.. Computation of medial axis and offset curves of curved boundaries in planar domain. Computer-Aided Design 2008;40(4):465–475.
  • Cao et al. (2009) Cao, L., Jia, Z., Liu, J.. Computation of medial axis and offset curves of curved boundaries in planar domains based on the cesaro s approach. Computer Aided Geometric Design 2009;26:444–454.
  • Barnhill et al. (1987) Barnhill, R.E., Farin, G., Jordan, M., Piper, B.R.. Surface/surface intersection. Computer Aided Geometric Design 1987;4:3–16.