Boundary conforming mesh to T-NURCC surface conversion

Computers & Graphics 82 (2019) 95–105

Contents lists available at ScienceDirect

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

Special Section on SMI 2019

Boundary conforming mesh to T-NURCC surface conversion Martin Marinov∗, Marco Amagliani, Peter Charrot Autodesk, Inc., Autodesk House, Castle Park, Cambridge CB3 0RA, United Kingdom

a r t i c l e

i n f o

Article history: Received 10 March 2019 Revised 20 May 2019 Accepted 21 May 2019 Available online 28 May 2019 Keywords: Mesh conversion T-spline and T-NURCC surfaces Boundary curves and conditions Transfinite interpolation Surface approximation Generative design

a b s t r a c t We present a method for converting an open triangle mesh, whose boundary is constrained to a set of smooth curves, to a T-NURCC surface approximating seamlessly both the mesh and the curve geometry. The curves represent boundaries of B-Rep faces adjacent to the mesh and include intersections, trimcurves, patch boundaries, etc. This conversion problem arises in generative design and topology optimization, where “organic” meshes are synthesized by a physics solver to complement an existing mechanical B-Rep model. It can also occur in reverse engineering, where a scanned mesh is converted to a collection of adjacent analytical and free-form surfaces. The core of our algorithm is a novel transfinite interpolation construction defining a C0 patch network interpolating the curves at the boundaries and coincident with the mesh in the interior. We convert this network to a T-NURCC surface by applying an approximation algorithm designed to produce a very accurate boundary fit and reproduce exactly all boundary corners. These properties allow a B-Rep modeling kernel to stitch the output T-NURCCs with adjacent B-Rep faces into a watertight, low tolerance B-Rep CAD model. We demonstrate our method’s results and applicability on several examples synthesized by topology optimization and generative design applications. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction To meet increased demand for more efficient and lightweight designs that can take advantage of recent advances in additive manufacturing, engineers are increasingly turning to physics simulation-based methods for CAD geometry definition. In particular, generative design and topology optimization [1] offer a quick and cost-effective way to obtain optimized geometry which fulfills the physical design requirements. As most physics solvers operate on volumetric representations such as level sets [2] and tetrahedral meshes [3], the optimized designs are exported in almost all cases as meshes. Hence, to integrate these designs in a CAD assembly, engineers typically have to remodel the synthesized shape from scratch using traditional modeling techniques at a significant time cost. CAD systems can help avoid the remodel step by providing tools to reconstruct a B-Rep model from the optimized design mesh. This model can be built by stitching optimized, organic mesh segments to residual surfaces extracted from the B-Rep input to the optimization process. Since mesh geometry is not sufficiently smooth for downstream CAD modeling operations, the mesh segments are converted to smooth surfaces using established methods such as [4,5]. In general, the boundaries of surfaces output ∗

Corresponding author. E-mail address: [email protected] (M. Marinov).

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

by existing mesh conversion methods will not meet accurately the B-Rep geometry. Hence, tolerant B-Rep edges and vertices [6] are created when the conversion surfaces are stitched to the residual trimmed B-Rep. If these tolerances are not sufficiently small, they can lead to downstream modeling failures. For instance, when offsetting large tolerance B-Reps, adjacent surfaces can pull apart too far from each other, leading to a model that cannot be made watertight. To solve this problem, we present a boundary conforming mesh conversion method that can make the tolerances extremely small: if  > 0 is the underlying B-Rep resolution, i.e., the smallest representable distance in the model, the gap from the conversion to the corresponding B-Rep boundary can be made < 10 (Fig. 1). In downstream modeling operations, such tiny tolerance B-Reps are practically equivalent to precise, non-tolerant models. Our input is an open, general topology mesh M where ∂ M is constrained to a set of smooth (analytical, procedural, spline) curves {Ci }. These curves represent boundary edges on the residual B-Rep surface and can include intersections, trim-curves, patch boundaries, etc. The conversion outputs a T-NURCC [7] surface S, where ∂ S approximates {Ci } to a very tight tolerance ξ > 0 and exactly reproduces all corner points using G1 boundary discontinuities. Therefore, the maximum tolerance from stitching the conversion surface to the residual B-Rep is < ξ . As in prior methods, our algorithm also bounds the approximation error between M and S to a user-specified tolerance δ : δ ≥ ξ > 0. Since

96

M. Marinov, M. Amagliani and P. Charrot / Computers & Graphics 82 (2019) 95–105

Fig. 1. Our conversion applied to an organic mesh (in beige) synthesized by topology optimization of a CAD solid model. Left: the mesh boundary is constrained to the highlighted trim curves (in blue) bounding the residual trimmed solid B-Rep surface (in transparent dark tint). Please note the gaps between the boundary mesh edges and the trim curves in the zoom-in. Middle: we construct a T-NURCC surface seamlessly approximating the trim curves to a very tight tolerance and the mesh interior to a 103 × larger tolerance. The largest boundary deviation is just 4.9 × 10−6 × the mesh bounding sphere diameter. The control T-mesh complexity (5287 patches) gradually decreases away from the boundary. Right: a watertight solid B-Rep stitched from the conversion output and the trimmed solid surfaces. T-spline patches are merged to produce a simpler final B-Rep face layout. The largest edge tolerance in the stitched B-Rep is just 3 , where  = 10−6 is the model resolution (smallest representable distance). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

ξ is dependent on the B-Rep model resolution, while δ is related to the coarse geometry representation of the physics solver that synthesized M, typically δ  ξ . As shown in Fig. 1, our method is able to seamlessly approximate both {Ci } and M, and produce an adaptively structured control T-mesh defining a surface that meets these vastly different boundary and interior approximation requirements.

1.1. Related work Surface reconstruction from mesh data is a well-studied topic underpinning modern reverse engineering applications. Since a full review of this research field is beyond the scope of our work, we refer to the recent surveys [8,9]. Instead, we focus on the most relevant topic of general topology mesh to smooth surface conversion. In [4], a conversion to a B-spline G1 patch network was presented, which bounds the maximum approximation error. In [10,11] hierarchical B-spline methods are proposed to reconstruct complex local surface details. Truncated hierarchical B-splines [12] were used more recently for CAD model reconstruction [13] and surface approximation [14]. At the same time, subdivision surface reconstruction approaches emerged, offering a significantly simpler management of the global surface continuity. Conversion to a globally smooth Catmull-Clark subdivision surface [15] was first proposed in [16]. Hoppe et al. [17] addressed the L∞ (S, M) < δ problem by constructing an adaptively refined piece-wise smooth Loop subdivision surface [18] with feature curves representing G1 surface discontinuities. Subdivision surface conversion techniques were further explored in [19–22]. Since subdivision surface representations lack local refinement methods that preserve the control mesh structure, these methods resort to empirical techniques to insert locally more control points where necessary to solve

L∞ (S, M) < δ . This can result in excessive control mesh complexity and poor output surface parameterization alignment with the input shape. In [5,23], the globally continuous parameterization algorithms [24,25] are applied to construct a T-NURCC surface domain aligned with the input mesh. Local T-spline refinement [26] allows these methods to solve L∞ (S, M) < δ without compromising the surface control structure as in prior techniques. More recent work in this area [27–29] explores improvements on performance and adaptivity. Boundary conditions are not supported by these methods, which limits accuracy when the conversion is integrated within a CAD model. Recognizing the need to represent boundary geometry accurately when embedding free-form surfaces within general models, specialized techniques were developed to conform to boundary conditions. In [30], combined subdivision surfaces were proposed to perform transfinite interpolation [31] of curves assigned to the boundary of a subdivision surface. In [32,33] and [34,35], multi-resolution subdivision surfaces [36] are used to accurately approximate boundaries and sharp features by assigning displacements on multiple levels of an adaptively-refined subdivision surface. These methods do not address the problem of seamless approximation of smooth curves and boundary conforming meshes. Also, these procedural surface representations are not supported in readily available CAD software. In contrast, our algorithm outputs a T-NURCC surface which can be easily converted to a G1 NURBS patch network, a representation universally used in the CAD industry. Another class of related algorithms [37–39] were developed recently to convert a NURBS CAD model with trim curves to a trimfree Catmull-Clark or T-NURCC surface. These methods can take advantage of the structure and accuracy of the input NURBS surfaces and therefore do not face the same challenges as unstructured mesh conversion techniques.

M. Marinov, M. Amagliani and P. Charrot / Computers & Graphics 82 (2019) 95–105

97

Fig. 2. Map size-field effect on our transfinite interpolation. The input dome-shaped mesh edges are constrained to a smooth curve (in blue on the far left). Next, from left to right, we show three different conversions. The target quad edge size is halved for each subsequent example. Correspondingly, the surface area influenced by the curve progressively decreases, but maintains alignment with the curve shape. The error-driven refinement of the surface domain is disabled for these conversions to illustrate the transfinite interpolation construction in isolation from the rest of our method. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

1.2. Contributions Seamless integration of meshes within CAD models is an important, but elusive goal. While some progress in this direction was made recently in [40], in practice, the lack of smoothness severely limits mesh usability in downstream modeling operations. Hence, automatic conversion of meshes to smooth surfaces, which can be integrated in watertight, low tolerance B-Rep models, is a key enabler for generative design and topology optimization CAD workflows. In this work, we present a boundary conforming mesh conversion method specifically designed to meet the requirements of these workflows. Our main contribution is a transfinite interpolation construction that seamlessly combines the input mesh with B-Rep boundary curves. A simple, but flawed solution to this problem is to replace every triangle with boundary edges with a triangular patch [41] which interpolates the corresponding boundary curve segments and interior mesh edges. This leads to a jagged, tessellation-biased area where {Ci } influence the output surface (Fig. 3(a)). Instead, in our method, we compute a global quad map of M and blend the map image with procedural quad patches interpolating {Ci } at the boundaries of the domain. Thus, our construction aligns the surface area affected by the interpolated smooth curves with the boundary shape (Fig. 3(b)). Moreover, the extent of the curve influence can be controlled by a user-specified sizefield independently of the mesh tessellation as shown in Fig. 2. Existing conversion methods use a global tolerance δ > L∞ (S, M) to bound the approximation error. This is impractical to meet the required accuracy in our applications: for the example in Fig. 1, enforcing the required boundary error bound globally required

2.8 × 105 patches. Hence, we introduce a boundary error bound ξ : δ ≥ ξ > 0. Since in practice δ  ξ , we propose techniques that balance the refinement to maintain stability while solving the nonlinear problems L∞ (∂ S, {Ci }) < ξ and L∞ (S, M) < δ . This strategy achieves the required boundary accuracy at acceptable patch density: the conversion surface in Fig. 1 consists of just 5287 patches. We also propose several modifications to the mesh parameterization at the boundary to model the input curves continuity more efficiently. 1.3. Overview Our input is a general topology, orientable 2-manifold mesh M with boundary and C1 curve segments {Ci (t): t ∈ [ti,0 , ti,1 ]} assigned to the boundary edges. {Vi }, {Ei }, {Ti } are the mesh vertices, edges and triangles. ∂ M is constrained to the curve segments, i.e., for each Ei = Vi,0 , Vi,1 ⊂ ∂ M, we have Ci (ti,0 ) = Vi,0 and Ci (ti,1 ) = Vi,1 . The conversion consists of two main stages (Fig. 4). In the first ˜ that seamlessly stage, we construct a C0 quad patch network M combines the mesh M with the B-Rep boundary curves {Ci }. The ˜ is extracted by computing a globally conquad patch layout of M tinuous map φ : M →  ⊂ R2 . Using φ , we map ∂  to the input curves {Ci } and define patches interpolating the displacement field {Ci } − φ −1 (∂ ). Adding these patches to M = φ −1 (), we obtain ˜ () interpolating {Ci } at ∂  and coincident with the C0 surface M M in the interior. In the second stage, we define a T-NURCC surface S(), and ˜ ) < ξ and proceed to solve the approximation problems L∞ (∂ S, ∂ M ˜ ) < δ, where δ ≥ ξ > 0. L∞ (∂ S, ∂ M ˜ ) < ξ is solved first using L∞ (S, M an iteration of least-squares optimization and error-driven adaptive refinement. The refinement also includes a balancing component that propagates new knots to the surface interior where necessary to ensure the overall process stability. Once the boundary conditions are met, we freeze the boundary control points and carry out a similar process in the interior. 2. Transfinite interpolation In this section we develop the parameterization and geometry ˜ which interpolates the input curves {Ci } at ∂ M of the C0 surface M and coincides with M away from the boundary. 2.1. Parameterization

Fig. 3. (a) Interpolating the smooth curve (in blue) using triangular patches replacing the input triangles adjacent to the boundary affects a surface area (in gray) with a jagged, tessellation-dependent shape. (b) In contrast, our transfinite interpolation construction leads to a boundary-shaped surface area (in yellow) influenced by the interpolated curve. The influence area is independent of the mesh tessellation shown with dotted lines. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

Unlike classical algorithms for mesh parameterization [42,43], globally continuous parameterization methods such as [24,25,44,45] do not require general topology meshes to be partitioned into segments homeomorphic to a disc. Instead, a set of cone singularities are specified on the mesh [46], and can be

98

M. Marinov, M. Amagliani and P. Charrot / Computers & Graphics 82 (2019) 95–105

Fig. 4. Overview. The first row illustrates the input data and the transfinite interpolation stage. The second row shows the flow of the approximation stage.

obtained automatically from a cross-field computation [47]. Each triangle is mapped to a single chart and continuity across the chart edges is ensured by adhering to transition function constraints [24]. Quantizing the singularity parameters to integer points enables the extraction of a quad partition of the parameter domain [45,47], in which the singularities correspond to the extraordinary vertices of the quad layout. For a comprehensive review of global quad parameterization techniques, we refer to the survey [48] and the tutorial [49]. As demonstrated in [5,23], global parameterizations are a powerful tool in mesh to T-NURCC surface conversion. We follow these works by applying the more recent Integer Grid Map (IGM) method [45]. This algorithm optimizes a convex mixed-integer quadratic problem, so the map computation is exceptionally robust. Singularity parameter quantization is done efficiently as in [50]. To reduce the map distortion and improve the interior singularity positions, we apply Laplacian smoothing of the computed parameter values that resamples M very densely, introducing negligible changes to the mesh geometry. The IGM φ : M →  ⊂ R2 conforms to pre-specified cross- and size-fields. In our conversion, these fields are provided by the application to control the map orientation and stretch. The cross-field sets a parameter frame for each triangle and, typically, is aligned with principal curvature estimates as in [47]. Since the map image  also defines the domain of the output T-NURCC surface S, we modify the cross-field at the mesh boundary to improve ∂ S. The effect of these tweaks becomes evident after the definition of S, so we describe them later in Section 3.1.3. The size-field adapts the quad partition edge lengths to the local feature sizes and defines the extent of influence area of the boundary curves {Ci } (Fig. 2). The map image  is partitioned in the unit quad layout Qi :  i + [0, 1]2 , i ∈ Z2 , {Qi | Qi = }. { i } and {ϒ i } denote the quad layout vertices and edges respectively. The IGM inverse map φ −1 :  → M plays a critical role in our construction: given Qi and (u, v ) ∈ [0, 1]2 , φ −1 (Qi , u, v ) retrieves the corresponding point on M. By storing for each Qi the set of overlapping triangles i = {T j | dim(φ (T j ) ∩ Qi ) = 2}, φ −1 is evaluated very efficiently by locating T j ∈ i : (i + (u, v )) ∈ φ (T j ).

(E0 , E1 , . . . , EN ) and E0 (tA ) = φ −1 (A ), EN (tB ) = φ −1 (B ), tA ∈ [t0,0 , t0,1 ), tB ∈ (tN,0 , tN,1 ]. Joining the edge curve segments C0 [tA , t0,1 ], C1 [t1,0 , t1,1 ], .., CN [tN,0 , tB ] together yields the geometry of  . To complete the map ϒ →  and avoid dependency on the parameterizations of {Ci }, we parameterize  by its arc-length normalized to [0, 1]. Where numerical integration is required to evaluate the t arc-lengths {si (t ) = t i,1 |Ci (t )|dt }, e.g., spline curves, an accurate i,0

B-spline approximation of si is precomputed [51]. 2.3. Patch construction Let D = {Qi | dim(Qi ∩ ∂ ) = 1} be the set of all quads with at least one boundary edge. For each Qi ∈ D we construct a patch ˜ i blending the input mesh geometry φ −1 (Qi ) with the composM ite curves { j } assigned to the boundary edges of Qi (Fig. 5). Let ∂ Qi = {ϒ0 , ϒ1 , ϒ2 , ϒ3 } be the IGM edges ordered anti-clockwise along the isolines v = 0, u = 1, v = 1, u = 0. Let also φi−1 (u, v ) = φ −1 (Qi , u, v ), (u, v ) ∈ [0, 1]2 . ∀ϒ j , if ϒ j ⊂ ∂ , define the vector field curve μ j : [0, 1] → R3 which represents the displacement from φ −1 (ϒ j ) to  j as follows:

μ0 (t ) = 0 (t ) − φi−1 (t, 0 ), μ1 (t ) = 1 (t ) − φi−1 (1, t ),

μ2 (t ) = 2 (1 − t ) − φi−1 (t, 1 ), μ3 (t ) = 3 (1 − t ) − φi−1 (0, t ).

Otherwise, define μ j = 0 ∈ R3 . Then, construct the displacement field blending patch:

i (u, v )=(1 − v ) · μ0 (u ) + v · μ2 (u ) + u · μ1 (v ) + (1 − u ) · μ3 (v )

2.2. Boundary map Since φ : ∂ M → ∂ , we map ∂  → {Ci } by composing a curve  for each ϒ = A B ⊂ ∂ . Let φ −1 (ϒ ) → E ⊂ ∂ M where E =

Fig. 5. An illustration of the transfinite interpolation patch construction.

M. Marinov, M. Amagliani and P. Charrot / Computers & Graphics 82 (2019) 95–105

99

Fig. 6. Boundary conditions. Solid lines represent the quad domains {Qi } and dashed lines represent 0-channel edges. Red dots represent the boundary control points which control ∂ S. Large red dots represent boundary C1 discontinuities. Green dots represent control points which control the tangent plane along ∂ S. (a) Regular vertex. (b) Regular vertex representing a boundary corner with a C1 discontinuity. Local knot insertion is used to create boundary control points defining both adjacent first derivatives along ∂ S. The T-joints are replicated on the interior 0-channel edges to allow tangent plane control. (c) Convex corner vertex corresponding to an index 1 singularity. (d) Concave corner vertex corresponding to an index -1 singularity. As in (b), local knot insertion is applied to provide control of the first derivatives. However, in this case the added boundary knots are propagated into the surface interior (red dashed lines) to prevent the interior basis function from affecting ∂ S. This configuration is also applied for higher valence boundary vertices. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

and finally, the transfinite interpolation patch:

˜ i (u, v ) = φ −1 (u, v ) + i (u, v ). M i In the most common case of a single boundary edge ∈ ∂ Qi , it ˜ i interpolates the boundary curve. is straightforward to see that M If there are two boundary edges on opposite sides of Qi , i simplifies to a linear blend of the corresponding displacements. To ˜ i interpolates the boundary curves when Qi has two see why M (or more) adjacent boundary edges, consider that the vertex  between these edges is a boundary singularity of φ . Since IGM singularities are placed on input mesh vertices by construction [45,47],  coincides with a vertex V ∈ ∂ M. As ∂ M is constrained to {Ci }, the boundary curves interpolate V. Hence, the corresponding displacements and i vanish at  . In particular, if ∂ Qi ≡ ∂ , the displacements vanish at all four corners and i is a Coons patch [52]. ˜ () from the patches M ˜ i ↔ Qi where M ˜ i ≡ φ −1 if We compose M i ˜ follows directly from the continuity of Qi ∈ / D. The continuity of M

φ −1 and { i }.

3. Approximation The second stage of our algorithm approximates the patch net˜ () with a T-NURCC surface S(). It consists of three steps: work M ˜ ) < ξ , freezing ∂ S and solving constructing S(), solving L∞ (∂ S, ∂ M ˜ ) < δ. L∞ (S, M 3.1. T-NURCC surface We extract the initial quad control mesh for the output surface S from the IGM quad partitioning {Qi } (Section 2.1). Since all IGM singularities are located on input mesh vertices, there is no need for local refinement to accommodate “in-triangle” singularities as in [5]. The initial quad layout of S is fully uniform, hence we set the knot interval value di = 1 for each quad edge ϒ i . As we refine S with local edge and quad splits [26], the knot interval values are updated accordingly. Splitting ϒ i → {ϒ j , ϒ k } creates the new knot intervals d j + dk = di . Correspondingly, quadrisecting Qi → {Qi, j,k } j,k=0,1 splits the quad side knot intervals di,u and di,v , equal to the sum of the edge parameter intervals along the u and v sides. 3.1.1. Correspondence From now on, we denote the quad domain partition of the refined T-NURCC surface with {Qi } and the initial quad partition with {Q˜i }. Initially, each patch S(Q˜i ) trivially corresponds to the

˜ i . When refining Q˜i , we maintain the correspondence from patch M each Q j ⊂ Q˜i using a parameter translation and scale map ηi, j : ˜ i (ni, j (u, v )). Later in Q j (u, v ) → Q˜i (u˜, v˜ ), so we have S(Q j , u, v ) → M ˜ i (u, v ) inthe paper we omit this detail and write S(Q j , u, v ) → M stead. 3.1.2. Boundary conditions To control the approximation at ∂ S we set up Bézier boundary conditions by applying the standard method of adding a 0-channel of quad faces along ∂  (Fig. 6). Where necessary, additional vertices are added along the boundary to interpolate boundary first derivatives and represent accurately input boundary corners. At regular boundary vertices and index 1 singularities (Fig. 6(a) and (c), this process is equivalent to boundary knot insertion in tensorproduct NURBS surfaces [53]. Using local knot refinement, we create G1 discontinuities on regular vertices (Fig. 6(b)) to represent boundary corners away from map singularities. However, concave boundary vertices on negative index singularities are significantly more complex to handle as discussed in Appendix A of [39]. If we use the C1 discontinuity interior propagation approach to represent them (Fig. 6(d)), we must add C1 constraints when solving the interior min L2 problem (see Section 3.3.1). 3.1.3. Cross-field boundary adjustments To use simpler boundary conditions where possible, we apply some adjustments to the input cross-field before computing the IGM as described in Section 2.1. First, we prevent boundary singularities everywhere except at G0 joints of the boundary curves by constraining the cross-field index of all non-corner boundary vertices to 0. Thus ∂ S is C2 away from any corners (example in Fig. 13), which reduces the surface complexity as there are no tangent control points on C1 joints. To avoid C1 interior discontinuities emanating from concave boundary corners, we adapt a similar strategy to Sederberg et al. [37], and relocate negative boundary cross-field singularities to nearby interior vertices. This replaces the corner boundary condition in Fig. 6(d) with the one in Fig. 6(b). Let G be the subset of boundary vertices assigned to negative index singularities. ∀Vi ∈ G ˆ = {V j } as shown in we select corresponding interior vertices G Fig. 7. Then, we swap the cross-field indices for all pairs Vi ∈ G ↔ ˆ and update the cross-field. During the subsequent IGM comVj ∈ G putation, we quantize the parameters ∀Vi ∈ G to force φ (Vi ) ∈ Z2 . Applying local knot insertion to the regular boundary control mesh

100

M. Marinov, M. Amagliani and P. Charrot / Computers & Graphics 82 (2019) 95–105

Fig. 7. Relocation of a negative boundary singularity. Overlaid dotted blue lines show the local quad partitions. Left: a singularity of index −1 (blue dot) is placed on a large interior angle α > π (concave) boundary vertex Vi . Right: the singularity is relocated to an adjacent interior vertex Vj such that the edge connecting Vi , Vj is as close as possible to the bisector of α . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article)

vertices {φ (Vi )} to set up the boundary condition in Fig. 6(b) allows us to represent the corresponding boundary curve corners at G. This technique has a few limitations worth noting: If the input mesh near a vertex Vi ∈ G is too coarse, a suitable relocation vertex may be absent or too far from the bisector of the interior angle α . We can overcome this situation by refining the input mesh locally to introduce more interior vertices. Also, if the accumulated interior angle is too large, e.g., α > 2π , and the singularity index < −1, the relocation may introduce large parametric distortion which can lead to poor surface quality. In such cases, we accept the interior C1 discontinuities instead.

Fig. 8. (a) An interior unit quad adjacent to a boundary 0-channel quad with a bisected boundary edge. (b) We replicate the split on the interior 0-channel edge, bisecting the 0-channel quad. This configuration is balanced. (c) If, during subsequent refinement, we bisect one (or both) smaller 0-channel quad, the configuration becomes unbalanced. We quadrisect the interior quad to rebalance the layout.

refinement when a parametric error bound is not required. An alternative error definition would be to use the Hausdorff dis˜ ) between corresponding boundary loop curves tance H (∂ S, ∂ M ˜ ) ≥ H (∂ S, ∂ M ˜ ), this geometric measure can [55]. Since L∞ (∂ S, ∂ M lead to less refinement. As we would like to explore this option in a future work, our results (Section 4) include the graph of ˜ ) for all examples. H (∂ S, ∂ M 3.2.4. Boundary refinement We bisect all {ϒi | ϒi ∈ ∂  ∧ L∞ ϒ ≥ ξ } by a local knot insertion i

3.2. Boundary iteration The boundary approximation iterates several steps: solving min L2 , estimating L∞ and refining where appropriate. Eventually, ˜ ) < ξ. the process converges when L∞ (∂ S, ∂ M L2

3.2.1. Boundary fit ∂ S consists of C2 NURBS curves except at the corner points. We  2 ˜ , ∂ S) = ˜ discretize L2 (∂ M ∂  (∂ M − ∂ S ) as follows:

˜ , ∂ S) = L2 ( ∂ M

 

ϒ ∈∂ 

≈



dϒ

0

(M˜ ϒ (t ) − Sϒ (t ))2 dt



˜ ϒ (tk ) − Sϒ (tk ))2 . wk ( M

ϒ ∈∂  tk ∈[0,dϒ ]

{(wk , tk )} are the weights and corresponding uniform parameters

of a Newton-Cotes quadrature formula [54] approximating the L2 integral over the interval [0, dϒ ] for each ϒ ∈ ∂ . The order of the formula is based on the continuity of the corresponding composite curve  . We solve a least squares system [53] to find the control ˜ , ∂ S ). point positions along ∂ S at min L2 (∂ M 3.2.2. Boundary corners Let  be a boundary quad vertex representing a C1 discontinuity as shown in Fig. 6(b)–(d). Let V ∈ ∂ M be the corresponding input mesh corner vertex. We constrain ∂ S at  by adding interpolation constraints to the least squares system [53] for the position of V and the first derivatives i (di ) and  j (0 ), where (ϒ i , ϒ j ) are its adjacent quad edges. These constraints force ∂ S to represent accurately the boundary shape at V. 3.2.3. Boundary error For each quad edge ϒ ∈ ∂  we apply non-linear numerical op∞ ˜ timization [54] to compute L∞ ϒ = maxt∈[0,dϒ ] |Mϒ (t ) − Sϒ (t )|. Lϒ is evaluated efficiently by restricting interval subdivision to bisection and sharing dyadic parameter samples with the Newton–Cotes quadrature formulas in Section 3.2.1. ˜ ) = maxϒ ∈∂  L∞ depends on the paHowever, since L∞ (∂ S, ∂ M ϒ ˜ rameterization of ∂ M (Section 2.2), it can lead to unnecessary

at the mid-parameter (Fig. 8(a)). This choice integrates with our L2 and L∞ discretization, as we can reuse the samples evaluated in these computations. To prevent highly non-uniform knot sequences which can lead to surface quality issues, we balance the boundary by bisecting all edges

   ϒ ⊂ ∂  ∧ ∃ ϒ j ⊂ ∂  ∧ ϒi ∩ ϒ j  = ∅ ϒi  i . di > 2d j > 0

The boundary edge splits are replicated on the corresponding 0-channel interior edges, bisecting the incident 0-channel quad (Fig. 8(b)). The additional vertex controls the first derivative across the boundary and can be used by G1 boundary conditions. As the boundary is refined, we may create configurations where smallsupport boundary basis functions overlap with larger interior basis functions. To avoid potential adverse effects of such non-uniform parameterization on the surface quality, we adaptively propagate the boundary refinement into the interior of the surface by quadrisecting any quads {Qi } adjacent to the 0-channel if di > 2dj , where di is the 0-channel quad side knot length, and dj is the knot length of the shortest 0-channel edge on this quad side (Fig. 8(c)). Further refinement into the interior may be required to rebalance the parameterization, which we carry out as detailed in Section 3.3.4. 3.3. Interior iteration From this point on, we fix the control points along ∂ S and restrict all modifications to the interior of S. Again, we iterate solving ˜ , S) < δ. min L2 , estimating L∞ and refining until L∞ (M 3.3.1. Interior L2 fit Let dQ,u and dQ,v be the parameter space intervals along the u ˜ Q (u, v ) − SQ (u, v ). and v directions of a quad Q, and RQ (u, v ) = M We discretize the L2 metric using Newton-Cotes cubature:

˜ , S) = L2 ( M ≈

 

(M˜ − S )2 d p =





 Q∈



Q∈ uk ∈[0,dQ,u ] vl ∈[0,dQ,v ]

dQ,u 0



dQ,v 0

R2Q (u, v )dvdu

wk,l R2Q ((uk, vl )).

M. Marinov, M. Amagliani and P. Charrot / Computers & Graphics 82 (2019) 95–105

101

continuity. These constraints are applied on the blue interior control points shown in Fig. 6(d).

Fig. 9. Quad domain samples grid, nesting and refinement. (a) Samples (blue dots) for the L2 metric discretization are defined on a 5 × 5 grid. Each sample consists ˜ and the S basis functions combination. Edge and vertex samof the position on M ples are reused for each incident quad. (b) When evaluating the L∞ error we double the sample grid resolution in both parameter directions. This “nests” the refinement samples (orange dots) within the existing grid. (c) If we subsequently refine a quad, the finer quad grid samples are immediately available, and only the basis function combinations must be updated to account for the added control points. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

{(uk , vl )} define a uniform sample grid over Q and {wk,l } are the ˜ is C0 , we corresponding cubature weights. Since the interior of M use the trapezoidal rule with a 5 × 5 grid (Fig. 9(a)). The least squares system is composed as in [5], except that the control point positions along ∂ S are fixed and there are no fairing terms added to the energy. 3.3.2. Interior discontinuities If there are any interior C1 discontinuities propagated from the boundary conditions setup in Section 3.1.2, we need to add derivative interpolation constraints across them to achieve C1 interior

3.3.3. Interior error While multi-variate optimization [54] can be used to compute accurately the parametric error L∞ = maxu∈[0,dQ,u ],v∈[0,dQ,v ] | Q RQ (u, v )| the computation cost is unacceptable in most applications. Instead, we refine the uniform sample grid as shown in Fig. 9(b) and compute L∞  maxuk ∈[0,dQ,u ],vl ∈[0,dQ,v ] |RQ (uk , vl )| on Q the finer 9 × 9 grid.

3.3.4. Interior refinement We apply quadrisection (Fig. 9(c)) to refine all {Qi |L∞ ≥ δ} and Q balance the parameterization by quadrisecting the quads





 ∃Q j ∧ Qi ∩ Q j = ∅ Qi  . di,u · di,v > 4 · d j,u · d j,v > 0

Extra care is needed when refining {Qi |Qi ∩ ∂  = ∅} in cases where the required edge splits are already present on the boundary, or when boundary edges are split further. In the latter case, we use the added boundary control points to improve the approximation ˜ , ∂ S ) (Section 3.2.1) after comalong ∂ S by solving again min L2 (∂ M pleting the interior iteration. We observe improved surface quality if we refine patches as required to maintain the “standard” T-spline basis property [7,56]. We apply a greedy approach to this issue that leads occasionally to excessive refinement. To reduce the number of refined patches, we are keen to explore the work in [57].

Fig. 10. A generative design organic mesh, ξ = 7 × 10−6 . Final |{Qi }| = 15441, run time: 8m.

102

M. Marinov, M. Amagliani and P. Charrot / Computers & Graphics 82 (2019) 95–105

Fig. 11. Generative design with 3 organic meshes. Metrics shown for the most complex mesh. ξ = 3 × 10−6 . Final |{Qi }| = 34532, run time: 42m.

4. Results Figs. 10, 11, 15, and 16 demonstrate our method on organic meshes produced by topology optimization and generative design applications. All mesh boundaries are constrained to procedural trim-curves constructed to lie exactly on the existing non-organic surfaces in the optimized B-Rep model. In all examples ξ is set to 100 × the B-Rep resolution. To make the tolerances meaningful in the context of the each conversion example, the value of ξ is given relative to the bounding sphere diameter of the input mesh. δ = 100ξ in all examples. Run times include the IGM computation. The synthetic example in Fig. 13 illustrates the boundary iteration in isolation. All results include a detailed chart of our approximation metrics at the iteration indicated on the horizontal axis. The error metric values are shown on the left vertical axis. Since we measure the ˜ , S) boundary metrics only for the boundary iterations and L∞ (M only for the interior iterations, these graphs may not cover the entire horizontal range. The black graph tracks the number of output patches |{Qi }| after each iteration and the values are shown on the right vertical axis.

ponent. For the example in Fig. 1, we refined the input mesh boundary as suggested above, increasing by 10% the original triangle number (2371). To meet the same error tolerances, the conversion produced a surface which has 80% more patches. Close examination of the outcome (Fig. 12) reveals that excess refinement occurs at the surface boundary as the approximation process compensates for the lack of smoothness of the input boundary geometry. Subsequently, the interior patch layout is denser as well due to balancing and standard basis refinement. At the same time, stitching the converted surface to the residual B-Rep resulted in 50 × larger tolerances. The interior approximation iteration in Section 3.3 is a solution for the general mesh to T-NURCC conversion problem. Several techniques differentiate this aspect of our method from the state of the art solution in [5]: the mesh parameterization method is IGM [45] instead of PGP [24], new control points are added at

4.1. Comparison to prior art Our method is designed to operate on boundary constrained meshes and approximate very tightly the smooth boundary curve constraints. To our knowledge, there is no prior art that provides such capabilities. However, a comparable approach can be naively implemented by refining the input mesh boundary (and its vicinity), while displacing the new boundary vertices onto the curve constraints. The refined mesh can be converted using a prior algorithm, e.g., [5]. As discussed in Section 1.2, this would be impractical: even moderately complex models require massive patch density to meet the required boundary accuracy in the interior. To render the approach practical, it must provide a separate boundary error tolerance. We can emulate the potential outcome of such alternatives using our implementation without the transfinite interpolation com-

Fig. 12. Transfinite interpolation benefit. The conversion on the left is output by our method as described in the paper. The conversion on the right is produced without the transfinite interpolation component. Even though we doubled the input mesh density along the boundary, the output surface is 80% denser and the maximum tolerance after stitching is 50 × larger.

M. Marinov, M. Amagliani and P. Charrot / Computers & Graphics 82 (2019) 95–105

103

Fig. 13. First row: an input mesh tessellating the unit quad is assigned a C2 curve (in blue) with a single convex corner. ξ = 10−4 . Approximation metrics: H < ξ at iteration 5, |{Qi }| = 710; 3 more iterations are needed for L∞ < ξ , |{Qi }| = 1255. Second and third row: S at each iteration. Our cross-field singularity constraints ensure that ∂ S is C2 everywhere except at the designated corner point. The corner shape is accurately matched from the very first iteration. The first 3 iterations cause notable refinement in the interior due to balancing and standard T-spline basis refinement. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

the same time for all patches violating the error bound, and no fairing terms are added to the L2 energy in Section 3.3.1. These choices complement well the boundary conforming components, and as we demonstrate in our results, do not compromise the output surface quality in the interior. To support further this observation, in Fig. 14 we include our conversion of a closed mesh processed in [5]. It is worth noting that our interior approximation is significantly more costly per iteration due to the global min L2 recomputation after refinement. This recomputation can be re-

placed by a faster refit technique [5,29,58] if quicker run times are required. 5. Conclusion Our method converges and performs robustly even for the most detailed and complex organic meshes synthesized by the topology optimization and generative design applications that employ it. Since the conversion surface boundaries conform to the boundary edges of the trimmed residual B-Reps, watertight CAD models with very small tolerances can be produced to meet the requirements of downstream modeling operations. On the other hand, the rapid increase of the patch density in the first few iterations to maintain the standard basis and balanced parameterization properties is unwelcome but necessary. Attempts to forgo these properties led to surface quality issues, which perhaps could be overcome by a more careful implementation, that is resistant to the numerical issues arising from rational basis functions and highly irregular spans. Declaration of interests

Fig. 14. The rocker arm mesh converted by our method at the same relative δ = 0.002 as in [5]. The surface has 4185 control points and 3993 patches. Approximation time: 59s; map time: 4m.

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.

104

M. Marinov, M. Amagliani and P. Charrot / Computers & Graphics 82 (2019) 95–105

Fig. 15. An organic mesh from topology optimization, ξ = 5.1 × 10−6 . Final |{Qi }| = 46173, run time: 79m.

Fig. 16. A topology optimized model with 3 organic meshes. Metrics shown for the most complex mesh. ξ = 4.75 × 10−5 . Final |{Qi }| = 10752, run time: 6m.

M. Marinov, M. Amagliani and P. Charrot / Computers & Graphics 82 (2019) 95–105

Acknowledgments We would like to thank the anonymous reviewers for their helpful feedback, Max Lyon, David Bommes and Leif Kobbelt for their support with integrating IGM, Siavash Meshkat, Nanda Santhanam, Karl Willis and Michael Smell for providing generative design and topology optimization test data, Jean Flower, Iain Henley and Tristan Barback for developing our input models, Dominic Prior for proof-reading the paper, Desislava Aleksandrova for helping with the composition of the results images, and our colleagues in the Autodesk modeling components and simulation groups for supporting the development of our algorithms. This work is funded by Autodesk, Inc. References [1] Brackett D, Ashcroft I, Hague R. Topology optimization for additive manufacturing. In: Proceedings of the SFF 11; 2011. p. 348–62. 1, S [2] van Dijk NP, Maute K, Langelaar M, van Keulen F. Level-set methods for structural topology optimization: a review. Struct Multidiscip Optim 2013;48(3):437–72. [3] Chen X, Peng D, Gao S. SVM-based topological optimization of tetrahedral meshes. In: Proceedings of the IMR 13; 2013. p. 211–24. [4] Eck M, Hoppe H. Automatic reconstruction of B-spline surfaces of arbitrary topological type. In: Proceedings of the ACM SIGGRAPH 96; 1996. p. 325–34. ISBN 0897917464. [5] Li W-C, Ray N, Lévy B. Automatic and interactive mesh to T-spline conversion. In: Proceedings of the SGP 06; 2006. p. 191–200. ISBN 3-905673-36-3. [6] Jackson DJ. Boundary representation modelling with local tolerances. In: Proceedings of the ACM SMA 95; 1995. p. 247–54. ISBN 0897916727. [7] Sederberg TW, Zheng J, Bakenov A, Nasri A. T-splines and T-NURCCs. In: Proceedings of the ACM SIGGRAPH 03; 2003. p. 477–84. ISBN 1-58113-709-5. [8] Geng Z, Bidanda B. Review of reverse engineering systems - current state of the art. Virtual Phys Prototyp 2017;12(2):161–72. [9] Buonamici F, Carfagni M, Furferi R, Governi L, Lapini A, Volpe Y. Reverse engineering modeling methods and tools: a survey. Comput-Aided Des Appl 2018;15(3):443–64. [10] Forsey DR, Bartels RH. Surface fitting with hierarchical splines. ACM Trans Gr 1995;14(2):134–61. [11] Forsey DR, Bartels RH. Hierarchical B-spline refinement. SIGGRAPH Comput Gr 1988;22(4):205–12. [12] Giannelli C, Jüttler B, Speleers H. THB-splines: The truncated basis for hierarchical splines. Comput Aided Geom Des 2012;29(7):485–98. [13] Kiss G, Giannelli C, Zore U, Jüttler B, Großmann D, Barner J. Adaptive CAD model (re-)construction with THB-splines. Gr Models 2014;76(5):273–88. [14] Giannelli C, Jüttler B, Kleiss SK, Mantzaflaris A, Simeon B, Špeh J. THB-splines: an effective mathematical technology for adaptive refinement in geometric design and isogeometric analysis. Comput Methods Appl Mech Eng 2016;299:337–65. [15] Catmull E, Clark J. Recursively generated B-spline surfaces on arbitrary topological meshes. Comput Aided Des 1978;10(6):350–5. [16] Halstead M, Kass M, DeRose T. Efficient, fair interpolation using Catmull–Clark surfaces. In: Proceedings of the SIGGRAPH 93; 1993. p. 35–44. ISBN 0-89791-601-8. [17] Hoppe H, DeRose T, Duchamp T, Halstead M, Jin H, McDonald J, et al. Piecewise smooth surface reconstruction. In: Proceedings of the ACM SIGGRAPH 94; 1994. p. 295–302. ISBN 0897916670. [18] Loop C. Smooth subdivision surfaces based on triangles. (Master’s thesis); 1987. [19] Krishnamurthy V, Levoy M. Fitting smooth surfaces to dense polygon meshes. In: Proceedings of the ACM SIGGRAPH 96; 1996. p. 313–24. ISBN 0897917464. [20] Ma W, Zhao N. Catmull–Clark surface fitting for reverse engineering. In: Proceedings of the GMP 00; 2000. p. 274–84. [21] Ma W, Ma X, Tso S-K, Pan Z. A direct approach for subdivision surface fitting from a dense triangle mesh. Comput-Aided Des 2004;36(6):525–36. [22] Marinov M, Kobbelt L. Optimization methods for scattered data approximation with subdivision surfaces. Gr Models 2005;67(5):452–73. [23] Wang Y, Zheng J. Adaptive T-spline surface approximation of triangular meshes. In: Proceedings of the ICICS’07; 2007. p. 1–5. ISBN 1424409837.

105

[24] Ray N, Li WC, Lévy B, Sheffer A, Alliez P. Periodic global parameterization. ACM Trans Gr 2006;25(4):1460–85. [25] Gu X, Yau S-T. Global conformal surface parameterization. In: Proceedings of the SGP 03; 2003. p. 127–37. ISBN 1-58113-687-0. [26] Sederberg TW, Cardon DL, Finnigan GT, North NS, Zheng J, Lyche T. T-spline simplification and local refinement. In: Proceedings of the ACM SIGGRAPH 04; 2004. p. 276. [27] Peng X, Tang Y. Adaptive reconstruction of T-Spline surfaces of arbitrary topological type. In: Proceedings of the CAIDCD 08; 2008. p. 599–603. ISBN 978-1-4244-3290-5. [28] Wang Y, Zheng J. Curvature-guided adaptive T-spline surface fitting. Comput-Aided Des 2013;45(8–9):1095–107. [29] Feng C, Taguchi Y. FasTFit: A fast T-spline fitting algorithm. Comput-Aided Des 2017;92:11–21. [30] Levin A. Combined subdivision schemes for the design of surfaces satisfying boundary conditions. Comput-Aided Geom Des 1999;16(5):345–54. [31] Sabin M. Transfinite surface interpolation. In: Proceedings of the mathematics of surfaces; 1996. p. 517–34. ISBN 0-19-851198-1. [32] Biermann H, Kristjansson D, Zorin D. Approximate Boolean operations on free– form solids. In: Proceedings of the SIGGRAPH 01; 2001. p. 185–94. [33] Biermann H, Martin IM, Zorin D, Bernardini F. Sharp features on multiresolution subdivision surfaces. Gr Models 2002;64(2):61–77. [34] Litke N, Levin A, Schröder P. Fitting subdivision surfaces. In: IEEE Visualization 20 01; 20 01a. p. 319–24. [35] Litke N, Levin A, Schröder P. Trimming for subdivision surfaces. Comput-Aided Geom Des 2001b;18(5):463–81. [36] Zorin D, Schröder P, Sweldens W. Interactive multiresolution mesh editing. In: Proceedings of thef ACM SIGGRAPH 97; 1997. p. 259–68. [37] Sederberg TW, Finnigan GT, Li X, Lin H, Ipson H. Watertight trimmed NURBS. In: Proceedings of the ACM SIGGRAPH 08; 2008. p. 1. ISBN 9781450301121. 27 [38] Shen J, Kosinka J, Sabin MA, Dodgson NA. Conversion of trimmed NURBS surfaces to Catmull–Clark subdivision surfaces. Comput Aided Geom Des 2014;31(7–8):486–98. [39] Shen J, Kosinka J, Sabin M, Dodgson N. Converting a CAD model into a non-uniform subdivision surface. Comput Aided Geom Des 2016;48(C):17–35. [40] Lai LM, Yuen MM. Blending of mesh objects to parametric surface. Comput Gr 2015;46(C):283–93. [41] Gregory JA, Charrot P. A C1 triangular interpolation patch for computer-aided geometric design. Comput Gr Image Process 1980;13(1):80–7. [42] Floater MS, Hormann K. Surface parameterization: a tutorial and survey. In: Dodgson NA, Floater MS, Sabin MA, editors. Advances in Multiresolution for Geometric Modelling. Berlin, Heidelberg: Springer; 2005. p. 157–86. [43] Sheffer A, Praun E, Rose K. Mesh parameterization methods and their applications. Found Trends Comput Gr Vis 2006;2(2):105–71. [44] Khodakovsky A, Litke N, Schröder P. Globally smooth parameterizations with low distortion. ACM Trans Gr 2003;22(3):350–7. [45] Bommes D, Campen M, Ebke H-C, Alliez P, Kobbelt L. Integer-grid maps for reliable quad meshing. ACM Trans Gr 2013a;32(4):1. [46] Kharevych L, Springborn B, Schröder P. Discrete conformal mappings via circle patterns. ACM Trans Gr 2006;25(2):412–38. [47] Bommes D, Zimmer H, Kobbelt L. Mixed-integer quadrangulation. In: Proceedings of the ACM SIGGRAPH 09; 2009. p. 77:1–77:10. ISBN 978-1-60558-726-4. [48] Bommes D, Lvy B, Pietroni N, Puppo E, Silva C, Tarini M, et al. Quad-mesh generation and processing: a survey. Comput Gr Forum 2013b;32(6):51–76. [49] Campen M. Partitioning surfaces into quad patches. In: Bousseau A, Gutierrez D, editors. EG 2017 – Tutorials; 2017. [50] Campen M, Bommes D, Kobbelt L. Quantized global parametrization. ACM Trans Gr 2015;34(6):1–12. [51] Peterson JW. Arc length parameterization of spline curves. Comput-Aided Des 2006;14(2):1–11. [52] Coons SA. Surfaces for computer-adided design of space forms. Technical Report; 1967. [53] Prautzsch H, Boehm W, Paluszny M. Bézier and B-spline techniques. Mathematics and Visualization. Berlin, Heidelberg: Springer; 2002. [54] Press WH, Vetterling WT, Teukolsky SA, Flannery BP. Numerical recipes in C++ : the art of scientific computing. Cambridge University Press; 2002. [55] Chen X-D, Ma W, Xu G, Paul J-C. Computing the Hausdorff distance between two B-spline curves. Comput-Aided Des 2010;42(12):1197–206. [56] Wang W, Zhang Y, Scott MA, Hughes TJR. Converting an unstructured quadrilateral mesh to a standard t-spline surface. Comput Mech 2011;48(4):477–98. [57] Scott MA, Li X, Sederberg TW, Hughes TJ. Local refinement of analysis-suitable T-splines. Comput Methods Appl Mech Eng 2012;213–216:206–22. [58] Lin H, Zhang Z. An efficient method for fitting large data sets using T-splines. SIAM J Sci Comput 2013;35(6):A3052–A3068.