A framework for adaptive width control of dense contour-parallel toolpaths in fused deposition modeling

3D printing techniques such as Fused Deposition Modeling (FDM) have enabled the fabrication of complex geometry quickly and cheaply. High stiffness parts are produced by filling the 2D polygons of consecutive layers with contour-parallel extrusion toolpaths. Uniform width toolpaths consisting of inward offsets from the outline polygons produce over- and underfill regions in the center of the shape, which are especially detrimental to the mechanical performance of thin parts. In order to fill shapes with arbitrary diameter densely the toolpaths require adaptive width. Existing approaches for generating toolpaths with adaptive width result in a large variation in widths, which for some hardware systems is difficult to realize accurately. In this paper we present a framework which supports multiple schemes to generate toolpaths with adaptive width, by employing a function to decide the number of beads and their widths. Furthermore, we propose a novel scheme which reduces extreme bead widths, while limiting the number of altered toolpaths. We statistically validate the effectiveness of our framework and this novel scheme on a data set of representative 3D models, and physically validate it by developing a technique, called back pressure compensation, for off-the-shelf FDM systems to effectively realize adaptive width.


Introduction
3D printing enables the fabrication of complex geometry under few design constraints compared to conventional fabrication techniques. Recent developments have seen a rapid growth in both the use and capabilities of desktop 3D printing systems. Fused Deposition Modeling (FDM) is one of the most common 3D printing techniques. It is widely used because of the versatility in the types of plastic which can be used and the relatively low running costs. FDM printers are used, for example, in showcasing scale models of buildings, casings for electronics, prototypes for blow molded parts, jigs and fixtures. Recent research adressed manufacturing complex volumetric structures such as microstructures [1,2,3] and topology optimized structures [4,5,6]. Many of these applications involve 3D models with detailed features within the order of magnitude of the nozzle size, which restrains the field of the process planning algorithms.
FDM printers extrude semi-continuous beads of molten plastic through a nozzle, which moves along a planned toolpath within each layer of a 3D object. A common strategy is to extrude along a number of parallel toolpaths which follow the shape of the contour of the layer and fill up the remaining area using parallel straight toolpaths. Contour-parallel toolpaths fit to the layer outlines more accurately, because the resolution of the positioning system is an order of magnitude smaller than the size of the hole in the nozzle. This paper is concerned with the generation of such contour-parallel toolpaths and addresses several issues which commonly occur in 3D models with narrow geometry.
The simple technique for generating the dense contour-parallel toolpaths of a layer consists of performing uniform inward offsets with the size of the nozzle from the outline shape. However, for geometrical features which are not an exact multiple of the nozzle size this method produces areas where an extrusion bead is placed twice: overfill areas; and areas which are not filled at all: underfill areas. See Fig. 1a. Overfills cause a buildup of pressure in the mechanical extrusion system, which can result in bulges or even a full print failure. Underfills, on the other hand, can cause a drastic decrease in the part stiffness or even for small features not to be printed at all. These problems are exacerbated for models with layer outlines with small features, because the over-and underfill areas are relatively large compared to the those features.
One promising direction to avoid over-and underfills is to employ toolpaths with adaptive width. Ding et al. developed a toolpath strategy for wire and arc additive manufacturing which produces a width variation typically lower than a factor of 3, but is far greater for some parts [7,8]. However, the range of bead widths manufacturable by FDM systems is limited. A nozzle of w = 0.4 mm will typically start to cause fluttered extrusion around lines narrower than 0.3 mm and lines will start to bulge upward if they are wider than the flat part of the nozzle, which is typically 1.0 mm.
The current state of the art of contour-parallel toolpath generation developed by Jin et al. employs a strategy which alters the widths of the centermost beads within a range of widths [0.25w, 1.8w] [9], which is similar to the strategy employed by the open source industry standard software package Ultimaker Cura [10]. See Fig. 1b. Still, controlling the extrusion width through movement speed changes or through volumetric flow control (e.g. linear advance) yields diminishing accuracy for deposition widths deviating more from the nozzle size, since process parameters such as nozzle temperature are optimized for beads with the nozzle size. Moreover, reducing the variation in width is beneficial for limiting the variation in mechanical prop- Illustration of different toolpaths for a shape showcasing a range of shape radii (black). These results can be read as a graph with feature size on the horizontal axis and its corresponding beading along the vertical axis. (a) Toolpaths using uniform offsetting results in large overfill (orange) and underfill (azure). (b) Toolpaths with adaptive width [9] where beads that are wider or narrower than the nozzle size are indicated in red and blue, respectively. (c) Our approach minimizes over-and underfill with less extreme widths. erties of the resulting product, meaning it conforms better to a simulation which employs a homogeneity assumption. We therefore reduce the bead width range by distributing the workload from the centermost bead over neighboring beads.
Our contributions are as follows: • A geometric framework allowing various adaptive bead width control schemes used to generate contour-parallel toolpaths which minimize under-and overfill.
• A specific beading scheme, which reduces the variation in the extrusion widths to within [0.75w, 1.5w].
• A back pressure compensation approach to accurately realize adaptive bead width on Bowden style hardware systems.

Related Work
Toolpath generation consists in generating a path in the a planar contour, representing the intersection of a plane and a 3D solid object. The nozzle is then instructed to move along the path while extruding material. Sites along the toolpaths are assigned several properties such as movement speed, but for this paper we will focus on the assigned width of the extruded bead. Toolpath generation is an integral part of process planning for 3D printing. For an overview of the processing pipeline, we refer to the survey by Livesu et al. [11]. For reducing printing time and material cost, sparse infill structures such as triangular and hexagonal patterns have been used to approximate the interior of 3D shapes. In this paper, we focus on the generation of toolpaths for dense regions, such as the boundary shell of 3D shapes. This is sometimes called 'dense infill' [11].
The toolpath has a direct influence on the printing time, material cost, and mechanical properties of the printed object [12,13]. FDM calls for toolpaths with several desirable properties. First, the extrusion path should be as continuous as possible. A discontinuous path requires to stop and restart material extrusion. For certain materials, such as TPU, this could lead to printing defects or even print failure [14]. Second, the toolpath is preferred to be smooth. Sharp turns require to reduce the movement speed of the nozzle, and so this prolongs the printing process. Third, the extruded path should cover the region of the contour without underfilling. Such underfill negatively influences the mechanical performance of the parts. Fourth, the extrusion paths should not overlap with one another. Such overfill causes a pressure build up in the mechanical system, which leads to overextrusion in later paths and in extreme cases cause print failure [14]. An analysis of under-and overfill from a vertical cross-section was presented in [15]. Our method is primarily concerned with minimizing under-and overfill within horizontal cross-sections.
Two basic strategies for dense toolpath generation are the direction-parallel strategy and the contour-parallel strategy. Direction-parallel (or zig-zag) toolpaths fill an arbitrarily shaped contour with a set of parallel, equally spaced line segments. These parallel segments are linked together at one of their extremities, to avoid discontinuous extrusion. Contour-parallel toolpaths typically consists of a set of equally spaced offsets from the slice boundary outline polygons. Steuben et al. presented a method for generating sparse infill toolpaths based on the isocontours of surface plots of some variable generated on each 2D contour [16]. In order to increase the continuity of contour-parallel toolpaths, a strategy to connect dense toolpaths into spirals was introduced by Zhao et al. [17] and later extended to also connect a mixture of dense and sparse toolpaths together [14]. Jin et al. discusses several approaches for connecting direction-parallel and contour-parallel toolpaths into continuous paths [18]. Spiral toolpaths have also been applied to (CNC) machining [19,20]. One of the problems with contour-parallel toolpaths is that it tends to leave gaps between the toolpaths (see Fig. 1a). This is due to the fact that the diameter of the part is not exact multiple of the (constant) deposition width in those regions. To avoid problems with such gaps, hybrid approaches that combine direction and contour-parallel are often used [21,22]. Close to the slice boundary, there are several contour-parallel curves, while the interior is filled using a zig-zag pattern. For complex shapes, the entire cross-section could be decomposed into a set of patches, and for each of them the basic strategies can be applied [23,18]. Alternative toolpath patterns, seen also in CNC machining, include space-filling curves [24,25,26].
Reducing under-and over-filling can be accurately achieved by making use of adaptive deposition width. Adaptive width can be used to locally match the nonuniform space between adjacent paths, and thus to ensure a better filling of the area. Kao and Prinz propose smooth adaptive toolpaths based on the medial axis skeleton of the boundary contour [27]. Their approach handles simple geometry where there are no branches in the medial axis. An extension was proposed by Ding et al. to handle complex shapes [7]. However, this extension inherits a problem in the original method: from any point in the skeleton to the boundary, the number of toolpaths is constant. Depending on the size of small and large features in the layer outlines, this strategy can require a range of bead widths beyond the capabilities of the manufacturing system. Jin et al. proposed a strategy of adding toolpaths with varying width along the center edges of the skeleton, while leaving other paths unchanged [9]. The resulting beads have widths within the wide range of [0.25w, 1.8w] (see Fig. 1b). In this paper we propose a novel scheme to distribute the width alterations throughout a region around the center, and thus limit the occurrence of extreme variation in width (see Fig. 1c).
Under-and over-filling issues have a large proportional impact for thin geometric features. Jin et al. proposed a sparse wavy path pattern for thin-walled parts [28]. Besides under-and overfilling, there are a few other robustness issues in toolpath generation for thin geometric features. Moesen et al. proposed a method to reliably manufacture thin geometric features using laser-based additive manufacturing techniques [29]. ¿¿ achange!! Behandish et al. presented a method to characterize local-topological discrepancies due to material under-and over-deposition, and used this information to modify the as-manufactured outcomes [30].
For ease of reference we have included a legend showing the terms employed in this manuscript in Fig. 4. These terms will be   further explained as they first appear throughout this paper.

Overview
Given arbitrarily shaped polygons which represent the 2D outline of a layer of a 3D model, our method generates extrusion toolpaths with varying width, i.e. a set of polylines where each site consists of a location and an extrusion width; in between the sites we linearly interpolate the position and extrusion width.
Our method starts with computing a graph which represents the topology of the input polygon: its skeleton. Our skeleton is based on the medial axis transform (MAT), a strategy that has been commonly used for generating contour-parallel toolpaths [31]. We visualize the skeleton as the union of cones (UoC) (Section 3.2), by raising each point in the domain to a height that equals the shortest distance from the point to the polygon boundary (Fig. 2c). Contour-parallel toolpaths with uniform width can be interpreted as the intersection of the union of cones with a set of horizontal planes at equally spaced heights (Figs. 2a to 2e).
As depicted in Fig. 2f, our method first identifies edges and nodes of the skeleton in the center of the polygon, which correspond to ridges and peaks in the UoC (Section 3.3). The heights b at these elements are then quantized to an integer number of beadsb. To ensure a smooth toolpath between regions with quantized integer heights that differ, we add new nodes in the skeleton with quantized heights and interpolate the heightsb in between (Section 3.4). The union of cones corresponding to the smoothed skeleton is then sliced at regular intervals to obtain toolpaths with varying spacing, which translates into varying width (Section 3.6). The video in the supplementary material provides an example animation of this approach.
This section explains how we generate toolpaths using our framework with uniform bead widths and evenly distributed locations between the center of the polygon and the outline. In Section 4 we describe how to apply the framework to different beading schemes and we show several such beading schemes.

Union of cones
The union of cones (UoC) is derived from a common skeletonization of the polygonal outline shape: the medial axis. By assigning each node in the skeleton a height equal to its shortest distance to the outline we obtain the shape of the UoC. Starting from the medial axis we further decompose the shape into simple fragments, so that the domain contains only quads and triangles. This decomposition constitutes an approximation of the UoC.
Medial axis transform. The medial axis is a representation commonly used to analyze a shape. It is defined as the set of positions where the inscribed circle meets the boundary in at least two locations [32,33]. The resulting skeleton consists of straight edges and parabolic edges. An example is illustrated in Fig. 3a. We call the set of points on the outline polygon P closest to a skeletal point v its support: The shortest distance for a point on the skeleton is called its feature radius, R(v). The medial axis along with the feature radius values along the skeleton form a complete shape descriptor, known as the medial axis transform (MAT). By vertically raising the center of an inscribed circle to a height that equals the center's feature radius, a cone is formed. The union of all such cones forms a 3D solid volume. The medial axis can thus also be interpreted as ribs of the surface of the union of cones [32].
Skeletal trapezoidation. Starting from the medial axis we decompose the input polygon into a set of quads and triangles, so that we can perform the slicing stage on simple shapes. We employ a shape decomposition similar to the one proposed by Ding et al. [7]. The basic idea is to add edges connecting each node v on the medial axis to each of its support points sup(v). The resulting skeleton decomposes the outline shape into trapezoids and triangles. Considering the fact that the concept of trapezoidation   conventionally allows for the degenerate case where a trapezoid resolves into a triangle [34,35], we call this shape decomposition the Skeletal Trapezoidation (ST). The edges generated by the MAT are classified into three types: 1. line-line edge -straight edge generated from two line segments in the outline polygons, 2. vertex-line edge -parabolic edge resulting from an outline vertex and a line segment in the outline, and 3. vertex-vertex edge -straight edge resulting from two outline vertices. The vertex-line and vertex-vertex edges are discretized into pieces with a length up to 0.2 mm, which gives an approximation error of only ±0.01 mm. This allows to approximate the feature radius between two discretized nodes v 0 and v 1 by linear interpolation. Again we connect the newly inserted nodes to their support, which results in vertex-line regions and vertex-vertex regions such as depicted in Fig. 4.
Approximation of union of cones. The skeletal trapezoidation (ST) provides a means to visualize the union of cones (UoC) approximated by a 3D surface mesh composed of quadrilateral and triangular patches. We assign each node in ST a (real number) height value measured in terms of beads, referred to as the bead count b. We define the bead count as the number of beads to fit along the diameter of the inscribed circle centered at node v, i.e.
where w * is the nozzle size. We divide the diameter rather than the radius as this allows to deal with an odd number of beads while using integer logic. Note that although the overview of the method was described geometrically in terms of the UoC, the actual toolpath generation relies on the two-dimensional ST; the use of the bead count as a height value is only a visualization aid. The significance measure can be simplified using Implementation. The medial axis of a polygonal shape is a subset of the Voronoi Diagram generated from the line segments and vertices of the shape [33]. The edges of the Voronoi diagram that fall outside of the outline shape are irrelevant for our purpose and are thus discarded. Note that besides the full medial axis, the Voronoi diagram also contains edges connecting to concave vertices in the outline shape (see Fig. 3b). These extra edges are a subset of the edges connecting a node to its support, so we keep them in. From the Voronoi diagram we add nodes to discretize parabolic edges and edges formed by two concave outline vertices, and then connect all nodes to their supports, forming a skeletal trapezoidation. We then assign each node the bead count values using Eq. (2). We compute the Voronoi diagram using the Boost C++ libraries [36], which implements the algorithm proposed by Fortune [37]. A half-edge data-structure is used to represent the Voronoi diagram (Fig. 3d).

Center classification
In order to prevent over-and underfill from occurring in central regions, parts of the ST are marked as being central. Our framework will decide on a beading at all the marked nodes in 'the center' and apply the beading outward to the unmarked nodes (Section 3.5).
A node in the ST is marked as as central if its feature radius is larger than that of all its neighboring nodes, i.e. a local maxima. An edge and its two nodes are also marked as being central if it is significant according to a significance measure.
Significance measure. We make use of the bisector angle as an indicator of significance which is commonly used in shape analysis. The bisector angle α is the interior angle ∠p 0 lp 1 ≤ 180 • , between any location l on an edge of the ST and its two supporting points {p 0 , p 1 } = sup(l) [38]. An edge is significant if the bisector angle on any location on the edge exceeds a prescribed α max . As illustrated in Fig. 5a, for a polygon with a pointy wedge area of an angle β, we have α = 180 • − β. This corresponds to overfill areas and underfill areas the size of 1 /4(w * ) 2 (tan(α/2) − α/2) when filled using the simple technique of uniform bead width w * . A too large α max may leave a lot of under-/overfill, while a too small value may introduce toolpaths to fill in negligibly small underfills. We therefore set α max = 135 • . Although significance measures are commonly used as a heuristic for finding the parts of a skeleton which are in some sense relevant [38,39], we use the bisector angle as an exact indicator of the amount of overfill and underfill in the uniform toolpaths of constant width.
To avoid evaluating the bisector angle at any location on all edges, we devise an efficient measure which operates only on the two nodes of an edge. Because all locations along a line-line edge have the same bisector angle we can evaluate whether the edge is significant by checking whether (see Fig. 5b). This ratio has a clear geometrical interpretation as the slope of the ridge in the UoC surface. For vertex-line edges and vertex-vertex edges only a portion of the edge is significant. We therefore introduce nodes at the boundaries of the significant portion during the discretization of such edges (see Appendix A).
The significance of all edges can then accurately be evaluated using Eq. (3).
Marking filtering. After initializing the marking at all edges and nodes, we filter out high frequency changes in the marking in order to ensure that the generated toolpath is smooth. The filtering is performed by additionally marking some unmarked elements, rather than the opposite since unmarking central regions reintroduces large over-and underfill areas. From each marked node v 0 with an upward unmarked edge attached we walk along the upward edges; if the total length traversed until we reach another marked node v 1 is shorter than some filter distance d unmarked max , we mark all edges encountered as being central. We use d unmarked max = w * in order to filter out high frequency oscillations in the order of magnitude of the nozzle size, while keeping close to the significance measure.

Central height adjustment
After the central regions have been identified, their heights are quantized. First, the initial bead countb is quantized into an integer bead countb at the marked nodes using a quantization operator q, then the locations along the edges where q makes a jump from one bead count n to another n + 1 are identified and then ramps are introduced to smoothly transition from n to n + 1 using fractional bead countsb along the smooth transition.
Quantization. We define a quantization operator q to map a feature diameter (d = 2R(v)) to a bead count: q : R → N. Because our quantization scheme should round to the nearest integer multiple of the nozzle size, we have q(d) = d/w * + 1 /2 . Alternative quantization schemes are discussed in Section 4. By applying q to the heights of central nodes we quantize the bead count: Transition anchors. For a marked edge which connects nodes v 0 and v 1 withb v 0 ≤ n <b v 1 , we determine the transition anchor locations at which the bead count transitions from n to n + 1. To this end, we introduce the function which gives the feature diameter d at which the bead count q transitions from n to n + 1. The location of the anchor v x is then computed by inversely interpolating An illustration of the anchors is shown in Fig. 6b. We perform a filtering step to prevent frequently changing the bead count back and forth within a short distance. For two consecutive anchors which transition to opposite directions, if the distance between them is smaller than some limit d transition max , the bead counts at all nodes in-between are set to the surrounding bead counts, and consequently these anchors are removed (See Fig. 6c). A value of d transition max = 1 mm seems to produce satisfactory results. This means that for some small regions we generate toolpaths with bead widths outside the typical range.
Smooth transitions. A sharp transition from n to n + 1 beads at an anchor location creates sharp turns in the toolpath (see Fig. 7 top). We introduce a transition length t(n) to ensure a smooth transition (see Fig. 7). The length of the transition is set to t(n) = w * and it is centered at the anchor, i.e. the distance from the lower end v 0 to the anchor position v x is set to where ∆ is the total distance along the edges between two nodes. The transition length t(n) ensures that the center beads don't overlap with the innermost transitioning beads, while keeping the amount of underfill low and the toolpath smooth. The transition anchor position t 0 (n) ensures that the transitions never overlap with each other or with locations where all beads have the preferred width w * .
We discard any transition anchor which is too close to the end of a chain of marked edges for the smoothed transition to fully fit within the marked region. In order to make the transition ramps robust against small perturbations in the outline shape which cause extra (support) edges in the skeleton, we modify the nodes v x which are between the two ends v 0 and v 1 of the transition by (re-)assigning them a fractional bead countb which is linearly interpolated between the two ends of the transition (see Fig. 6e Note that although the ST is not stable against noise in the boundary shape, the distance field itself is, so by designing our algorithms such that they are stable against changes in the topology of the skeleton our method is stable against small perturbations in the outline. Finally we update the ST by adding support edges at the transition ends. As shown in Fig. 6f, the marked regions in the UoC mesh have become horizontal at integer multiples of 1 /2w * for long stretches with ramps in between.

Beading
Now that we know how to determine the bead counts in the marked central regions the question is how the unmarked regions are handled. Determining bead count values for the unmarked nodes and interpolating linearly along the unmarked edges would mean that toolpath sites would be distributed evenly along each unmarked bone; while that would suffice for the evenly distributed beading scheme, it wouldn't allow for more sophisticated, non-linear schemes. Instead we determine the radial distance to the boundary at which each bead should occur from the boundary to the center. Each central node is associated with a sequence of radial distances L which control the locations of the beads, starting from the outer bead and ending in the center. Together with a sequence of bead widths W, these form what we call a beading B. For our distributed beading scheme we compute the beading for a central node v with n = b v beads and a diameter r = R(v) as: B(n, r) = (W(n, r), L(n, r)) = w 0 . . . w n/2 −1 , l 0 . . . l n/2 −1 w i = r/n for all i ∈ N : i < n/2 where w i and l i are the width and location of the ith bead, respectively, counting from the outline inward. Example beadings for  an odd and even bead count with arbitrary widths are visualized in Fig. 8a.
Beading interpolation. The beading is defined in terms of an integer number of beads, while we have assigned a fractional bead count to nodes within a transition region. In order to generate a beading for a node v with n < b * v < n + 1 we linearly interpolate the bead widths and locations between a beading B 1 based on n and a beading B 2 based on n + 1 (see Fig. 8). Such interpolation is also used to deal with beading conflicts (see Fig. 9). There we also apply beading interpolation from a marked node v m upward along unmarked bones, and interpolate between v m and the beading at the top of the slope over some distance t beading from the lower marked node, which we set to t beading = w * , so that the transition is not too swift.
Beading propagation. The beading information is then broadcast throughout the ST from central regions outward, so that each unmarked node v knows the beading of the marked node on top of the ramp on which v is placed. We first broadcast the beading information upward from all marked nodes, so that we can then deal with beading conflicts in a downward phase. The downward phase makes sure that all nodes have a beading associated with it, so that the slicing algorithm can efficiently slice the edges leading up to a marked or unmarked node.

Toolpath extraction
Once each node has been assigned a beading, we proceed to generate the toolpath sites along the edges of the ST. A site S consists of a location v a width w and an index i, which are computed for an edge v 0 v 1 from the beading B of the upper node v 1 : Fig. 10. We store all sites of an edge in a mapping from edge to a list of sites.
We then generate extrusion segments for each trapezoid by connecting together the sites of the same index. See Fig. 11. If the amount of sites on both sides of the trapezoid is not the same then this trapezoid is in a transition and we leave one inner site unconnected.  The beading conflict is resolved by gradually interpolating between the two beadings. The ramp to the upper ridge doesn't line up with the lower ridge, which means that the toolpaths (dashed) resulting from the beading propagated from above doesn't align with the beading from the thin outline feature (highlighted in red).
(c) Separate beadings at either top node Because the bead count is defined in terms of the feature diameter rather than the radius, only some of the bead count valuesb in a central region coincide with a slicing height. When the bead countb is even, the ridge is sliced as normal; the intersection between a slicing plane and the mesh surface results in a polyline on both sides of the ridge, which are connected together into a polygonal toolpath. When the bead countb is odd, the ridge will coincide exactly with a slicing height, which results in a single polyline toolpath being generated along the middle of the feature. In that case we should prevent the algorithm from generating the center extrusion segment twice from the trapezoids on either side of that segment. We therefore use some arbitrary condition to decide which one of the two to include based on the ordering of the coordinates of v 0 and v 1 : All trapezoids in the ST are assigned to separate domains, corresponding to which boundary polygon they are connected to (see Fig. 11a) [7]. By traversing the trapezoids per domain in order we can efficiently connect all segments into polylines. See Fig. 11. In a final step we connect the ends of polylines together, so that the final toolpaths contain both polygons and polylines.
Around the transition locations and around nodes with odd bead count and more than two marked edges attached there will be intersections in the toolpaths. Such intersections cause overfill because the nozzle passes the location multiple times. We deal with this special case by forcing a new polyline when traversing the trapezoids, and in the final polyline connection step we greed- (b) Extrusion segment chaining Fig. 11. Generating toolpaths on a part of the test outline shape by chaining together extrusion segments along each polygon domain. Each edge is assigned toolpath sites (yellow) which are connected together as shown in the singled out trapezoid. By following the trapezoids along the domain (cyan) of a single outline polygon, the extrusion segments can efficiently be connected into existing polyline toolpaths (light and dark gray). ily connect the first two polylines ending in the same location and retreat all other polylines ending in that same location in order to prevent the overfill. In order to retreat a polyline which ends in a site S , we remove part of the polyline paths up to the intersection by a distance of w S d intersection max . We set d intersection max = 75 % in order to slightly favor overfilling over underfilling. This ratio effectively deals with the balance between overfill and underfill generated at that location after the retreat has been applied. See Fig. 12.

Beading schemes
A critical component in toolpath generation is how to distribute the beads over the feature radius. While the framework presented in the previous section takes evenly distributed beads as an example, it allows to apply different beading schemes to configure the bead distribution to cater for specific requirements from the application, 3D printer or material. Definition 4.1. A beading scheme is defined by the quantization operator q and the beading operator B: {q(d), B(n, r)}. The beading function B(n, r) consists of (W(n, r), L(n, r)), which provides sequences of n bead widths and of n distances from the outline to fill up a radial distance r.
For the smoothness and continuity of toolpaths we require that W n is monotonic and continuous at each bead index n for constant bead count c: 0 ≤ ∂W(c,r) n ∂r ≤ 1. We further ensure that beads don't overlap, that beads are extruded from the center of where they end up and that odd bead counts produce a single polyline toolpath exactly in the center by determining the bead locations from the widths: We introduce several beading schemes which determine the bead count and their widths in various ways. We can emulate a variety of toolpath generation methods from related literature by defining new beading schemes. We also introduce new beading schemes which produce toolpaths with less extreme widths compared to techniques from existing literature.
Uniform beading scheme. We can define a beading scheme which emulates the uniform width offsetting technique by disabling the marking of edges, so that we never employ transitioning. We can simply set α max = 180 • and supply a simple beading scheme given by Table 1a.
Outer bead. We can emulate the method from Moesen et al. by carefully choosing how the beading scheme functions deal with the outermost bead. Also we turn off the reduction of toolpaths near 3-way intersections d intersection max = 0 %, so that the polygonal toolpaths emulate the remaining area to be filled by another path planning technique similar to their technique. We don't need transitioning, so we also set t(n) = 0. See Table 1b.
Constant bead count. We can emulate the method from Ding et al. by dividing the feature radius over the widths of a constant number of beads. Additionally in order to emulate their definition of "branches" we mark all ST edges (α max = 0 • ) and we unmark the outer edges connected to the outline shape in a separate algorithm. Note that this deviation from the proposed framework violates the robustness against small perturbations in the outline polygon, since this marking depends on the topology of the graph of the ST. See Table 1c.
Centered. We can emulate the method from Jin et al. by transcribing how they deviate from the uniform width toolpaths. We therefore base the beading scheme on the bead count q − (d) defined by the uniform beading scheme. Jin et al. replace two beads from the uniform toolpaths by a single one when the distance between the center of those beads falls short of d min = 0.8w * , which gives us w max = d min + 2 · 1 2 w * = 1.8w * . Conversely, they place an extra bead when the distance exceeds d max = 1.25w * [9], which gives us w min = d max − 2 · 1 2 w * = 0.25w * [9, p. 72]. We emulate the rounded polygonal path rerouting they define by supplying a transition length t(n) = 1 2 w * which results in a discretized version of their rounded polygon segment. See Table 1e.
Evenly distributed. By taking the advantages of the above two schemes we can define a beading scheme which constitutes a novel toolpathing technique. We can evenly divide the local feature diameter over the widths of all beads, but choose a local bead count better matching the local feature size. We determine the local bead count by dividing the diameter by the preferred bead width and rounding to the nearest integer. This reduces the demands on the system and deviation from mechanical properties caused by beads with extreme deviations from the preferred width. See Table 1d.
Inward distributed scheme. The evenly distributed scheme can be conceptualized as calculating the total discrepancy E between the actual feature diameter d and the total preferred width nw * , dividing the total discrepancy by the number of beads and setting the width of each bead to w * + E/n. However, depending on the application we might want a different distribution of widths. We therefore supply a beading scheme which supports an arbitrary distribution of the discrepancy. The distribution is determined by some weighing function ω(n, r), which defines the portion of the discrepancy to distribute to each bead. See Table 1f. For example, we can choose an ω which distributes the discrepancy over the innermost N beads, and distribute most of it to the inner beads. See Fig. 13. That way we limit the region of impact of the distributed scheme to a central region and have the preferred bead width w * in regions farther away. This limits the impact of transitioning regions so that transitions keep the toolpaths smooth farther away from the central regions.
Widening meta-scheme. Complementary to any of these schemes we can enforce a minimum feature size and minimum bead width in our framework. Regions where the model is narrower than some r min can be printed with a bead width w min larger Shell meta-scheme. The industry standard of FDM is to generate only a limited contour-parallel perimeters and to fill the remainder using a direction-parallel strategy. We therefore provide a meta-scheme to generate adaptive bead width toolpaths only in narrow regions and generate the limited number of perimeters M using the preferred width in regions which are wide enough. We also take care not to leave gaps which are too small to be filled using the direction-parallel strategy: These meta-schemes introduce non-linearities in the quantization function. Because the beading is only evaluated at nodes in the skeleton, we need to make sure that there are nodes at the locations along the skeleton where the non-linearities happen. We therefore insert extra nodes along with their ribs at locations v with a radial distance R(v) = r min for widening and at R(v) ∈ Mw * , q −1 (M), q −1 (M) + 1 /2w * for the transition from narrow shell to unconstrained shell. Combining all metaschemes functionality we can generate results such as depicted in Fig. 14c.

Fabrication
In order to accurately manufacture adaptive width toolpaths using an off-the-shelf 3D printing system, we need a model which relates the required width to process parameters such as movement speed and filament extrusion speed. A different approach might be appropriate depending on whether the filament feeder is mounted directly on the print head (a.k.a. direct drive) or the filament fed from the back of the printer to the print head via a Bowden tube. Because Bowden style 3D printing systems have the filament feeder relatively far away from the nozzle, changing the internal pressure in the system requires a large amount of filament movement, which requires a prohibitive amount of time.

Back pressure compensation
Because changing the internal pressure is difficult in our setup, we keep the internal pressure constant, and vary the movement speed instead. One approach would be to keep the filament inflow f (in mm 3 /s) constant by varying movement speed accordingly [40]. However, that doesn't result in the intended filament outflow variation -see Fig. 15a. We conjecture that the filament outflow is related to the total pressure in the system, which depends not only on the amount of filament in between the feeder wheel and the nozzle (which we keep constant), but also depends on the back pressure that the previous layer exerts on the filament protruding from the nozzle. The amount of back pressure is most likely monotonically related to the requested line width. We compensate for the back pressure using a simple linear model: where v(w) is the movement speed as a function of requested bead width w, f (w) is the filament outflow, f 0 is a constant reference flow, w 0 is a constant reference bead width and k is the amount of back pressure compensation. Our back pressure compensation method effectively changes the speed to realize adaptive width, but this approach is limited, since the movement speed is constrained by acceleration considerations near bends in the toolpath [41]. Moreover, as the layer height is decreased the back pressure becomes larger compared to the internal pressure, which might cause the back pressure compensation method to demand prohibitively slow movement speeds. Furthermore, the shape and filling of the previous layer might influence the amount of back pressure. Accurate flow control can be further enhanced by using a direct drive hardware system and by employing pressure advance algorithms which dynamically change the internal pressure [42]. Conversely such a setup might benefit from some form of back pressure compensation as well.

Print results
Using increments of 0.1 we established that using a factor of k = 1.1 yields satisfactory bead width variation for our setup where we use f 0 = v 0 w 0 h with v 0 = 30 mm/s, w 0 = 0.4 mm and h = 0.1 mm. See Fig. 15b. The fact that the printed lines are wider than intended is compensated for using a flow reduction to 90 %. Test prints were performed on an unmodified Ultimaker S5 system, with a standard 0.4 mm nozzle and PLA filament. The printing order is determined greedily by choosing the closest point of a polygonal extrusion path, or the closest of either end point in case of an open polyline extrusion path. Because the machine instructions file format G-code doesn't natively support adaptive width beads, we discretize adaptive width extrusions into 0.2 mm long segments of the average width. The print results can be viewed in Fig. 16.
In Fig. 16b the underfill problem of the naive uniform offset approach is most prevalent for the Ultimaker word mark, which negatively impacts the visual quality and the stiffness of the part. Moreover, in the case of the spatially graded honeycomb, there are several fully disconnected hexagons, which means the object falls apart when picked up. The honeycomb print is also missing all parts which are slightly more thin than the preferred bead width w * . Figure 16c still shows some underfill, but considerably less than the uniform approach. These prints also exhibit dark regions where the translucency of the layer is less because the bead is higher. This can be explained by inaccuracies in the back pressure compensation method, which arise for bead widths which deviate from the preferred width by a large amount. Figure 16d diminishes the underfill nearly completely and the visual quality of these prints is more homogenous than those of the other methods. Moreover, the absence of dark regions signifies that our proposed method is more robust against inaccuracies in the deposition system. However, both the centered and inward distributed approach introduce transitions to a different bead count in the word 'Delft', which reduces the dimensional accuracy on the outline around those locations.

Results and discussion
We evaluate the proposed framework and the various beading schemes on a set of different types of 3D models, ranging over various applications and various types of geometry. The data set is described in Appendix B. We sliced all models in the data set and selected 300 random slices for analysis. Toolpaths of these 300 outline shapes are generated using the uniform technique as implemented by Clipper [43] -a state-of-the-art polygon offset library, and by our framework using four beading schemes, i.e. the constant bead count scheme with a bead count of C = 4, the centered, the evenly distributed, and the inward distributed beading scheme using N = 2, all with a preferred bead width of w * = 0.5 mm and using the widening meta-scheme to enforce a minimum printed feature size of w min = 2r min = 0.3 mm.
The tests were performed on a desktop PC equipped with an Intel Core i7-7500U CPU @ 2.70 GHz (a single core is used) and 16.3 GB memory. We report on the total statistics summed over the whole data set, because averaging would be biased.

Accuracy
We first evaluate the accuracy of different beading schemes in terms of the relative amount of the overfill and underfill. We   construct the over-and underfill area by comparing the shapes covered by each extrusion move with each other and with the total shape of the boundary polygons. (For implementation details see Appendix C.) This results in polygonal shapes such as visualized in the top half of Fig. 17: there are orange shapes where the beads overlap and azure shapes in the voids in between the beads. We compare the total area in mm 2 of these overfill and underfill shapes to the total area of the boundary for each sample in the data set and report the average percentages in Fig. 18a. The inward distributed scheme has a calculated overfill of 0.30 % and an underfill of 0.24 %. This is lower compared to the uniform scheme, which results in 1.63 % overfill and 1.62 % underfill in the data set.

Uniformity
We visualize the bead widths resulting from the different schemes in the bottom of Fig. 17. We binned the toolpaths into width bins at 0.01 mm increments and determine the total toolpath length pertaining to each bin. From these statistics we calculate the mean and standard deviation and report them in Fig. 18d. We found that the mean width of the inward and evenly distributed schemes is close to the preferred bead width of 0.5 mm, while their standard deviation is lower than for the centered and constant bead count scheme. These results show  that, while causing less overfill and underfill, inwards distributed and evenly distributed schemes deviate less or less often from the preferred bead width compared to the other schemes.

Print time
The total time it takes to print a part is influenced not only by our back pressure compensation scheme, but also by the geometry of the toolpaths. In order to separate these effects we report on the total print time when using back pressure compensation and when using a constant (maximum) movement speed in Fig. 18b. We estimate print times using a simulation of the Marlin firmware using the default movement settings of the setup described in Section 5.2. While the idealized print time is predominantly determined by the total toolpath length, the print time using back pressure compensation is predominantly determined by the occurrence of wide beads, because they have a reduced the flow in mm 3 /s. Because of acceleration constraints imposed by the hardware the maximum movement speed is not reached near sharp corners. We therefore also report on the angles of the bends in the toolpaths in Fig. 18e. Furthermore, the print time is negatively affected by discontinuities in the extrusion process. Between extrusions the printer needs to stop extrusion, travel to the next extrusion path and restart the extrusion process, which may introduce defects and incurs extra print time. For closed polygonal toolpaths we can start anywhere within the path, so we can optimize the starting location so as to minimize the travel time. We therefore report both on the open and closed path count in Fig. 18c. Fig. 18f plots the computation time against the vertex count of the layer for the full data set, comparing the uniform technique implemented using Clipper [43] to our framework with the inward distributed scheme. For polygonal shapes with as many as 10 4 vertices, the computation for both approaches is less than 1 second, with our method being approximately five times that of the uniform technique. These results could be improved upon by utilizing the locality inherent in our algorithms for parallelization on the GPU.

Computational performance
The computational complexity is limited by the generation of the Voronoi Diagram, which is O(n log n), where n is the number of vertices in the input shape. The other steps in our framework have a complexity of O(m), where m is the number of elements in the ST. Therefore, the total running time of our algorithm is O(n log n). Results in Fig. 18f confirm that both our framework and the uniform technique have an expected running time of approximately 5 × 10 −6 n log n seconds.

Comparison of beading schemes
We can see from Fig. 17a(top) and 18a that the uniform technique causes a lot of overfills and underfills: on average 1.6 % of the total target area is covered by underfill and likewise for overfill. To our knowledge, the uniform beading scheme, as well as the outer beading scheme, is of little use to FDM printers.
The constant bead count scheme effectively deals with underfills, but generates orders of magnitude more overfills compared to the other schemes. Also, the scheme comes at the cost of greatly varying bead widths and an average bead width that is not close to the preferred bead width. Note that most overfill areas occur near regions of alternating bead width. While the scheme results in short toolpaths, as indicated by the idealized print time, it also results in a wide range of bead widths, which cause the back pressure compensation print time to be very large. See Fig. 18. For an input outline shape which contains both very small and very large features, the constant bead count scheme produces bead widths which can fall outside of the range of manufacturable bead widths. Moreover the centrality marking is not robust against small perturbations in the outline; adding a small chamfer in a corner causes the unmarked ST to be very small at that location, which results in tiny bead widths. See top right of Fig. 17b.
In Fig. 17c we can see that the centered beading scheme effectively deals with overfill and produces desired bead widths in all locations, except for the extrusion paths in the center, where the bead widths range between 0.25w * and 1.8w * . However, it does produce some narrow underfill regions. Compared to the uniform technique the centered technique increases the (open) path count, but considerably reduces over-and underfill and decimates the number of toolpath angles below 45 • . See Fig. 18.
However, according to Fig. 18d the centered scheme exhibits a wider range of bead widths than the distributed schemes: the standard deviation of the bead widths in the centered scheme is approximately 53 µm, while that of the distributed schemes is approximately 23 µm. 1 Moreover, because the quantization operator rounds to the nearest number of beads, in the worst case where we switch from a single to two beads the widths switch from 0.75w * to 1.5w * , which is a considerably smaller range than in the centered scheme. We therefore conclude that the distributed schemes exhibit a lower bead width variation and lower (open) path count compared to the centered scheme. Figures 13, 17d and 17e show that in the inward distributed scheme the outer toolpaths have the preferred width more often than in the evenly distributed scheme, which means that the outline accuracy of the inward distributed beading is less affected by inaccuracies in the adaptive width control system. Furthermore, we find that compared to evenly distributed, the inward distributed scheme produces less corners with angles above 130 • and less overfill, because the area of influence that bead count transitions have is limited in the inward distributed scheme. Thus the inward distributed scheme prevents over-and underfill, generates smooth toolpaths with more homogeneous width and affects smaller more centered parts of the print than the other schemes, while incurring little to no extra print time.

Limitations
Because the performance of the various toolpathing techniques depends on the geometry of a model, they have ramifications for the practice of design for additive manufacturing. Because the naive method produces under-or overfill for parts of an outline with a constant diameter d 2iw * it is best practice to design a model such that horizontal cross-sections have a feature diameter of an even integer multiple i of the bead width. For the centerscheme and for the current state of the art one should only avoid parts for which (2i + 1.8)w * < d < (2i + 0.25)w * in order to avoid underfill. For the distributed schemes however, there is no diameter at which the framework produces under-or overfill for a part with a constant diameter d. The design consideration therefore reduces to limiting the diameter of your parts to be within the range [w min , ∞), where w min is a configurable parameter when using the widening meta-scheme.
The default limit bisector angle α max = 135 • ensures that we don't employ transitioning in shallow wedge regions, which would result in a lot of short odd single bead polylines, which would break up the semi-continuous nature of polygonal extrusion paths; α max = 135 • corresponds to w * / cos 1 /2α max ≈ 0.4 mm long segments and under-/overfill areas of 1 /4(w * ) 2 (tan(α/2) − α/2) ≈ 0.05 mm 2 . However, future work might be aimed at reducing under-/overfill in regions with a low bisector angle without the introduction of short single polyline extrusion segments. If the over-/underfill problem is also solved for non-significant regions we might be able to increase α max and reduce the discontinuity introduced by short extrusion segments.
Another limitation of our method is that in a location v with locally maximal R(v) ≈ (i + 1 /2)w * the odd bead count will result in a single polyline extrusion segment consisting of only a single point. This can be viewed in the bottom right of Fig. 17e for example. In order to print such a dot, we make it into a 10 µm long extrusion segment, with an altered width such that the total volume remains correct. A more principled way of dealing with such a situation remains future work.
Finally it should be noted that although our framework can accurately emulate the constant bead count approach by Ding et al., its emulation of the centered approach by Jin et al. is imperfect. The transitions resulting from out framework introduce sharper corners and there is more width variation in those corners. Whereas the width of the connecting segment in the approach by Jin et al. is the preferred width w * , the bead widths closer to the center resulting from our framework will be twice the local radius, which is larger than w * . However, this inflated bead width variation is expected to have an insignificant impact on the measured bead width variation.

Applications
Toolpaths with varying width is particularly meaningful for narrow parts, since there the negative effect of under-and overfill is more pronounced than in wide parts. In extreme cases, thin features will not be filled at all. Therefore, our framework, while working for wide parts as well, shows most of its potential for objects which contain thin parts. Figure 19 collectively shows the application of the proposed inward distributed scheme for various types of 3D model, including both thin parts (architectural models, casings, embossed text, gears and microstructures) and wide parts (Fig. 19b) and organic shapes (Fig. 19c)).
For architectural models and casings, preventing over-and underfill is expected to make them stronger. For embossed text, preventing underfill reduces the various holes in the top surfaces, which is detrimental to the visual quality of those top surfaces. For gears and similar mechanical parts that are designed with finite element analysis, the less variation in extrusion widths is closer to the assumptions under fast analysis (e.g. using homogenization [44]).
Of particular interest are microstructures that could be uniquely fabricated by 3D printing. For example, topology optimized bone-like structures [45] contain filaments of varying thickness that follow a varying stress direction (Fig. 19g). An angled Gyroid structure with uniform thickness also results in outline shapes with varying width (Fig. 19h). These structures are accurately densely filled using our framework. Another class of microstructures consists of parameterized patterns with varying thickness to achieve functional gradation. Figure 19i shows the contour-parallel toolpaths with varying width of a hexagonal grid neatly switches between different bead counts over the volume, preventing the jagged moves a direction-parallel toolpaths would create for such a case [1].

Conclusion
In this paper we have introduced a framework for computing contour-parallel toolpaths employing adaptive bead width in order to minimize underfill and overfill areas. We introduced beading schemes which improve on the state of the art, and we have introduced a back pressure compensation method for accurate fabrication of adaptive width.
Our framework is flexible, demonstrated by the several beading schemes which emulate existing techniques. The computation times of our framework are on par with the state-of-the-art library for performing offsets of non-adaptive bead width. Our framework is stable: small local changes in the outline shape cause only small changes in the toolpath.
Compared to the state of the art, the inward distributed beading scheme reduces the amount of beads with a width deviating extremely from the preferred bead width by changing the width of several beads near the center instead of only the center-most bead. It is therefore expected to limit the impact of varying the bead width in terms of production accuracy and homogeneity of material properties, which in turn is helpful to efficiently simulate an FDM manufactured part.
The proposed beading scheme greatly improves the process planning for parts with thin contours, which often occur for example in architectural models, prototypes for casings or microstructures. Meanwhile it leaves most of the toolpaths the same as the uniform width technique in large features, meaning that existing studies which relate process parameters with mechanical properties of the print are still applicable. Compared to the naive approach of constant width toolpaths our beading scheme is expected to improve the stiffness, dimensional accuracy and visual qualities of the manufactured model. It is expected that as distributed beading schemes are implemented in commercial software packages and bead width variation control become commonplace, the practice of design for additive manufacturing can disregards some of the nozzle size considerations.  Fig. 19. Visualization of the widths for the output toolpaths of the inward distributed beading scheme (N = 3) applied to various example application objects. From left to right and top to bottom: a house, a case for electronics, a statue, two common logos, a gear, a topologically optimized bone structure, a tilted homogeneous gyroid structure and a heterogeneous thickness hexagonal grid.

C. Accuracy calculation
In order to estimate the overfill and underfill, we need to accurately calculate the area covered by a single extrusion path. If we would simply use an isosceles trapezoidal area, we would get overfill artifacts at corners in the toolpath (Fig. 20a). We therefore use a semi-circle (Fig. 20b) with a diameter equal to the starting width in the one end of each segment, and exclude it at the other end, because it will be included in the next segment. For polyline extrusion paths which are not closed, we also include the semi-circle of the destination location (Fig. 20c).
Using boolean operations we can obtain the polygonal regions for overfill and those for underfill. In order to deal with rounding errors we perform a morphological close of 5 µm, before calculating the total area in mm. We also calculate regions which are covered thrice by different extrusion segments and add twice its area to the total overfill area amount.