1. Introduction
Hexahedral and hexdominant volumetric meshing of 3D shapes is a well investigated yet still open research topic. At their core, hexahedral meshing algorithms balance fidelity to the input surface geometry against element quality, and seek to compute meshes with well shaped, or as cuboid as possible, elements whose outer surface closely aligns with that of the input model (Section 2). To achieve high surface fidelity and to keep the element budget low, users prefer meshes whose external quads align with the surface curvature directions. Users also often strongly prefer meshes whose edges interpolate geometric or semantic surface feature curves. We propose Loopy Cuts, a robust new meshing algorithm that is specifically designed to satisfy these preferences.
Current attempts at automatic hex mesh generation fail to provide the combination of robustness, feature and curvature alignment we seek (Section 2). Methods based on adaptive grids [Marechal2009] or PolyCube maps [Livesu:2013:PolyCut] do not effectively handle surface features that do not align with the global frame, and cannot easily support alignment to a curvature field (Figure 2). Methods that align mesh elements with a volumetric frame field generate meshes that support feature and field alignment [liu2018singularity; Li:2012] but lack robustness, failing to automatically mesh more complex inputs (Section 6, Figure 3).
Our method complements these existing approaches by providing the desired combination of alignment and robustness. It achieves this goal by employing a meshing strategy that mimics semimanual block decomposition. While this popular semimanual approach is commonly used in industry to create quality hex meshes, existing attempts to replicate the process automatically have failed to produce similar quality results. We successfully automate this process, producing results onpar with those obtained using manual decomposition on a wide range of inputs. Our method produces allhex meshes on the vast majority of inputs (Figure 13), and introduces highly localized nonhex elements when necessitated by the input feature curve networks (Figure 17). As an additional benefit, the user can be optionally involved in the meshing process, customizing the decomposition via a simple user interaction (Figure 18).
Our method is based on two simple observations. We first note that we can obtain a decomposition that produces wellshaped hexahedral elements by forming a set of cutting surfaces bounded by cutting loops distributed strategically across the input model’s surface. We further note that by selecting a set of cutting loops that are aligned with a surface crossfield and which interpolate the input features, we can produce a decomposition that respects the geometric characteristic of the input model.
Following this observation, our meshing method computes cutting loops and cutting surfaces bounded by them that satisfy the combination of the criteria above. Starting from a crossfield which is aligned to a set of semantic or geometric feature curves (Figure 1a), we extract a set of welldistributed fieldaligned loops on the object surface (Figure 1b). We facilitate formation of simple low curvature loops, and corresponding wellshaped cutting surfaces, by balancing field alignment against loop geodesicity. We generate and dynamically replenish a pool of candidate cutting loops. We repeatedly select the most impactful cutting loops from this pool, ones farthest from previously selected ones. We form cutting surfaces that interpolate these loops and adapt to the shape of the input model by using level sets of a smooth volumetric field; we simplify cutting surface geometry by enabling formation of cutting surfaces bounded by multiple loops. We use these surfaces to decompose the volume into coarse simple polyhedral blocks, terminating the process once the blocks satisfy our quality requirements (Figure 1c). We refine the resulting blocks via midpoint subdivision producing a coarse hexdominant mesh that can then be used to produce high quality hexahedral meshes across a range of resolutions (Figure 1d).
We validate our framework by testing it on multiple inputs of varying complexity, exhibiting a range of geometric and user prescribed features (Section 6). We were able to generate featurepreserving fieldaligned meshes across all tested models, and to obtain allhex meshes on the great majority of them. The minimal element quality of the meshes we produced across all models was 0.33, comfortably above the minimal threshold of 0.2 typically expected by industrial users [Pebay]. We highlight the advantages of our method by comparing our results to those produced using a range of existing strategies.
Our overall contribution is a robust and fully automatic method for block decomposition based hexahedral (or hexdominant) meshing. We generate comparable or better quality meshes than previous approaches, while guaranteeing featurepreservation and surfacefield alignment. This contribution is made possible by our use of fieldaware surface loop computation, generation of valid, smooth cutting surfaces that interpolate these loops, and the loop evaluation procedure that leads to robust detection of most effective loops to embed into the block structure.
2. Background and Related Works
Hexahedral meshing
The generation of highquality hexahedral meshes is a well researched open problem, which is regarded as extremely difficult among practitioners [Shepherd2008; blacker2000meeting; Owen]. Quality meshes are expected to satisfy two core requirements: mesh elements have to be maximally close to cuboids, and the outer surface of the mesh has to closely approximate the input model. Additional requirements include [blacker2000meeting]: geometric matching, which requires meshes to include all input surface feature curves in the mesh edge network; boundary sensitivity, which implies strong preference for meshes aligned with surface curvature directions; and orientation insensitivity, which requires the mesh to be independent of the model orientation. Unfortunately, none of the existing methods can robustly satisfy the union of the five requirements above on general inputs.
Traditional approaches to hexahedral meshing can be roughly characterized into a few categories [Shepherd2008; blacker2000meeting]: primitivebased, indirect, advancingfront, gridbased, skeletonbased, and decomposition. Primitive based methods and indirect methods are only capable of producing acceptable quality meshes on narrow sets of geometries [Shepherd2008].
Advancing front methods [tautges1996whisker; Plastering1; Kremer2014]
seek to propagate an existing quad surface mesh inward. Recent research on fieldbased hexahedral meshing can be seen as an extension of this approach: the methods first propagate a tensor field from the surface into the interior of the model and then use this volumetric field to compute a mesh
[Nieser2014; Li:2012; jiang2014frame; kowalski2016smoothness; liu2018singularity; solomon2017boundary]. A core advantage of both recent and traditional advancing front approaches is the ability to support all three of the additional requirements above. Unfortunately, both types of approaches do not robustly generalize to models with complex topology. As noted by Liu et al. [liu2018singularity] existing volumetric field methods require manual intervention to obtain desirable results, even for relatively simple shapes (see Figure 13 in their paper). Similar to these fieldbased methods, we use a surface crossfield as a starting point, but instead of propagating it into the model’s interior, we use it to guide the formation of cutting loops and surfaces, sidestepping the challenge of processing interior singularities, and achieving all the three goals above.Grid [Lin:2015:QGA:2828638.2828677; Schneiders1996] and octreebased methods [Marechal2009; ito2009octree] use a regular or adaptive grid to mesh the interior of the input models, and use different strategies to connect this mesh to the model’s surface. While these methods can robustly handle the core hexmeshing requirements, the meshes they produce are orientation dependent and are not aligned with curvature directions. Capturing surface features using these approaches requires excessive local refinement (Figure 2 middle), and can result in feature loss.
Polycubebased meshing methods map generic 3D shapes to orthogonal polyhedra (or polycubes [tarini2004polycube]), mesh this polycube using a regular grid, and then map the resulting mesh back to the original shape [Gregson:11; fu2016efficient; huang2014; Livesu:2013:PolyCut; Fang:2016]. These methods produce gridconnectivity meshes in the interior of the model and can be seen as a generalization of the gridbased approach. Similar to gridbased approaches, these methods are orientation dependent and may not match surface features, even after mesh refinement (Figure 2 left).
Automated decomposition techniques seek to cut the object into parts, which can then be meshed conformingly using existing algorithms. Insideout skeleton [LAPS17; Livesu2016] and medialaxis based decomposition approaches [li1995hexahedral; Sheffer1999; Quadros] fail to generalize to complex shapes. Methods that start from a dense hexmesh and derive a coarse block decomposition from it [gao2017simplification; Cherchi; gao2015hexahedral] may fail to align with features that were not present in the input mesh. Surfacedriven blockdecomposition techniques use surface features to guide cutting surfaces. Cuts can be used to define both the primal [cooper; BlackerMulti; RuizGironesMultisweep; Liu97automatichexahedral] or the dual structure [GAO2018]. Dual methods hardly support feature alignment, because cuts do not directly produce mesh edges, and geometric snapping should be used in post processing to conform with the input features. Primal methods are orientation independent, and frequently satisfy both geometry matching and boundary sensitivity. Unfortunately, both the established [Shepherd2008] and more recent [WANG2017; kowalski2012fun] automatic blockdecomposition methods are limited in the set of geometries they can be applied to.
They thus are highly reliant on the sharp feature networks on the model surfaces, and cannot process freeform natural shapes, or shapes with smooth and rounded features. Our method follows the blockdecomposition approach popularized by these techniques, and inherits their core advantages. However, contrary to those methods, it does not rely just on the feature curve networks on the input models, and can robustly handle generic freeform models with both smooth and sharp features (Figure 3). In particular, the torus and the wave in the inset are examples which [GAO2018] and [WANG2017], respectively, list as failure cases.
Hexdominant meshes.
A range of methods [LevyLP; hybridHexa; ray2017hexahedral; Sokolov:2016:HM:2965650.2930662] support computation of mixed element meshes aiming to minimizing the percentage of nonhex elements while satisfying mesh quality requirements. While recent hexdominant methods are very robust in terms of the inputs they handle, they still produce a significant percentage of non cuboidal elements [hybridHexa; ray2017hexahedral]. Our approach produces nonhex elements only incident to singular (non valence3) surface feature vertices (Section 6).
Crossfields and fieldcoherent loops.
Generating and tracing direction fields on surfaces, or other spatial domains, is becoming a fundamental preprocessing step in numerous applications in computer graphics and geometry processing [vaxman2016directional; de2016vector]. Paths traced using most existing methods are not designed to be closed, and are typically terminated when approaching a singularity or another similarly directed path. A range of recent methods seek to connect crossfield singularities with short, field aligned paths [boiermartin2004ptm; carr2006rectangular; daniels2009semi]. We trace closed fieldcoherent loops away from singularities following the approaches of [campen2012dual; Pietroni:2016], which both rely on the formalism introduced in [kalberer2007qsp]. We use the discrete graph based structure of [Pietroni:2016] to efficiently trace such loops and compute fieldaware geodesic distances.
3. Overview
The input to our method is a 3D model described by a closed triangle mesh , together with a set of line features demarcated as chains of edges in . We compute a crossfield on the surface using stateoftheart methodology [bommes2009mixed; ComplexRoots:Diamanti:2014]. We constrain the field to to follow the line features in the input and the main curvature directions elsewhere. We incrementally decompose the model into blocks using a sequence of cuts, each of which cuts through the entire model, producing conforming blocks with shared surfaces. Our choice of cutting surfaces is guaranteed to produce blocks whose boundary edges incorporate the input feature curves and are aligned with the input crossfield. The assembly of such blocks, namely their combinatorial structure, will be referred as a metamesh ; this structure is updated after each cut and it is formed by cells that correspond to the individual blocks.
Once our incremental cutting results in blocks that are simple enough, we can consider an actual geometrical representation of the metamesh composed from polyhedral cells whose faces correspond to portions of the original surface or of the cutting surfaces and whose edges and vertices correspond to the intersections between these different surfaces (note that cell faces need not be planar).
3.1. Decomposition Goals
To obtain the desired output hex (or hexdominant) mesh quality we aim for the metamesh to satisfy the following topological and geometrical criteria. We strictly enforce the Boolean criteria when possible and optimize the continuous ones.
 t1.:

Each cell must have genus zero.
 t2.:

Each cell face must be bounded by at least three edges.
 t3.:

Each vertex of each cell should have valence three.
 g1.:

Each cell should be convex.
 g2.:

Each cell should be well shaped, i.e. have planar faces and orthogonal edges.
 g3.:

Each cell should approximate its corresponding block within a given accuracy.
Criteria t1 and t2 guarantee that subdividing each cell would produce a polyhedral mesh. Criterion t3 guarantees that each cell can be split into hexahedra with a single step of midpoint subdivision. However strictly satisfying it while still incorporating all surface line features may in practice be impossible – see, for instance, the Schneider’s four sided pyramid [Shepherd2008] (inset) – when the combination of sharp and user drawn feature curves results in an input for which no practical hexmesh exists. We therefore require our method to produce cell vertices of valence three away from irregular featurecurve network vertices, but place no constraints on the valence of cell vertices placed at irregular network vertices. Consequently, absent such vertices our method is guaranteed to produce an allhex mesh. If such vertices are present the nonhexahedral elements are constrained to the cells containing them.
Criteria g1 and g2 control the quality of the metamesh cells and indirectly determine the shape of the elements in the final hex mesh. As noted by Owen [Owen] a single nonconvex element in a hexahedral mesh makes a simulation using this mesh suspect. Finally, criterion g3 ensures that the meshes produced by refining the meta mesh and projecting the new vertices to the corresponding block surfaces produces meshes with comparable quality to that of the cells. We address such criteria indirectly, as explained in the following subsection.
3.2. Algorithm
Informally speaking, we incrementally slice the model by applying cuts, which are defined by loops that strictly adhere to the input line features and follow the input cross field. We update the metamesh during this incremental process, and stop the process once the decomposition satisfies the above goals. Once this decomposition is complete we refine it using midpoint subdivision, converting all cells with valencethree vertices into hexahedra. We then refine the mesh to the user designed density and optimize the mesh quality, using offthe shelf optimization code [Livesu:2015].
Generating ordered loops. Section 4
In order to obtain an even, fieldaware, distribution of the cuts that respect the input line features, we strategically specify the order in which we choose the cutdefining, or cutting, loops. For this purpose we assemble a preliminary ordered queue of loops that we later scan to generate the cuts that define our decomposition.
We determine the order as follows: first, we insert the loops that are needed to interpolate the input features (Sec. 4.1); we then add loops that, respecting the crossfield, progressively cover the surface in a uniform way (Sec. 4.2). We achieve this increasingly denser and denser uniform distribution by defining a distance over the space of loops, and applying a biased furthestfirst sampling procedure using this metric. We bias this sampling procedure to satisfy criterion t2 by making sure that each loop in the queue is intersected by other loops at three or more points (see inset: the red loops have only two intersections; adding one more loop takes them all to three intersections, giving a configuration sufficient to produce a valid metamesh). Thanks to the alignment with the cross field, the selected loops intersect (near)orthogonally, satisfying criterion g2.
Processing the loop queue, Section 5
We build our collection of cutting surfaces by repeatedly extracting the top loop from the ordered queue and using it as a starting point for construction of the cutting surface.
We note that in many instances, such as when forming genus zero blocks starting from a higher genus surface, or when processing models with deep concavities, it may be necessary or simply better to form cutting surfaces bounded by more than one loop (see Figure 7). We therefore develop a loop grouping technique that, given a loop extracted from the queue, finds additional loops such that the combined set jointly bounds a well formed cutting surface. The resulting set of one or more loops partitions the input surface into two disjoint charts. We compute a harmonic field in the interior of the model whose isosurfaces smoothly interpolate these charts. We use the isosurface approximately equidistant from these two charts as our cutting surface (Sec. 5.2).
Once a surface is computed, we use the intersections between it and the current set of blocks to refine the block structure and the metamesh by defining new cells, faces, edges and vertices. To satisfy criterion t3, the cutting process constrains new metamesh vertices to lie at the intersections of three surfaces, ensuring they have valence three in all their incident cells. Thus all cells away from singular featurenetwork vertices in our metamesh are split into hexahedra during subdivision. Moreover, since surfaces intersect nearly orthogonally, we are able to warrant that criterion g1 will be satisfied.
We continue to cut the model using new loops from the queue until all cells satisfy our topological requirements and provide a good approximation of the outer shape (criterion g3). The quality of approximation is evaluated simply by measuring the difference of area between each external face of the metamesh and its corresponding patch of surface on mesh .







(a)  (b)  (c)  (d)  (e) 
4. Computing Cutting Loops
Given a cross field over a surface , we trace fieldcoherent geodesic paths w.r.t. , as in [Pietroni:2016]. A formal definition of fieldcoherent paths and loops rests upon the stratification of as defined in [kalberer2007qsp] and briefly summarized in Appendix A.1. Informally we define a fieldcoherent geodesic loop through a point to be a closed curve that goes through , is loosely following one of the directions of and is as short as possible according to the anisotropic distance defined in Equation 3 (Appendix A.1). Roughly speaking, fieldcoherency forces a loop to approximately follow the underlying direction field until it gets back to its origin. Requiring a curve to follow exactly is unlikely to produce closed loops, as tracing the filed directions exactly would usually result in a spiral effect; allowing the curve to drift from while remaining as short as possible allows us to use as a guide to extract consistent loops.
Each point on which is away from a field singularity can be crossed by exactly two fieldcoherent geodesic loops, which are orthogonal to one another at (disregarding the traced curve’s orientation). We consider the set representing the space of all loops for any point .
We incrementally build a queue of loops , which will be scanned to generate the cuts, as follows:

We initialize with loops that incorporate a subset of the input curve features (Sec. 4.1);

We extend using a sampling process that aims to select loops that intersect orthogonally and are evenly distributed on . This step is based on furthest sampling on and requires the definition of a suitable distance on the space of loops (Sec. 4.2).
Since the space of loops is infinite, the sampling process is discretized by forming an extensive pool of loops from which we sample, and which is built dynamically during the selection process (Section 4.3).
4.1. Generating Loops Incorporating Line Features
The extraction and classification of sharp input features is out of the scope of this paper and they are considered part of the input. Our initial feature set consists of both open and closed curves. We form a feature curve network by placing vertices at feature curve intersections, sharp corners, and dangling endvertices of open curves. We break the curves into segments bounded by vertices and treat each segment as a separate curve in the process below. We classify each feature curve as
concave, convex or flat based on the average dihedral angle along it. We employ this classification to determine the number of cutting loops we want to incorporate each feature curve into: we expect each flat feature to belong to a single cutting loop and expect each concave feature to belong to two cutting loops. We do not perform cuts along the convex features, but use them to initialize the metamesh (Section 4.4). This cutting loop formation strategy fosters the formation of asrightas possible dihedral angles for all blocks incident to a given feature (Figure 9).We initialize the cutting loop queue as follows:

We add to all loops formed by closed nonconvex features.

We extend each open concave feature into two complete loops by loop tracing and extend each open flat feature into one loop (Appendix A.1).

For each convex feature that has one dangling endpoint we try to extend it into a loop or, if not possible, we trace a loop through that endpoint in the direction orthogonal to the feature curve, adding it to the queue. This step helps eliminate vertices of valence 1 or 2 on the subsequent metamesh which would induce nonhexahedral blocks.
Figure 4 shows some examples of concave line features and their extension to form loops. Note that loops should not pass through singularities (see Appendix A.1), while the endpoints of an open feature are likely to lie at singularities. We get around this issue by tracing two loops and that run parallel to the feature along its two sides, and infinitesimally close to it. The portions of the loops and corresponding to the feature are constrained to follow it. After completion these portions are snapped to , while the remaining portions are left free to follow the field elsewhere. Since the endpoints of might be singularities, the loops and may take different routes (see Figure 4.cd). We favor the formation of loops that encompass more than one line feature (see Figure 4.e), by introducing a bias during loop tracing, which reduces the length of curves running along curve features (see Section 4.3 and Appendix A.2 for details).
We discard self intersecting loops. When tracing loops we prevent any new loop from tangentially intersecting, or partially overlapping with, another loop or feature curve (other than the one it lies on). See Appendix A.1 for a formal definition of tangential and orthogonal intersections. To this end we form a set of constrained edges initialized with all convex feature curves. For any newly formed loop we prevent it from tangentially intersecting curves in or .
4.2. Sampling Loops
We augment the queue with additional loops in order to obtain a fairly regular distributions of loops over and to allow for formation of blocks that satisfy the topological constraints in Section 3.1. We define an arc as a portion of a loop that lies between two consecutive intersections of with other loops of . To satisfy criterion t2 we require each loop to consist of at least three arcs.
4.2.1. Loop Space Distance
In order to evenly sample loops from the space , we introduce a notion of (nonsymmetric) distance between loops, as the average over one loop of the shortest distance from each of its points to the other loop. Let and be two loops, then we define:
(1) 
where is the length of the shortest fieldcoherent geodesic path joining a point of to . In order to get an intuition for this distance, consider that two roughly parallel loops on are close to one another, while loops that are either intersecting orthogonally, or wind about different handles of an object (of nonnull genus) are usually far apart.
Now considering a set of loops , and a loop not belonging to , we can generalize
(2) 
hence the notion of farthest loop from is welldefined as
4.2.2. Loop Insertion
We extend the queue incrementally through a process of farthest point sampling. We form a discretization of the infinite space of loops , which we refer to as the pool, which is built dynamically and used to sample loops to build as follows. We define the set as the set that contains all loops in together with the set of all convex curve features. We form the initial pool as follows.

We sample all curves of at fixed intervals and, for each sample, we trace orthogonal loops that we add to . Traced loops are constrained to avoid tangential intersections with the elements of ;

We perform a Poisson point sampling on [Corsini2012] to obtain a set of seeds and, for each seed , we trace the two orthogonal loops, each constrained by the elements of as above, and we add such loops to .
Each time we add a new loop to , we replenish the pool with new loops obtained by sampling as described in item (1). When adding loops to , we prevent them from tangentially intersecting any loop in .
We add loops from the pool to the queue using the following priority rule. Let be the subset of made of loops that are formed by less than 3 arcs. We give higher priority to those loops in that intersect at least one loop in ; among them, we select the loop that maximizes its distance from all loops in and we add to . Next, we remove from all loops that intersect tangentially, and retrace them from their sources in , constrained to the updated .
This step can be repeated at will and fosters the formation of a queue that contains loops with at least three arcs and are uniformly distributed over . Figure 5 shows a few steps of loop generation (see Section 4.3 for details). Note in Figure 5 how both the collection of loops in and the pool become progressively more dense as the process goes on. The generation of loops is stopped as criteria t2 and g3 are both satisfied (where geometric fidelity g3 is tested on the patches intercepted by the network of loops on mesh ); in order to speedup the process, criterion g3 is not tested until t2 is satisfied at all loops.
4.3. Discrete Loop Tracing
To efficiently and robustly implement the tracing procedure and furthest sampling process described in the two previous sections, we trace fieldcoherent geodesic loops with the discrete graphbased approach of [Pietroni:2016], which brings to the stratified structure (Appendix A.1) the technique of [Lanthier:2001] to evaluate geodesic paths and distances. In short, in [Lanthier:2001] shortest paths and distances are found by a Dijkstra search on an extended graph , which is built over edges and vertices, plus Steiner points sampled on edges and arcs connecting vertices of and Steiner points across each triangle. In our case, four pointarrows are generated per vertex and per Steiner point, which are properly arranged on , and just fieldcoherent arcs connecting them are considered. The graph is built once, and used in all subsequent processing. One important advantage of this method is that crossings and overlaps of paths can be handled in a robust, combinatorial way that does not involve numerics: it is possible to precisely identify whether two paths intersect orthogonally or tangentially by simply comparing arcs that belong to the same triangle of and checking their underlying direction fields. Thus during tracing we prevent any new loops from tangentially intersecting loops in by blocking Dijkstra propagation along arcs in that are orthogonal to those arcs belonging to paths already in .
Loop Space Distance
Given a queue of loops , we set a source for each node of graph traversed by each loop and we run a Dijkstra propagation. The distance from any other loop is computed easily by collecting distances at all nodes traversed by after propagation. Note that a single Dijkstra propagation is sufficient to set distances at all nodes of , thus it is sufficient to evaluate distances from all loops.
Extending curve features to loops
At the beginning of the process, all open line features must be extended to loops, as described in Section 4.1. Since most such features join singularities of , we cannot directly extend them from their endpoints, which are not represented in graph . Therefore, we proceed by tracing loops that run parallel to each feature, and we snap such loops to the features later on. This procedure is described in more details in Appendix A.2.
(a)  (b) 
4.4. Convex Features and MetaMesh Initialization
We initialize the metamesh as being formed of a unique cell, with surface , embedding on it all convex features as edges, and their corners and endpoints as vertices; faces of are the components of that are disconnected by this network of convex features (possibly, there can be a single face with cuts). For each convex feature that has one dangling endpoint we extend the feature along its direction to form a loop. We add the edges along this loop to the metamesh, updating the face structure accordingly.
5. Model Cutting
Our cutting method receives as input the current top loop in the loop priority queue, and uses it to form a cutting surface that partitions the current model into two disjoint parts, refining the intersected blocks and their corresponding metamesh cells. Each cutting surface is bounded by a set of loops that includes . The cutting process involves three major steps: loop grouping, cutting surface computation, and block set and metamesh refinement.
5.1. Loop grouping
While a cut bounded by a single loop is enough to halve a genus zero object, additional loops may be necessary to cut in half higher genus shapes. Moreover, using a single loop may produce highly nonplanar cuts in the presence of deep cavities. Using surfaces with multiple boundary loops in such cases allows for cuts with less geometric distortion (Figure 7). We address both scenarios by forming cutting surfaces with multiple boundary loops, producing better shaped blocks. We group loops together using a combination of topological and geometric criteria.
We observe that, from a cutting perspective, we need to distinguish between three different types of loops, which can be roughly classified according to the interaction between the outward loop curve normals and their coincident surface normals (Figure 8):

Type I loops conceptually lie outside the object, and have outward normals well aligned with the surface normal. These loops may bound a well shaped, or near planar, cutting surface on their own, and form the majority of our cutting surface boundaries.

Type II loops have normals that point in the opposite direction to the surface normals, and are typically located inside tunnels or cavities. To obtain a well shaped cutting surface, these loops always require a cutting mate of type I, which defines the outer boundary of the same cutting surface, and may occasionally be grouped with other loops of same type, generating a surface with multiple holes.

Type III loops have normals that lie close to the surface tangent plane. To provide a nondegenerate partition of the object, they need to be paired with a cutting mate of the same type, to jointly form an approximately cylindrical cutting surface.
Note that this characterization is not formal, but rather a tool to predict the shape of the surfaces that would result from using each loop or a set of loops as a boundary of a cutting surface. In fact, a cut that interpolates two loops of type III and a cut that interpolates one loop of type I and one loop of type II are topologically equivalent (they are both generalized cylinders with two disjoint boundaries); but while in the former case the cylinder has a well defined inner axis, in the latter case it may be completely foreshortened with two coplanar boundaries.
Loop assessment
To assign a unique type to each loop we study the relation between its normals and the surface normals (Figure 8). We compute loop normals using the curve’s Frenet frame, oriented according to the winding order that provides normals that largely point away from the loop’s center of mass. We integrate the dot product between loop and surface normals, and perform the classification by applying a symmetric threshold centered at 0 to discriminate between the three types. Values higher than 0.3 denote type I loops; values lower than 0.3 denote type II loops; and values in between denote type III loops.
Note that if a loop spans a sharp geometric feature there will be a mismatch between the surface normals on one side of the loop and the normals on the opposite side. As explained in Appendix A.2, our loop tracing method always produces two loops along each concave feature on its opposite side. Each loop is infinitesimally close to the feature, and is parallel to it without ever crossing it. Consequently each such loop has a different surface normal and induces different cuts. At the same time, in geometric space the portions of the loops along the feature completely overlap. Thus along closed features such as the two concave rings of the Bearing Plate, the loops generate two cuts orthogonal to one another emanating from the same geometric loop (Figure 9).
Grouping
Starting from a loop , we wish to find the smallest set of complementary cutting loops that jointly with this loop cut the surface in half and produce a cutting surface whose shape is most reflective of the loop’s type. If the loop is of type I or II, we seek to obtain a cutting surface that is as planar as possible, whereas if the loop is of type III we wish to obtain a cylinderlike cut. We reflect this preference using a loop pairing metric defined as follows:
where the terms and measure the similarity between the loops spanning planes or cylinders, respectively. Specifically, is if and span the same approximating plane, and grows proportionally to the angular and euclidean distance between the planes they span. is if and span the same approximating cylinder, and grows proportionally to the angular distance between their axes and the ratio between their cylinder radii. Both metrics are defined in the range and implemented using Gaussian functions. The details on the computation of the approximating planes and cylinders and the metric used to assess their similarity are provided in Appendix B.1.
Given an initial loop , we determine the optimal set of cutting loops by computing a binary labeling of the surface that minimizes
and is constrained to have different labels along the opposite sides of to guarantee that the boundaries of the bipartition include it. We compute the solution using a mincut formulation applied to the dual of mesh . Cuts along mesh edges (arcs in the dual graph) that do not belong to any loop receive costs, whereas a cut along a mesh edge that belongs to some loop contributes to the energy proportionally to the ratio between its length and the total length of the loop it belongs to
The minimizer of automatically provides the smallest set of additional loops that globally halve and are geometrically best aligned with . We compute the solution using the mincut implementation provided by [boykov2001fast; kolmogorov2004energy; boykov2004experimental]. Note that minimizing at local (per edge) level may produce inconsistent global results, where only a subset of the edges participating in some loop is selected for cutting (Figure 10). Empirical observations indicate that this happens only if there are no good geometric matches for in the loop queue . We recover from these pathological situations by using a backup strategy, which consists in minimizing the total length of the cuts, as detailed in Appendix B.2. The backup strategy is purely topological, and essentially cuts along all topological handles that are necessary to globally halve the shape with as short as possible cuts. Since loops are geodesic, cutting along a whole loop in general produces a shorter boundary than partially cutting along multiple loops, and therefore a valid cut set is usually found. If no valid cutting loop set can be found also with the backup strategy, the initial loop is discarded.
Handling Cavities
The loop grouping as described above introduces complementary loops in the cut set only if does not bipartition the surface of alone. If is a separating loop no complementary loops will be found, as they would increase the energy. As a result, loops of type II located inside cavities will not receive a matching loop of type I, producing highly non planar cutting surfaces that induce a poor block decomposition (Figure 7, left). Thus we first locate all type II loops in our queue and apply the grouping process from them to other loops. If one of those loops is subsequently pulled from the queue we use its previously computed group as the cutting surface boundary. To provide a quality decomposition of complex objects containing a variety of tunnels and cavities, several inner loops of Type II may need to match to the same outer loop of Type I at once, producing a single surface cut with many holes inside (Figure 7, right). We obtain the desired effect by imposing specific boundary conditions on the mincut algorithm. For each loop of Type II, we first compute its best geometric match of Type I using . We then cluster together all loops of Type II that match the same Type I outer loop, and each time any of these loops is used for cutting, we force the mincut step to cut through all of them. Similarly, for each loop of Type III that does not surround a handle (i.e. it is contractible), we find its best geometric matching according to , and force mincut to include both such loops in the bipartition.
5.2. Solid cut
Given the set of cutting loops that bipartitions the surface of the object, we extend the cut to the interior by defining it as the level set of a volumetric harmonic function that aligns with . Harmonic functions satisfy the maximum principle, which ensures that if the function is non flat, then maxima and minima arise only at the boundaries of the domain and never in the interior. By prescribing Dirichlet boundary conditions on the surface , we can therefore control the topology of the cutting surface, and guarantee that it interpolates all the cutting loops in , is manifold, and does not interpolate any other points on the model’s surface.
Considering the bipartition of the surface induced by the loop set , we denote as the vertices of these loops, as the vertices to the left of , and as the vertices to the right of . We then solve for the volumetric harmonic function , subject to the following Dirichlet boundary conditions
where the 0.5 level set of defines the cutting surface bounded by our set of loops (Figure 11); and denotes the distance between each surface vertex and its closest loop in the cut set , measured using a Dijkstra search on the graph of edges of the volumetric mesh . The normalization factor is used to bound the function to the range , and is the maximum distance from a surface vertex to the set of loops . We implement the operator as a simple combinatorial Laplacian, using uniform weights for all vertices in the one ring. The choice of combinatorial Laplacian was dictated by the need for solution robustness. Our cutting surfaces are embedded into the connectivity of the tetrahedral mesh and cut across mesh edges. This embedding process progressively worsens the shape of the tetrahedral mesh elements as more cuts are introduced. While the popular cotangent Laplacian is geometryaware and produces smoother level sets on well shaped meshes, our experiment show that using it on poor quality meshes results in illconditioned and numerically unstable linear systems. In contrast, the combinatorial Laplacian is robust to meshing defects, and always yields a wellconditioned positive semidefinite matrix. While in theory it is not geometryaware, in practice the cutting surfaces we obtained were always good enough for our purposes, even when the tetrahedral elements were poorly shaped.
5.3. Metamesh update
The metamesh is initialized as a single volumetric cell, which incorporates on its boundary all the input feature lines as edges, their endpoints as vertices, and the surface patches of that are disconnected by them as faces. Solid cuts are progressively embedded into the connectivity of a tetmesh by splitting all the edges traversed by the 0.5 level set (Figure 12). As a result, each cut is a collection of triangular facets of that separates pairs of tetrahedra incident at them at both sides. Using these facets as boundaries and performing exhaustive region growing on the tetmesh we produce and maintain a labeling of tetrahedra defining the volumetric cells of our decomposition. We extract the faces, edges and vertices of the associated metamesh from this labeling. Each face of the tetmesh at the interface between two different labels (including the null label conventionally assigned to outer space) is labeled with such pair of labels, and all tetmesh faces with the same pair of labels are gathered to form a metamesh face; then edges of the metamesh are obtained by gathering chains of edges of the tetmesh that have the same set of two or more incident faces of the metamesh. Finally, tetmesh vertices incident on three or more edges of the metamesh are classified as metamesh vertices. In case topological inconsistencies arise while intersecting a cutting surface with the current metamesh, we fix them as explained in Appendix B.3.
6. Results
We validate our method on a range of models of different complexity, both mechanical and organic, and showcase 31 fully hexahedral meshes and just one hexdominant mesh throughout this paper. We report numerical statistics in Table 1. Loopy Cuts produces meshes with a complex singularity structure, enabling them to conform to the input features and at the same time provide extremely high perelement quality, often much higher than alternative hexmeshing techniques (e.g. Figure 15). We demonstrate our ability to conform to a variety of features, both geometric (Figure 13) and synthetic (Figure 16). Our pipeline creates pure hexahedral meshes given surface feature networks with all valence three vertices. Our midpoint subdivision introduces non hexahedral elements only on metamesh cells that are incident at vertices with a valence different from three (Figure 17). These configurations seldom occur in practice, and in our experiments we produced only one hexdominant mesh containing four hybrid elements (less than 1% of total elements count). As shown in Figure 14, we generate full hexahedral meshes even when state of the art hexdominant meshing techniques introduce a significant amount of hybrid elements. In general, these techniques consistently produce meshes with a much lower percentage of regular elements (72% to 91% hexahedra for [hybridHexa], and 33% to 95% hexahedra for [Sokolov:2016:HM:2965650.2930662]), and contain complex hybrid elements that may have up to 40 facets (see Table 1 in [hybridHexa]).
Customization
Loopy Cuts is fully automatic, but can optionally be used in interactive mode to allow the user to customize the block decomposition. In Figure 18 we show a result obtained with manual interaction, where we guided the loop grouping and cut sequence to obtain a perfectly symmetric mesh with high anisotropy along one direction. We provide this control by enabling three interactive operations that can be activated via a trivial pointandclick user interface. The user can:

Prescribe a customized cut sequence. Each time a loop is selected for cutting, the system automatically finds complementary loops, and the volumetric cut is performed;

Manually pair loops for cutting. Note that this could also be use to accumulate multiple cuts into one (e.g. the cuts that separate the dents from the pinion wheel in Figure 18 were made all together, while the automatic approach would have processed them in sequence);

Mark a loop as a new convex feature. Such a loop will then be incorporated in the meta mesh as a chain of surface edges that split faces, but the meta mesh will not be cut through it.
An additional indirect form of customization can be obtained by editing the guiding surface field with userinteractive tools such as the one proposed in Instant Meshes [jakob2015instant]. Such tools can be used to prescribe additional soft constraints on the meshing process promoting alignment to secondary features or symmetries (Figure 19).
Implementation details.
We implemented Loopy Cuts as a single threaded C++ application, using Tetgen [si2015tetgen] for tetrahedralization, and Eigen [eigenweb] for numerics. Cross fields aligned to line features and surface curvature were computed using MIQ [bommes2009mixed] and the PolyVector Fields method [ComplexRoots:Diamanti:2014]. Sharp creases were automatically detected by thresholding dihedral angles, whereas other features were manually marked. Note that both field computation and feature detection are external to the algorithm; alternative techniques may also be used, as Loopy Cuts is completely agnostic to how such information are computed. We produce our output meshes by applying one step of midpoint subdivision [li1995hexahedral] to each block of the meta mesh, and optimizing the geometry with edgecone rectification [Livesu:2015]. Our pipeline runs in minutes; the global running time depends directly on the number of cuts, and indirectly on the complexity of the shape. On average, the generation of the cutting loops takes between 10 and 120 seconds on meshes in the range of 4K50K tris (often obtained by featurepreserving remeshing from highresolution meshes); while 3 to 8 seconds per cut were necessary to cut models between 100K and 200K tets.
Model  Cuts  Metamesh  Output mesh  

H  P  O  H  SJ  
3holes  19k  13  64  16  0  1K  .50/.87 
anchor  11k  16  44  28  9  0.7K  .33/.85 
anti backlash nut  13k  34  77  14  16  0.7K  .43/.84 
bamboo pen  36k  20  76  38  0  0.8K  .72/.93 
bearing plate  10k  16  43  12  16  0.7  .84/.98 
beetle  70k  34  154  38  0  1.6K  .50/.95 
blech  60k  42  132  30  0  1.3K  .64/.96 
block  28k  15  33  32  10  0.6K  .41/.82 
cap thing  57k  15  68  18  6  0.8K  .56/.92 
cube minus sphere  4k  3  17  4  1  0.2K  .62/.97 
cube 2 colors  2k  11  56  16  8  0.8K  .50/.96 
spiky cube  6k  22  136  72  56  1.8K  .57/.93 
delta arm base  37k  22  52  38  12  0.9K  .50/.77 
double torus  32k  11  0  88  4  0.6K  .49/.81 
fandisk  7k  17  126  11  1  2K  .48/.85 
hand  20k  29  31  82  37  1K  .33/.83 
hinge  19k  13  56  14  2  0. 6K  .52/.94 
impeller  24k  31  96  96  0  1.7K  .36/.88 
joint  9k  12  56  12  0  0.6K  .77/.97 
knob  14k  10  88  12  0  0.7K  .82/.97 
lever arm  9k  23  36  32  4  1.2K  .54/.92 
lock  42k  22  124  58  18  1.7K  .38/.85 
motor tail  15k  45  177  53  0  2K  .51/.95 
pinion  43k  40  230  30  0  2K  .41/.95 
pinion (manual)  43k  12  40  0  0  0.3K  .66/.95 
rod  9k  15  108  46  4  1.2K  .72/.95 
sculpture  32k  25  80  8  0  0.7K  .77/.97 
shaft  30k  10  20  32  8  0.6K  .44/.92 
torus  4k  12  65  52  0  0.8K  .82/.93 
toothbrush holder  13k  26  43  38  0  0.7K  .62/.93 
trebol  13k  10  0  15  17  0.2K  .47/.77 
wave  3k  11  40  16  16  0.5K  .60/.93 
7. Conclusions
We present a new approach to blockdecomposition for hexahedral meshing. The method combines tracing of strategically placed crossfield coherent cutting surface loops with computation of fair cutting surfaces interpolating one or more such loops. As shown by our results the method outperforms prior work in providing a combination of robustness, feature interpolation, and curvature alignment.
Limitations and Future work
Our framework is not guaranteed to compute a valid decomposition. In particular we envision the following sources of failure. First, our loop formation strategy relies on the underlying crossfield. On surfaces where the crossfield directions change multiple times, the resulting loops may be too complicated or the tracing may not be able to close loops properly, avoiding selfintersections. Second, our cutting surface computation depends on the existence of suitable groups of loops that jointly bound a a desired cutting surface. Sometimes, some loops find no mate to pair with, just because such loop cannot be found on the sole basis of the existing line features and the cross field. Finally, our grouping strategy is heuristic and may fail even when suitable loops exist. We did not exploit user interaction at its full potential, since in our current implementation this is allowed only during the pairing and cutting phase; enabling simple user interaction during loop selection may resolve with a few clicks most of the issues described above, for instance by forcing or deviating some loop, or changing the flow of the crossfield. Addressing all the above aspects is an interesting avenue for future research.
References
Appendix A Computing cutting loops
a.1. Fieldcoherent loops
Following [kalberer2007qsp], the four components of a cross field can be separated on a stratification of manifold into four sheets, defined as follows (inset). For every point of , except the singularities of , consider four copies , , and , each consisting of together with one of the four directions of at , such that and . We call each such copy of a pointarrow meaning that it incorporates both a position on and a direction on its tangent space. Space consists of four sheets, each corresponding to less the singularities of , such that the pointarrows defined above are distributed among the layers to form a smooth direction field (see inset). Generally speaking, if has singularities, the direction field on turns about such singularities, thus sliding between different sheets, and consists of a single connected component. Space is the quotient space of obtained by identifying pairs of pointarrows and , thus consisting of two sheets, each endowing a line field. Manifold can be also seen as a quotient space of , by identifying the four pointarrows at each point .
Following [Pietroni:2016], a smooth (oriented) curve on is said to be a fieldcoherent path if its tangent direction at all points does not differ for more than from the underlying direction field on (pink wedges in the inset). With abuse of notation, we denote by also the corresponding curves on and , regarded as quotient spaces of . Two fieldcoherent paths and are said to intersect tangentially if they intersect in ; while they intersect orthogonally if they intersect in but they do not intersect in .
The drift of a fieldcoherent path w.r.t. comes from the angle between the direction field and the tangent of at each point along it. Following [Pietroni:2016] we adopt an anisotropic metric on that increases the length of a path proportionally to its amount of drift:
(3) 
where
is a tangent vector at
, is its Euclidean norm, is the reference direction on , measures the unsigned angle between a pair of vectors, and is a parameter that tunes the amount of penalty for the drift. A fieldcoherent geodesic path between to pointarrows on is a fieldcoherent path joining them that is shortest according to the above metric. We define a fieldcoherent geodesic loop to be a nonnull fieldcoherent geodesic path that starts and ends at the same point.(a)  (b) 
(c)  (d) 
a.2. Extending features to loops
In order to force loops to run along line features, we modify the graph as follows:

Given a line feature , we create two corridors, each made of a strip of triangles of incident at , one for each side of (see Figure 20.a);

For each face in a corridor, we consider all Steiner nodes of that are coherent with the direction of and that lie on edges crossing the corridor (see Figure 20.b);

We reduce the weight of arcs connecting pairs of such nodes in the corridor (green arcs in Figure 20.c);

We inhibit all arcs that connect such nodes with nodes at the boundary of the corridor (red arcs in Figure 20.d).
For each line feature , we create one seed node per side of and we trace a set of candidate loops from all such nodes. Note that, in the modified graph, each path that enters a corridor is forced to traverse it totally, and paths traversing several corridors (i.e., joining or bridging different line features) are favoured because of their reduced cost (see Figure 6.a). Note that each feature may be traversed by multiple loops in the set of candidates. In the process of generating the loops that extends line features we select loops in a greedy manner, preferring the ones that span the largest length of open features.
Appendix B Model cutting
b.1. Similarity metrics for loop pairing
We formalize here the similarity metrics for planes and cylinders used to perform loop pairing. Both metrics are defined as penalty metrics in the range , meaning that lower values denote higher similarity. Given two loops of type I, II that are centered at and span two planes with normals , respectively, we define their plane similarity as follows
The first term weighs angle similarity between plane normals; the second term weighs the maximum distance between the centroid of and the plane spanned by , and vice versa. Vector is defined as .
Given two loops of type III, that are centered at , have radius and span two oriented lines , respectively, we define their cylinder similarity as follows
The first term measures angle similarity, and promotes pairing between loops laid on oppositely oriented surfaces. The second term matches the radii of the cylinders, and serves to find the best geometric matching in presence of concentric loops.
b.2. Backup pairing strategy
The backup pairing strategy is used when the geometric loop pairing does not provide a globally consistent bipartition of . It computes a bipartition that assigns to each triangle in the mesh either 0 or 1, and minimizes the total length of the edges having opposite labels at their sides
where are the mesh edges, are the labels associated to the triangles incident at it, and the per edge cost energy is defined as
We compute the solution using the mincut implementation provided by [boykov2001fast; kolmogorov2004energy; boykov2004experimental]. To guarantee that the initial cut is used to bipartition, we set boundary conditions that impose that the triangles at its two sides must receive opposite labels.
b.3. Topological cleaning
Since the final metamesh reveals itself cut after cut (Figure 12), extracting it when none or just a few cuts have been performed may result in a number of topological inconsistencies that affect its cells, such as: dangling edges, islands, nonmanifold vertices and faces with less than two vertices (Figure 21). We make the extraction of the metamesh robust against all these defects by using a set of topological operators that collapse all the topologically illegal entities, as depicted in the bottom line of the figure.