Efficient slicing of Catmull–Clark solids for 3D printed objects with functionally graded material

Computers & Graphics 82 (2019) 295–303

Contents lists available at ScienceDirect

Computers & Graphics journal homepage: www.elsevier.com/locate/cag

Special Section on SMI 2019

Efficient slicing of Catmull–Clark solids for 3D printed objects with functionally graded material Thu Huong Luu a,b,∗, Christian Altenhofen a,b, Tobias Ewald a,b, André Stork a,b, Dieter Fellner a,b,c a b c

Interactive Graphics Systems Group, TU Darmstadt, Fraunhoferstr. 5, Darmstadt 64283, Germany Fraunhofer IGD, Fraunhoferstr. 5, Darmstadt 64283, Germany Institute of Computer Graphics and Knowledge Visualization, Graz University of Technology, Inffeldgasse 16c, Graz 8010, Austria

a r t i c l e

i n f o

Article history: Received 19 March 2019 Revised 22 May 2019 Accepted 23 May 2019 Available online 28 May 2019 Keywords: Subdivision volumes Catmull–Clark solids Functionally graded materials Additive manufacturing Slicing

a b s t r a c t In the competition for the volumetric representation most suitable for functionally graded materials in additively manufactured (AM) objects, volumetric subdivision schemes, such as Catmull–Clark (CC) solids, are widely neglected. Although they show appealing properties, efficient implementations of some fundamental algorithms are still missing. In this paper, we present a fast algorithm for direct slicing of CCsolids generating bitmaps printable by multi-material AM machines. Our method optimizes runtime by exploiting constant time limit evaluation and other structural characteristics of CC-solids. We compare our algorithm with the state of the art in trivariate trimmed spline representations and show that our algorithm has similar runtime behavior as slicing trivariate splines, fully supporting the benefits of CCsolids.

1. Introduction Functionally graded materials (FGMs) or heterogeneous materials have gained great popularity and research interest in the past decades. Unlike conventionally engineered parts that are made from one homogeneous material, functionally graded material objects refer to objects with spatially varying material compositions created by mixing two or more materials inside the object [1]. This allows for introducing material heterogeneities into the design domain, e.g. anisotropic properties. The locally varying material distribution is usually designed to (optimally) achieve some overall physical properties satisfying a certain design objective. FGM objects are believed to possess superior properties to those that can be achieved by objects made from a homogeneous material. In a recent report [2], Kevin Eckes showed that graded material transitions greatly improve the stability of multimaterial objects compared to discontinuous changes in material. The increased availability of multi-material 3D printers is leading to an increased interest in manufacturing FGM objects. Current 3D printers, such as the Stratasys Objet Connex series [3] or the J735 and J750 printers [4], can produce different FGM objects

∗ Corresponding author at: Interactive Graphics Systems Group, TU Darmstadt, Fraunhoferstr. 5, 64283 Darmstadt, Germany. E-mail address: [email protected] (T.H. Luu).

https://doi.org/10.1016/j.cag.2019.05.023 0097-8493/© 2019 Elsevier Ltd. All rights reserved.

© 2019 Elsevier Ltd. All rights reserved.

using PolyJet 3D printing technology. The ProJet [5] and the Jet Fusion [6] are other examples for printers that use a similar inkjet approach. The state-of-the-art capabilities of free-form modeling tools are based on boundary representations (B-Rep), describing only the outer surface of a three-dimensional object. Thus, traditional computer aided design (CAD) systems using B-Reps only support assigning one material to each object or topological lump. Therefore, the current representation schemes for CAD objects are incapable of supporting FGM objects without introducing additional information to describe the internal structure of the model. Many different strategies have been proposed for the specification, design, and manufacturing of FGM objects [7]. Voxels are perhaps the most straightforward way of representing a FGM object, as they correspond very well with the input some 3D printers expect, e.g. a series of bitmaps [3]. Another option for modeling FGM objects are piecewise linear volume meshes [8]. For instance, Jackson et al. [9] model FGM objects by subdividing them into simple cells (finite element meshes), with material properties defined analytically within each cell. Trivariate spline/NURBS approaches [10,11] are also a well-known representation and can be regarded as a direct extension of parametric volumes, with additional material information incorporated in each control point. Being a continuous representation, it can be evaluated at any point in space depending on the resolution of the 3D printer so that no discretization nor conversion is required.

296

T.H. Luu, C. Altenhofen and T. Ewald et al. / Computers & Graphics 82 (2019) 295–303

Another continuous and volumetric representation to store and model material gradients are subdivision volumes [12] – a volumetric extension of subdivision surfaces. Extending the famous Catmull–Clark surfaces introduced in [13] results in Catmull–Clark solids (CC-solids) which exhibit the same appealing properties as its two-dimensional equivalent: Being able to handle topologically complex objects without the need for trimming, while a smooth parameterization across element borders allows for representing seamless material gradation. Following this promising approach, we utilize CC-solids as presented by Joy and MacCracken [14]. However, to our knowledge no one has dealt with the problem of directly slicing CC-solids, which is a crucial step in fabricating objects with additive manufacturing. Recently, a method for limit point evaluation was presented by Altenhofen et al. [15] that efficiently retrieves the Euclidean position and the material composition for a location in the parametric domain. However, to be able to evaluate the slices we have to solve the inverse problem called point inversion: Given a point in Euclidean space we have to determine the parametric values for that point in order to evaluate the material property value assigned to it. In this paper, we tackle the challenge of efficiently slicing CC-solids enabling a fast printing process. We believe Catmull–Clark solids to be a promising representation for FGM design and fabrication. Therefore, realizing fast processing of material evaluation and slicing is one step further to promote the use of CC-solids. The main contributions presented in this paper are summarized as follows: •

•

An efficient algorithm for point inversion in CC-solids optimized for massive data evaluation during slicing. An approach to quickly and efficiently slice the boundary of volumetric subdivision models.

This paper is organized as follows: In Section 2, we provide an overview of the major approaches for current representation schemes for FGM objects. In Section 3, we introduce our algorithms developed to efficiently slice and evaluate CC-solids. In Section 4, we demonstrate the effectiveness of the algorithms using some complex examples and perform a comparison with a slicing approach on trivariate splines. Finally, we conclude the paper in Section 5. 2. Related work Modeling of FGM objects and their properties is a major topic for research. Over the years a lot of work has been done in the areas of modeling and representation of FGM objects, generating a variety of representation forms for functionally graded materials [16,17]. Many different strategies for specification, design, and manufacturing of FGM objects have been proposed. However, no single strategy ever emerged as dominant [7]. In terms of exactness and compactness, existing data representations can be classified into two categories: evaluated and unevaluated representations. Evaluated representations model the object in a discrete and approximate manner, generating deviations from the original geometry. A well-known example is the voxel representation. Voxelbased representations model objects using a three-dimensional voxel grid. Each voxel represents a small cube in space with material values assigned to it. Thus, it is suitable for representing highly heterogeneous objects whose material distributions are strongly irregular [18]. In addition, this representation highly correlates with the manufacturing concept of modern 3D printers, such as the Stratasys J750 printer, which allow for material variation on a pervoxel basis. While being a descriptive representation with a simple data structure, memory and computational requirements make voxel-based representations costly to handle for large objects when

a high resolution is needed [19]. As one solution, Martin and Cohen [20] introduce a sparse hierarchical representation on voxels rather than a dense grid representation. Richards et al. [21] present another approach by describing functionally graded materials using volumetric gradient patterns that drive a process of probabilistic material deposition. Two examples of tools currently on the market that allow for volumetric material editing on a per voxel basis are GrabCAD Voxel Print [22] and Autodesk Monolith [23]. Contrary to the evaluated models, unevaluated models utilize exact geometric data representations and exact functions to represent material distribution. They are generally resolution independent, more compact, concise and mathematically more rigorous providing higher accuracy for geometry and material distributions. One subcategory is control point based representations, where FGM objects are modeled as extensions of parametric curves, surfaces or volumes with additional material information attached to each of the control points. In 2016, Massarwi and Elber [24] introduce a framework for modeling trivariate B-spline models as V-Reps (Volume Representation analogously to the commonly used B-Rep). On the basis of Massarwi and Elber’s representation, Ezair et al. [11] present an efficient method that enables the direct slicing of functionally graded material objects described as V-Reps and demonstrated their capabilities by fabricating objects using a modern multi-material 3D printer. Dokken et al. [10] studied and compared trivariate hierarchical B-splines, T-splines, and locally refined (LR) B-splines for representing 3D models in CAD and AM. However, trivariate splines require a tensor product topology consisting of regular grids connected through compatible interfaces. This drawback restricts the class of input geometry to purely hexahedral meshes with regular topology. Complex objects featuring holes or requiring irregular topologies can only be represented by either trimming the trivariate cell with a trimming surface or by splitting the domain into multiple blocks as shown by Dokken et al. [10]. As a result, these trivariate spline-based representations only provide C0 continuity across cells [24]. Higher continuity between elements can be constructed algorithmically, but implies high computational effort, especially when modifying the geometry. In contrast, subdivision, the de facto standard for computer animation [25], is able to model a smooth shape with arbitrary topology using just a single control mesh, thus providing a water-tight representation. Volumetric subdivision methods, such as the ones presented by Joy and MacCracken [14], Schaefer et al. [26], Bajaj et al. [27] and Chang and Qin [28] transfer the concept of surface subdivision to volumetric meshes. For Catmull–Clark solids, the underlying basis functions inherently provide higher continuity across element borders – C2 away from extraordinary locations, C1 at extraordinary vertices and edges1 . This smoothness and shape property of the basis functions is reflected in two characteristics of CC-solids relevant for 3D printing of FGM objects: It enables the representation of smooth material distributions and facilitates the modeling of organic shapes. Altenhofen et al. [12] already suggested a novel approach for representing material gradients in volumetric models using Catmull–Clark solids, supporting continuity and providing ways for interactive modeling with locally varying properties. They also present a first printing pipeline for fabricating FGM objects described by subdivision solids. However, when transferring the model to the printer, they perform a conversion to a tetrahedral mesh causing again the same problems known for discretized representations, such as the trade-off between number of elements

1 To the best of our knowledge, there is no proof for C1 smoothness covering all extraordinarities. However, we have never experienced any breakdown of regularity in practice.

T.H. Luu, C. Altenhofen and T. Ewald et al. / Computers & Graphics 82 (2019) 295–303

297

Fig. 1. Slicing pipeline: We start with a model and slice its boundary (meaning its outer surface) by intersecting it with a given slicing plane. The results are one or multiple contour lines in the current slice. Material information is then calculated for the slice by rasterizing the plane and computing material properties for each grid point inside the volume.

and precision. An efficient direct slicing algorithm for CC-solids is yet unpublished. 3. Efficient slicing of Catmull–Clark solids Following Altenhofen et al.’s approach [12], we model a FGM object as a CC-solid by assigning not only its 3D position but additionally a vector of material properties to each control point in the volumetric subdivision mesh. The mesh is defined as a cell complex consisting of vertices, edges, faces, and cells. Since CC-solids are generated by an iterative process refining the control mesh successively, every point in the limit volume V ⊂ R3 carries material information. We exclude the case of self-intersecting volumes since each location must have a unique material property. The overall slicing process consists of the following three main steps: 1) Slicing of the limit surface of CC-solids to extract outer contours (Section 3.1) 2) Rasterizing the slicing plane with a given printing resolution and using the outer contours to determine occupancy information, meaning whether a voxel is inside or outside the volume (Section 3.2.1) 3) Determining material information for every voxel inside the volume via point inversion (Section 3.2.2) Fig. 1 gives an overview of the slicing pipeline: For every single step, we show the results when processing an engine block model. For step 3), we need to directly evaluate the CC-solid limit. For an efficient implementation, we use the limit evaluation approach presented by Altenhofen et al. [15], which enables us to evaluate sample points in cells with extraordinary constellations in constant time. The alternative would be to subdivide the mesh until the sample point is covered by a regular configuration and utilize trivariate B-Spline evaluation. But with voxels having a dimension of 0.04 millimeters, a common printing resolution for current 3D printers, this might require up to approximately 10 to 15 subdivision steps depending on the model size, which cannot be handled efficiently. In order to make use of the constant-time evaluation, we make two assumptions on the initial control mesh: After a finite number of subdivision steps, •

•

the subdivided control mesh consists only of hexahedral cells and each cell of the subdivided control mesh contains at most one extraordinary vertex (EV) and all adjacent extraordinary edges (EE) must meet there. For CC-cells, an EE is defined as an edge incident to less or more than four hexahedra. An EV is defined as a vertex where more than two EEs meet.

Denoting the subdivided mesh by M and following [15], we parameterize every cell Ci in M in its local coordinate system

Fig. 2. Each cell Ci has a geometric transition function geoi from its local parameter space (u, v, w ) ∈ [0, 1]3 to its limit volume defined in Euclidean coordinates (x, y, z). For simplicity, only the 2D case is depicted.

(u, v, w ) ∈ [0, 1]3 . This provides a geometric transition function geoi : [0, 1]3 → V for each cell, see Fig. 2. In the same manner, we obtain a material transition function mati for every cell, again parameterized over (u, v, w ) ∈ [0, 1]3 . We collect all these parameterizations in (global) transition functions geo and mat mapping (i, u, v, w ) to geoi (u, v, w ) = (x, y, z ) and mati (u, v, w ) = m respectively, where i indexes the cell, u, v, w ∈ [0, 1] are parameter values on the ith domain cell, x, y, z are 3D coordinates of the corresponding limit point in the volume V, and m is a vector describing material properties at the respective limit point. It has to be noted, that for cells containing an extraordinary vertex, the coordinate system is chosen such that the EV has the parametric value (0, 0, 0). Thus, the orientation of the parameter space may change from cell to cell, which has to be considered during the slicing process. 3.1. Slicing the limit surface of Catmull–Clark solids In CAD, intersection algorithms are widely applied to many problems such as surface trimming and path generation. In most cases, there are no analytic solutions to such intersection problems. Therefore, we developed an efficient numerical algorithm for the demands of our scenario. The first step of our efficient slicing approach for Catmull–Clark solids is to compute the contours of each slice, i.e. the intersection curves between the subdivision limit surface of the model to be sliced and the given slicing plane. The calculation of the contours is based on the ideas presented in the papers of Zhu et al. [29] and Severn and Samavati [30]. Both papers present similar approaches for calculating an intersection curve between two subdivision surfaces. Since the limit surface of a Catmull–Clark solid is equivalent to the limit of a Catmull–Clark surface calculated from the boundary control points of the CC-solid, we can reduce the problem to calculating the intersection curve between a Catmull–Clark limit surface and a plane. In contrast to the approaches by Zhu et al. [29] and Severn and Samavati [30], our approach allows for different subdivision

298

T.H. Luu, C. Altenhofen and T. Ewald et al. / Computers & Graphics 82 (2019) 295–303

Fig. 3. Slicing contours for different objects. Depending on the shape of the object, our approach creates multiple independent contours (b). The slicing contours are shown in green. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4. Limit patches are subdivided locally until they fulfill the flatness criterion and can be approximated by a discrete quadrilateral. The number of local subdivision steps varies, based on the second derivative of the limit patch. The limit surface is shown in yellow. The intersection curve is shown in green. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

levels for each individual face. This way, each face is only subdivided adaptively as often as needed, which reduces the number of patches which need to be considered during calculations and improves the performance of the algorithm, especially for models with highly varying curvature (see Fig. 4). A second difference to both papers lies in the form of data representation and tracking. At the start, all faces are given as the boundary faces of the volumetric subdivision model. Since, for each face, every subdivision step is performed locally (each resulting in four sub-faces), topological information gets lost during these subdivisions. For example, it is not clear, which intersected faces lie in the neighborhood of which. To keep track of this information, we assign each subface a unique patch index. Its construction will be described later on in more detail. Please note, that although our approach creates discrete intersection curves, we ensure that all intersection points are on the limit surface and store them in u, v coordinates of their corresponding surface patch. The discretization error results from connecting neighboring intersection points with straight lines instead of exactly following the limit surface. The amount of error can be controlled by the flatness criterion introduced in Section 3.1.1. Fig. 3 shows slicing contours for an engine block model as well as for a coyote model. If the given subdivision model has concave regions or holes, multiple independent contours might result from the intersection with the slicing plane as shown in Fig. 3(b). Our approach for calculating the slicing contour consists of the following two steps: 1) Determining all limit patches intersected by the plane 2) Calculating the intersection curve(s) across these limit patches 3.1.1. Step 1: Determining intersected limit patches In analogy to the approaches presented by Zhu et al. [29] and Severn and Samavati [30], we first gather all limit patches that will potentially intersect with the slicing plane. As a Catmull–Clark limit patch is bounded by the convex hull of its local control mesh,

Fig. 5. Numbering scheme for sub-patches to capture neighborhood information. Indices are concatenated when preforming local subdivision steps.

we select potentially intersecting limit patches by testing this convex hull for intersection with the slicing plane. For every potentially intersected patch, the local control mesh is subdivided recursively until the limit patch of the subdivided mesh meets a certain flatness criterion: A limit patch is considered flat, if the norm of its second derivative in u and v direction is below a given threshold. As the derivatives of Catmull– Clark limit patches diverge at (u, v ) = (0, 0 ) for irregular topologies, we evaluate the second derivatives across the diagonals of the patch at (u, v ) = (0, 0.5 ), (u, v ) = (0.5, 1 ), (u, v ) = (1, 0.5 ), and (u, v ) = (0.5, 0 ). Because of the smoothing behavior of Catmull– Clark subdivision, each subdivision step increases the flatness of a patch. Therefore, after a finite number of subdivision steps, every patch meets the flatness criterion. Once a locally subdivided limit patches is considered flat enough, it is approximated by a discrete quadrilateral defined by the patch’s limit points at (u, v ) = (0, 0 ), (u, v ) = (0, 1 ), (u, v ) = (1, 0 ) and (u, v ) = (1, 1 ). To determine if this locally subdivided limit patch actually intersects with the slicing plane, its corresponding quadrilateral is tested for intersection. Fig. 4 visualizes these quadrilaterals together with the limit surface and the resulting intersection curve. During the process of local subdivision, the relative position of a subdivided patch inside its parent patch is recorded. As every patch subdivides into four sub-patches, we employ a numbering scheme of k = 0, 1, 2, 3 in analogy to Stam’s numbering scheme for sub-patches presented in [31]. Concatenating these indices k over all local subdivision steps results in a combined patch index as shown in Fig. 5. Binary encoding of these indices allows for efficiently determining neighboring patches as explained in Section 3.1.2. Finally, every locally subdivided limit patch that fulfills the flatness criterion and actually intersects the slicing plane is aggregated for calculating the actual intersection curve(s) (see Section 3.1.2). 3.1.2. Step 2: Calculating the intersection curves In the second step, the intersection curves are calculated on the basis of the intersected patches from step 1. This is done by calculating the intersections of each patch with the slicing plane and “stitching” these separate intersections together to obtain a set of intersection curves. As stated above, all limit patches collected in step 1 fulfill our flatness criterion and can therefore be substituted by their corresponding quadrilateral. In general, the local intersection line between a quadrilateral and a plane is a straight line defined by two intersection points on two different edges of the quadrilateral. Therefore, for every edge of the quadrilateral, we calculate the intersection point with the slicing plane. As this intersection point is not guaranteed to lie on the actual limit surface, we use the point as a starting point for iteratively tracing the exact intersection point on the corresponding limit patch. Using Newton’s method, we advance along each edge of the limit patch in u, v space to find the limit point on the

T.H. Luu, C. Altenhofen and T. Ewald et al. / Computers & Graphics 82 (2019) 295–303

plane closest to the intersection point. Since we move along an edge, only one parameter u or v is varied (depending on which edge we are on), while the other one is fixed. For a given patch, every intersection point found on one of its edges is at the same time an intersection point of the neighboring patch that shares this specific edge. For moving from one locally subdivided limit patch to its neighbors, the information encoded in its patch index is used. Patches that share an edge on the same subdivision level also share one bit in the binary representation, diagonal patches do not. When moving to another patch on a different subdivision level, the u, v coordinates of the intersection point have to be scaled accordingly. When traversing across the borders of the original unsubdivided patch, the binary patch indices have to be rotated such that the u, v coordinate spaces of the current patch and the next one are aligned. Starting from an arbitrary patch, the intersection curves are calculated by traversing from the patch over its neighbors in a depth-first manner. Finally, the u, v coordinates of all intersection points have to be transformed from the u, v space of their locally subdivided limit patch to the u, v space of the original unsubdivided limit patch they belong to. Here, again, the information encoded in the binary patch index is used. Ultimately, an intersection curve between the Catmull–Clark limit surface and the slicing plane is defined as an ordered sequence of discrete intersection points, each one being defined by a patch of the CC-surface control mesh and the parameters (u, v ) of the limit point.

3.2. Evaluating the slicing plane Given the slicing contour the next step is to evaluate the area inside the contour to derive actual material slices. This evaluation consists of rasterizing the slice with a regular grid specified by the printing resolution and retrieving material information for every voxel inside the model. We will describe both procedures in more detail in the following sections. Please note, although we store and calculate the information in a 2D representation namely a grid of pixels, each pixel corresponds to a voxel in the volumetric object. Thus, in the following sections we will use the term voxel even in the context of a 2D slicing plane in order to keep awareness of the volumetric nature of the calculated information.

3.2.1. Rasterization The rasterization process generates an occupancy bitmap that states which voxels in the slicing plane are inside the volume and which are not. By first determining the occupancy information for a given slice, we avoid attempting to determine material information for voxels outside the object. We determine the size of the generated bitmaps based on the extent of the mesh’s axis-aligned bounding box. Dimensions in x and y direction influence the number of voxels in either direction. For a given layer thickness dz , the dimension in z sets the number of slices needed to fabricate the whole object. Given the bitmap size, we sample the slicing plane in a regular grid. As each sample point P depicts the center of a voxel, the distance between sample points is given by the voxel’s dimensions. To determine whether a sample point is inside or outside the object we use an efficient winding number approach [32], where the slicing contour obtained in Section 3.1 is interpreted as a planar polygon with closed boundary edges. The implementation is as fast as the simple boundary crossing algorithm, which counts the number of times a ray starting from the point P crosses the polygon boundary edges. Furthermore, our algorithm gives the correct answer for non-trivial polygons (i.e., with self intersections), where the boundary crossing algorithm would fail.

299

3.2.2. Point inversion As a result of the rasterization procedure, we get a set of voxels described as points in Euclidean coordinates located in the volume for which the corresponding material properties have to be computed. Our approach calculating the material properties consists of: •

•

Point inversion: Translate each point P at position (x, y, z) back to its (unique) parametric value (i, u, v, w ) via inversion of the geometric transition function geo. Evaluate the material transition function mat (i, u, v, w ) = m to determine the material properties m for P.

Since explicitly and exactly inverting geo is impossible, an established solution for point inversion is to approximate (i, u, v, w ) by estimating an initial parameter value and improving it recursively via numerical tracing approaches, e.g. Newton’s method. This procedure was already used for B-spline curves and surfaces in [33], where much effort was spent for determining a good parametric value already located in the correct parameter segment. A similar procedure can be used for the three-dimensional case of CC-solids. However, since the search for the correct parameter segment is very time-consuming in those algorithms, we propose an optimized procedure. The idea is to quickly generate a list of cells potentially containing the correct parameter value, rather than spending much time finding the correct cell immediately (see line 1 of Algorithm 1). In Algorithm 1 Pseudo-code explaining the approach for point inversion. 1: Generate a list of cells C j potentially containing the unknown value (i, u, v, w ) 2: for each cell C j do Compute a starting value p0 := (u0 , v0 , w0 ) 3: k←0 4: 5: while |(x, y, z ) − geo j (uk , vk , wk )| > thres do pk+1 ← NewtonStep(geo j , pk ) 6: if StoppingCriteria( pk+1 , pk ) then 7: Continue for-loop with next cell 8: 9: end if k←k+1 10: end while 11: return ( j, uk , vk , wk ) 12: 13: end for addition, by treating each cell separately we do not have to consider the orientation change of the locally defined parameter space when transitioning from one cell to another. These cells are then traversed to find the cell actually containing the parameter value to P. The supposed extra effort of investigating several cells is bypassed, if the algorithm successfully terminates on one of the first entries in the list. To place the correct cell at one of the first entries, we sort the list according to some heuristics. How to generate and arrange the list of cells, how to compute the starting values (line 3) and details on the Newton method (line 5) will be explained in the following sections. In addition, we describe an optimization procedure which utilizes spatial coherence to further speed up the slicing process. Cell list generation and computation of starting parameters. Generation of a list with potential cells Cj is a crucial step for the overall performance of the point inversion algorithm. For a cell Cj we denote the set of control points belonging to Cj ’s one-ring neighborhood by CPj , see the green colored points in Fig. 6. A feature of CC-solids is that the limit of Cj lies within the convex hull of CPj . Hence, checking if the point P is inside the convex hull of CPj generates a list of potential cells Cj guaranteed to contain the correct

300

T.H. Luu, C. Altenhofen and T. Ewald et al. / Computers & Graphics 82 (2019) 295–303

Fig. 6. 2D profile of a control mesh, where each quad represents a hex cell. The black dot depicts the point P at position (x, y, z). Green: Control points CPj of a potential cell Cj with bounding box as dashed green line. Blue: Control points of Cj with inner bounding box as dashed blue line. Potential cell list L highlighted in orange and sublist l highlighted in blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

cell. To reduce runtime, we choose to check against axis aligned bounding boxes of CPj (shown as a dotted green rectangle in Fig. 6) and call the corresponding cell list L (highlighted in orange). Since the axis aligned boxes comprise the convex hulls, the list L is also guaranteed to contain the correct cell. However, we observed that a large portion of the limit of Cj is already contained in the bounding box of control points only incident to Cj . In Fig. 6, those control points are shown as blue framed circles and their corresponding bounding box is shown as a dotted blue rectangle. Thus, we additionally filter a sub-list l out of L with respect to P being located in the inner bounding box of Cj (consisting only of the blue cell in Fig. 6) to have a more narrowed list. Given the list of cells, we now need to determine starting values for each cell. In line 3 of Algorithm 1, the starting value p0 := (u0 , v0 , w0 ) for each cell Cj is computed as follows: We sample the parameter domain [0, 1]3 with a uniform 4 × 4 × 4 grid {q1 , . . . , q64 }, map it via geoj to geometry space and choose p0 as the grid point with minimal Euclidean distance to (x, y, z)

d j := |(x, y, z ) − geo j ( p0 )| = min

r=1,...,64

|(x, y, z ) − geo j (qr )|.

Since dj should be minimal for the correct cell, a cell list sorted w.r.t. increasing dj should have the correct cell early on. Based on this observation, we propose a heuristic to reduce the number of transition function evaluations by first computing starting values merely for l. Only if point inversion with l sorted w.r.t. increasing dj fails to compute the correct parameter value, we calculate starting values for the whole list L and run the point inversion algorithm with L sorted by dj . Newton’s method. Newton’s method starts in line 5 which is a numerical check if ( j, uk , vk , wk ), with k being the current Newton iteration, maps to (x, y, z) up to a tolerance thres. In our implementation we set thres = 0.001 millimeters. In general, it is sufficient to set thres to half of the minimal voxel dimension in x or y direction, because we do not need to approximate more precisely than the printer resolution. However, we decided to choose a more restrictive threshold since Newton’s method converges very quickly, thus giving us a more precise evaluation result. The actual update pk+1 = pk + k in line 6 is given by solving the system of linear equations

Jgeo j ( pk ) · k = (x, y, z ) − geo j ( pk ),

where Jgeoj (pk ) is the Jacobian of geoj computed at the parametric value pk . We make use of the open-source C++ library Eigen [34] to solve the system. Since [0, 1]3 is the parameter domain of geoj , we clip the output pk+1 of NewtonStep to lie in this domain. The boolean function StoppingCriteria in line 7 then examines if the Newton iteration reaches a stationary setting by checking the step size both in parameter and geometry space. In addition, to avoid an endless loop during point inversion, we restrict the maximal number of iterations to 10. This is a conservative estimate since experimental results showed that given the correct cell, the algorithm usually converges in less than five iterations. Hence, terminating each cell-loop by the StoppingCriteria is the only possibility left for Newton’s method to stop without any result. Because of the clipping in NewtonStep, this can either happen if (i, u, v, w ) is not in the current cell or if the Newton iteration converges to a wrong parameter value due to local extrema. However, we have never experienced the latter case when using CC-solids. We like to point out, that apart from using PI to determine material information for voxels lying inside the object, it can also be used to determine if a voxel is outside the object. If for (x, y, z) no potential cells are found or all cell-loops are terminated by the StoppingCriteria, we conclude that the point is outside the object. Optimization by utilizing coherence. Algorithm 1 describes the overall approach how to determine the material information for one point if no prior knowledge is given. The most time-consuming step in the algorithm is line 1 where we have to iterate over all cells and determine their bounding boxes. However, since we utilize a numeric tracing approach, we can avoid the costly step in line 1 by exploiting information from neighboring voxels which were already evaluated. In order to accelerate the computation within one slice, we implemented a tracing optimization method similar to the one by Elber et al. [11], where any adjacent voxels on the current or previous line (including diagonally placed voxels) are used as an initial value for tracing. This reduces the list of cells Cj to a list with only one entry, namely the cell of one previously calculated voxel. Additionally, we can skip the computation of a starting value in line 3 by using the already determined (u, v, w ) parameters of the voxel as starting point. Only if no evaluated neighboring voxels are given or the point inversion with all neighboring voxels as starting points failed because we iterated in wrong cells, the whole point inversion algorithm is invoked. We will show in Section 4 that this optimization significantly reduces the number of executions of the whole algorithm and thus overall runtime. 3.2.3. Output As output, we generate slices with material information for every single voxel. This information could consist of visual properties, such as color, or physical material properties, such as elasticity parameters or porosity. However, this kind of information cannot be used as-is as input for 3D printers. Many printers, especially the ones focused on PolyJet 3D printing technology, only accept a single material for each voxel. Therefore, we have to convert the property information into a form usable as printer input. For colored voxels, we utilize the generic printing driver software Cuttlefish [35] as described in Altenhofen et al.’s approach [12] to perform this conversion. The colors are half-toned using a 3D blue noise mask, which results in a concrete assignment of primary materials. Alternatively, we can apply error diffusion dithering algorithms such as [36] or [37] for every slice to perform this conversion. For material properties such as elasticity/stiffness, the property values have to be mapped to existing printing materials or material mixtures available to the AM process of choice. This mapping can be derived from printed probes and measurements which is a topic of future research.

T.H. Luu, C. Altenhofen and T. Ewald et al. / Computers & Graphics 82 (2019) 295–303

301

Table 1 Statistics on evaluating the slice highlighted in Fig. 9(a) with and without utilizing coherence.

Fig. 7. Visualization of all test models. From top left to bottom right: topology optimized, grid, coyote, engine block, ring.

4. Experimental results To show the efficiency of our algorithms, we performed experiments with different functionally graded material objects modeled with CC-solids. Additionally, we compare our algorithm with the state-of-the-art approach in trivariate trimmed spline representations [11]. We implemented the algorithms in C++ and ran our tests on an Intel Core i7-7700K (4 Cores, 4.2 GHz) Windows 10 machine. In the measurements, we set the sampling resolution to the printing resolution of the Stratasys J750, which is 600 × 300 DPI (0.0423 × 0.0846 millimeters) in xy direction, and 0.03 millimeters in z direction. For representing Catmull–Clark solids, we use a halfface data structure similar to OpenVolumeMesh [38]. Fig. 9 shows fabricated objects that were created using the methods outlined in this paper. These objects were printed on a Stratasys J750 printer with their GrabCAD Voxel Print interface [22]. We start by investigating how much our point inversion (PI) approach benefits from using spatial coherence as described in Section 3.2.2. In Table 1, we compare the algorithm using full PI all the time with the optimized version exploiting coherence when applied to the slice highlighted in Fig. 9(a) consisting of 364,768 voxels: We examined the time (overall and separated by the three main steps) needed to determine materials for all pixels in that

Optimization Portion of full PI [%]

None 100%

Utilizing coherence 0.2273 %

Process step

Time [s]

Percentage [%]

Time [s]

Percentage [%]

Contour slicing Rasterization Material evaluation Total

0.064 5.62 50.70 56.81

0.11 9.89 89.24

0.064 5.62 9.38 15.08

0.42 37.26 62.20

slice and the portion of voxels for which we had to invoke full PI. As expected, utilizing coherence provides a noticeable speedup, reducing runtime for material evaluation to a fifth and the amount of full PI calls to 0.2273% of all voxels. We like to point out that the time saved is not directly proportional to the reduced amount of full PI calls. This can be explained by the effect when a voxel and its predecessor are not located in the same cell forcing the coherence optimized algorithm to perform a multitude of Newton iterations in a wrong cell before calling full PI. In the present example, there were 17,645 point inversion calls in wrong cells corresponding to 4.6% of all calls. We assess the quality of our material retrieval results by looking at the average deviation of the point inversion result to the voxel center (x, y, z). Since the same basis functions are used for material and geometry modeling, we assume that a small geometric error implies a small error in determining materials. We performed an exemplary evaluation on the coyote model. In our experiment, the point inversion method was performed successfully for all pixels with an average deviation of 0.0 0 06 mm, which is even less than our threshold thres = 0.001 mm and far below half of the diagonal of a voxel. Thus, we conclude that our algorithm can voxel-accurately retrieve material information from CC-solids. In order to investigate the overall runtime of the coherence optimized algorithm, we introduce some test models shown in Fig. 7. The ring model, see Fig. 8(c), exhibits a regular structure without any EEs or EVs possessing a graded material distribution from the inner ring surface to the outer ring surface. All other figures are

Fig. 8. Cut through different models. The slicing colors may map to a physical property like stiffness or could be interpreted as a mixture between two materials. Slicing contours are shown in green. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

302

T.H. Luu, C. Altenhofen and T. Ewald et al. / Computers & Graphics 82 (2019) 295–303

Fig. 9. Fabricated objects with graded materials based on volumetric models described as Catmull–Clark solids. Table 2 Statistics for the models shown in Fig. 7. Times are measured with a parallel implementation utilizing four threads. Model

#Initial cells (#Cells in M)

#EE / #EV

Bitmap size

#Slices

#Voxel in model [ × 106 ]

Total time [h]

Average #Voxel per slice

Average time per slice [s]

Coyote Engine block Ring Topology optimized Grid

48 (2832) 301 (19,328) 512 (512) 163 (1304) 10 0 0 (80 0 0)

177 / 8 360 / 8 0/0 76 / 7 5040 / 1539

(1494 × 991) (1090 × 547) (1975 × 987) (3397 × 478) (1335 × 357)

738 688 576 142 866

218.87 138.31 506.77 57.12 156.01

1.17 0.68 1.65 0.68 2.34

295,608 200,902 876,759 400,615 180,022

5.71 3.54 10.31 17.08 9.73

more complex in the sense of a larger amount of extraordinary edges and vertices. The coyote model in Fig. 9(a) even contains a few wedges and tetrahedral elements next to hexahedral cells. The same complexity holds for the engine block model (Fig. 1) which also has six cylindrical holes with the cylinder liners having a different material property than the rest representing hardened material to withstand higher stresses. The topology optimized model in Fig. 9(b) exhibits a material distribution of one material in the center fading out to the sides. The grid model (Fig. 8(d)) exhibits the highest geometric complexity featuring a fine-grained regular lattice structure with a multitude of holes. Table 2 shows some statistics for the presented models. As explained in Section 3, we assume that each mesh only consists of hex cells. Thus, we subdivide the initial mesh resulting in a mesh M with more cells and control points. Table 2 shows the number of cells before and after refinement together with the number of EEs and EVs in the subdivided mesh M for every model. In addition, we list runtimes of the coherence optimized slicing algorithm. All slices are calculated independently, meaning they do not rely on information from previous slices. This allows for parallel processing of multiple slices. The times shown in Table 2 are measured with a parallel implementation utilizing four threads. The time needed to generate the bitmaps generally depends on the size of the model, more specifically on the number of printer voxels inside the model and the number of mesh cells. The more cells the model has, the higher the percentage of full PI calls and the lesser we can profit from spatial coherence, see the number of full PI calls for the engine block model. In addition, our experiments revealed that for complex contours the rasterization process starts to dominate the runtime as suggested by the timings of the topologically optimized model. Next, we compare the performance of our algorithm with the state-of-the-art slicing approach for trivariate trimmed spline representations by Ezair et al. [11]. Since their example models are not freely accessible and most of our models are too complex to be modeled straightforwardly as a set of trivariate cells, we can only perform rough comparisons between models with approximately similar complexity. Such a comparison is hampered even more since a cell in the case of Ezair et al. corresponds to a trimmed trivariate tensor product spline with an arbitrary number of parameter segments, while a CC-solid cell corresponds to a single parameter segment. Thus when comparing models in different representations, the number of CC-solid cells and trimmed spline

Whole PI [%] 0.4398 1.2775 0.0755 0.3586 2.1766

Table 3 Comparison of statistics of the models fabricated in [11] and our test models. Model

#Cells

Average #Voxel per slice

Coyote Engine block Ring Topology optimized Grid Chalice [11] Utah teapot [11] Duck [11]

2832 19,328 512 1304 80 0 0 5 6 1

295,608 200,902 876,759 400,615 180,022 235,110 255,924 921,864

Average time per slice [s] 5.71 3.54 10.31 17.08 9.73 5.75 3.14 10.4

segments should match, but the latter one is unknown for the models of Ezair et al. [11]. Nevertheless, Table 3 shows some statistics for the models sliced and fabricated in [11] next to our models. First, we compare our ring model with the duck model from [11]. Both are of comparable simplicity since the ring model can be represented by only one trivariate spline just like the duck model. In average, Ezair et al. create 90 0,0 0 0 material entries in the bitmap representing a slice in about 10 s. Our approach creates for the ring model 870,0 0 0 bitmap entries for which we also need 10 s. Thus, our slicing algorithm performs slicing in the same order of magnitude for regular cases. Comparing the other models directly is more difficult because they have different – if not incomparable – complexities in case of extraordinary constellations. However, we consider the number of voxels per slice and the associated slicing duration as a comparison criterion. The chalice, the teapot, the coyote and the engine block model exhibit a similar number of voxels per slice ranging from 20 0,0 0 0 to 295,0 0 0. The same accounts for their slicing time which ranges from 3.14 s to 5.75 s which shows that our slicing algorithm yields similar runtimes as the approach of Elzair et al. To establish an even more direct comparison between both technologies, we use the total number of voxels in a model and the overall runtime of the slicing algorithm to compute the average time needed for evaluating a single voxel. Liberating from the influence of the model complexity, we average this runtime per voxel over several models to account for a value only based on the algorithm. Using the Chalice and the teapot model, the algorithm of Ezair et al. takes around 1.8363 × 10−5 s to perform one point inversion, while our approach approximately needs 1.8468 × 10−5 s for the coyote and the engine block. Assuming

T.H. Luu, C. Altenhofen and T. Ewald et al. / Computers & Graphics 82 (2019) 295–303

a reasonably sized object consisting of 500 million voxels, the slicing approach by Ezair et al. would take 9181.47 s, whereas our approach needs 9234.16 s. Since the difference of 52.69 s (corresponding to 0.56% of the overall slicing time) is neglectable, we conclude that we perform the process of slicing in the same order of magnitude as the state-of-the-art approach by Ezair et al. We deliberately exclude the topology optimized model and the grid model from the analysis since their complexity is not matched by any example provided in [11]. Thus, including the models would result in a distorted comparison. Despite the complexity, our algorithm is able to slice the models in a reasonable amount of time. 5. Conclusion We presented an efficient solution for slicing of Catmull–Clark solids to realize the calculation of slices for additive manufacturing. This involves an efficient model-plane intersection method to determine slicing contours and an algorithm to perform point inversion for CC-solids retrieving parametric values for Euclidean points in the slicing plane. We showed in experiments that we can efficiently slice functionally graded material objects described as CCsolids with varying geometric constellations and complexities. By exploiting CC-solid characteristics such as the convex hull property and exploiting spatial coherence, we achieve similar runtime behavior as a state-of-the-art approach for trivariate splines while fully supporting the benefits of CC-solids. In future, we plan to perform a quantified error analysis to assess the quality of the queried material information of Algorithm 1. In order to speed up the slicing process, we plan to implement socalled cell marching to be able to determine the neighboring cell the sample point marches into during point inversion. In addition, combining the winding number approach with a flood-fill algorithm would result in an additional speed-up for rasterization calculations, which are a bottleneck for models with complex slicing contours. Another idea to speed up computation is to use results from a previous layer as initial values for contour slicing and material tracing. Transferring individual building blocks of our approach to the GPU would also significantly accelerate calculations. To accomplish this, we will make use of efficient GPU data structures for volumetric meshes (e.g. the one presented by Mueller-Roemer et al. [39]). Declaration of interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgement This work was partially funded by the Fraunhofer-internal project “futureAM” (project number: 600288), as well as by the German Federal Ministry of Education and Research (BMBF, funding number: 01IS17047) and by the European Community under grant agreement 825030 (QU4LITY). References [1] Birman V, Byrd LW. Modeling and analysis of functionally graded materials and structures. Appl Mech Rev 2007;60(5):195–216. [2] Aerosint. How to make cheap, scalable multi-material 3d printing a reality. 2018. https://medium.com/@aerosint/how-to-make-cheap-scalablemulti- material- 3d- printing- a- reality-d298fa1f76da. [3] Doubrovski E, Tsai EY, Dikovsky D, Geraedts JM, Herr H, Oxman N. Voxel-based fabrication through material property mapping: a design method for bitmap printing. Comput Aided Des 2015;60:3–13. [4] Stratasys Ltd.. Stratasys j750. 2019. https://www.stratasys.com/3d-printers/ j735-j750. [5] 3D Systems, Inc.. Projet mjp 5500x. 2017. https://www.3dsystems.com/ 3d- printers/projet- mjp- 5500x.

303

[6] HP Development Company. Jet fusion. 2017. http://www8.hp.com/us/en/ printers/3d-printers.html. [7] Qin Y, Qi Q, Scott PJ, Jiang X. Status, comparison, and future of the representations of additive manufacturing data. Comput Aided Des 2019;111:44–64. [8] Ziegler T, Kraft T. Functionally graded materials with a soft surface for improved indentation resistance: layout and corresponding design principles. Comput Mater Sci 2014;86:88–92. [9] Jackson T, Liu H, Patrikalakis N, Sachs E, Cima M. Modeling and designing functionally graded material components for fabrication with local composition control. Mater Des 1999;20(2–3):63–75. [10] Dokken T, Skytt V, Barrowclough O. Trivariate spline representations for computer aided design and additive manufacturing. Comput Math Appl 2018 in press. doi:10.1016/j.camwa.2018.08.017. [11] Ezair B, Dikovsky D, Elber G. Fabricating functionally graded material objects using trimmed trivariate volumetric representations. In: Proceedings of SMI’2017 fabrication and sculpting event (FASE), Berkeley, CA, USA; 2017. [12] Altenhofen C, Luu TH, Grasser T, Dennstdt M, Mueller-Roemer JS, Weber D, et al. Continuous property gradation for multi-material 3d-printed objects. In: Proceedings of the solid freeform fabrication symposium, 29; 2018. p. 1675–85. [13] Catmull E, Clark J. Recursively generated b-spline surfaces on arbitrary topological meshes. Comput Aided Des 1978;10(6):350–5. [14] Joy KI, MacCracken R. The refinement rules for Catmull-Clark solids; 1999. [15] Altenhofen C, Müller J, Weber D, Stork A, Fellner DW. Direct limit volumes: constant-time limit evaluation for Catmull-Clark solids. In: Fu H, Ghosh A, Kopf J, editors. Pacific graphics short papers. The Eurographics Association; 2018. ISBN 978-3-03868-073-4. doi: 10.2312/pg.20181285. [16] Miyamoto Y, Kaysser W, Rabin B, Kawasaki A, Ford RG. Functionally graded materials: design, processing and applications, 5. Springer Science & Business Media; 2013. [17] Kou X, Tan S. Heterogeneous object modeling: a review. Comput Aided Des 2007;39(4):284–301. [18] Chandru V, Manohar S, Prakash CE. Voxel-based modeling for layered manufacturing. IEEE Comput Graph Appl 1995;15(6):42–7. [19] Aigner M, Heinrich C, Jüttler B, Pilgerstorfer E, Simeon B, Vuong A-V. Swept volume parameterization for isogeometric analysis. In: Proceedings of the IMA international conference on mathematics of surfaces. Springer; 2009. p. 19–44. [20] Martin W, Cohen E. Representation and extraction of volumetric attributes using trivariate splines: a mathematical framework. In: Proceedings of the sixth ACM symposium on Solid modeling and applications. ACM; 2001. p. 234–40. [21] Richards DC, Abram TN, Rennie AEW. Designing digital materials with volumetric gradients. In: Proceedings of the 15th rapid design, prototyping & manufacturing conference (RDPM2017); 2017. [22] Stratasys Ltd.. Grabcad voxel print. 2016. http://www.stratasys.com/industries/ education/research. [23] Autodesk. Monolith. 2019. http://www.monolith.zone/. [24] Massarwi F, Elber G. A b-spline based framework for volumetric object modeling. Comput Aided Des 2016;78:36–47. [25] DeRose T, Kass M, Truong T. Subdivision surfaces in character animation. In: Proceedings of the 25th annual conference on computer graphics and interactive techniques. Citeseer; 1998. p. 85–94. [26] Schaefer S, Hakenberg J, Warren J. Smooth subdivision of tetrahedral meshes. In: Proceedings of the Eurographics/ACM SIGGRAPH symposium on geometry processing. ACM; 2004. p. 147–54. [27] Bajaj C, Schaefer S, Warren J, Xu G. A subdivision scheme for hexahedral meshes. Vis Comput 2002;18(5–6):343–56. doi:10.1007/s003710100150. [28] Chang Y-S, Qin H. A framework for multi-dimensional adaptive subdivision objects. In: Proceedings of the ninth ACM symposium on solid modeling and applications. Eurographics Association; 2004. p. 123–34. [29] Zhu X-P, Hu S-M, Tai C-L, Martin RR. A marching method for computing intersection curves of two subdivision solids. In: Mathematics of Surfaces XI. Springer; 2005. p. 458–71. [30] Severn A, Samavati F. Fast intersections for subdivision surfaces. In: Proceedings of the international conference on computational science and its applications. Springer; 2006. p. 91–100. [31] Stam J. Exact evaluation of Catmull-Clark subdivision surfaces at arbitrary parameter values. In: Proceedings of the 25th annual conference on Computer graphics and interactive techniques. New York, NY, USA: ACM; 1998. p. 395–404. [32] Sunday D. Inclusion of a point in a polygon. 2012. http://geomalgorithms.com/ a03- _inclusion.html. [33] Ma YL, Hewitt WT. Point inversion and projection for NURBS curve and surface: control polygon approach. Comput Aided Geomet Des 2003;20(2):79–99. [34] Jacob B, Guennebaud G. Eigen: a c++ template library for linear algebra. 2008. http://eigen.tuxfamily.org/. [35] Fraunhofer IGD. Cuttlefish. 2019. https://www.cuttlefish.de/. [36] Floyd RW. An adaptive algorithm for spatial gray-scale. In: Proceedings of the society for information display, 17; 1976. p. 75–7. [37] Cho W, Sachs EM, Patrikalakis NM, Troxel DE. A dithering algorithm for local composition control with three-dimensional printing. Comput Aided Des 2003;35(9):851–67. [38] Kremer M, Bommes D, Kobbelt L. Openvolumemesh–a versatile index-based data structure for 3d polytopal complexes. In: Proceedings of the 21st international meshing roundtable. Springer; 2013. p. 531–48. [39] Mueller-Roemer JS, Altenhofen C, Stork A. Ternary sparse matrix representation for volumetric mesh subdivision and processing on GPUS. In: Computer graphics forum, 36. Wiley Online Library; 2017. p. 59–69.