Watertight conversion of trimmed CAD surfaces to Clough – Tocher splines ✩

Article history: Received 22 April 2015 Accepted 1 June 2015 Available online 10 June 2015


Introduction
The common representation of shape available to all Computer-Aided Design systems, and therefore the standard form for exchanging CAD models between systems, is a B-rep (boundary representation); (Stroud, 2006). In this paradigm, a connected wireframe of edges and vertices describes the global topology and sharp features of a model. Faces are bounded by a collection of edges and have their internal geometry provided by a smooth 2-manifold surface: usually one of the natural quadrics or a non-uniform rational B-spline (NURBS) surface (Piegl and Tiller, 1997) for more general freeform shapes. The embedding surface typically extends beyond the bounding edges and the portion that is meaningful is therefore a "trimmed" region of the whole. This representation of shape has been successful in realising the broad range of applications that CAD systems support today, for example in visualisation and manufacturing.
However, when we are interested in the region of space bounded (or excluded) by a shape, for example when performing an engineering simulation (Cottrell et al., 2009;Jaxon and Qian, 2014) or computing a medial axis transform (Bucklow, 2014), our task is often made difficult by the deficiencies of such a B-rep. For example, the intersection curve of two ✩ This paper has been recommended for acceptance by Gerald Farin. NURBS surfaces is not, in general, a NURBS curve, leading to unavoidable gaps in trimmed NURBS models (Skytt and Vuong, 2013). Although the gaps can be made arbitrarily small, the result is still a discontinuous representation of shape. This is a significant barrier for the direct use of CAD models in most types of engineering simulation, which require at least continuity of position. Additionally, it is often desirable in downstream applications (such as slicing for 3D printing) that all surfaces use the same representation, continuity, and degree. This simplifies interrogation considerably, and allows for algorithms that use data-level parallelism (e.g. on a GPU) to improve performance.
There is therefore a need to convert B-reps, and particularly the trimmed surfaces that they comprise, into a form that removes the gaps between adjacent faces and allows for a homogeneous representation of shape. This paper investigates the performance of cubic C 1 Clough-Tocher splines for this purpose. These splines give a way of smoothly interpolating surface samples with gradients and can therefore • remove gaps between B-rep faces (giving surfaces that are C 1 within a B-rep face or C 0 when stitched together); • provide a simple, homogeneous representation of freeform shapes; • allow for local adaptivity.
Cubic C 1 Clough-Tocher splines are described in the next section. Various constructions of this type are summarised in Section 3, including previously unexplored variants. Section 4 addresses the problem of generating input data from trimmed NURBS surfaces that are then sampled and interpolated by a Clough-Tocher spline. Our experiments and a comparison of various constructions are given in Section 5. The paper concludes in Section 6 with a short summary and suggestions for future research.

Preliminaries
We start by introducing our notation and some basic concepts including Bézier triangles, the cubic C 1 Clough-Tocher spline space and its reduced variant. More details can be found in the papers by Farin (1986), Lai and Schumaker (2001), Alfeld and Schumaker (2002), and Speleers (2010).
Any polynomial p of total degree at most d defined on T can be expressed in the Berstein-Bézier form with i = (i, j, k), |i| = i + j + k, and i, j, k ≥ 0. These polynomials span the linear space denoted by d . In the cubic case, the Bézier ordinates P i , associated with the barycentric coordinates i/3, are shown in Fig. 1, left. Moving from the above functional setting to the parametric one, a Bézier patch is given by with control points P i ∈ R 3 forming a triangular control net.

Clough-Tocher spline space
Let be a conforming triangulation of a domain ⊂ R 2 with one or more polygonal boundary loops and let n v , n e , and n t be the number of vertices, edges, and triangles of , respectively. Additionally, let E denote the set of edges in and let the vertices U i of have Cartesian coordinates (also called parameter values) given by (u i , v i ).
In the Clough-Tocher construction, each of the triangles in , which we call macro-triangles, is partitioned into three micro-triangles (Clough and Tocher, 1965); see Fig. 2. For each triangle T ∈ , a split point Z is chosen and connected to the three vertices of T by new edges, giving rise to the Clough-Tocher refinement of , denoted by .
The cubic C 1 spline space defined on is called the Clough-Tocher spline space Its dimension is equal to 3n v +n e . Note that the dimension of the space S 1 3 ( ) over a general triangulation is still an open question, so the extra freedom granted by a Clough-Tocher refinement is useful to construct C 1 cubic splines in practice.
We will call the restriction s| T the macro-patch over T .
The three degrees of freedom corresponding to each vertex in are typically fixed by assigning the position and gradient at each vertex U i , i = 1, . . . , n v . That is, in the parametric setting, sampled at the vertices of from a parametric surface f defined over . This leaves a degree of freedom left for each edge in . This degree of freedom can be used to force the directional derivative of the spline s along an edge ε to be linear, instead of the general quadratic. More precisely, the reduced Clough-Tocher spline space is defined as (Speleers, 2010) where v ε is a (unit) vector not parallel to ε. This reduced space has dimension equal to 3n v and provides a unique solution to the interpolation problem of (9). A discussion on choosing the split point Z can be found in Schumaker and Speleers (2010). We consider three options for the split point Z : the usual barycentre of T , the incentre of T , and a previously unexplored option, the point in T that is given by the barycentric coordinates of the incentre of the corresponding triangle in 3D given by f i . These three options will be referred to by Bary, Inc2, and Inc3, respectively. We discuss the consequences of these options in Section 5.

Cubic C 1 Clough-Tocher spline constructions
An excellent survey of various C 1 cubic Clough-Tocher spline constructions was conducted by Kashyap (1996). That survey also included a new construction and a discussion of iterative methods.
We now briefly recall the most prominent C 1 cubic constructions and their novel variants, and, for convenience, provide formulas for computing all the necessary Bézier ordinates (control points).
From the discussion of C 1 continuity conditions (see Section 2.2 and Farin, 1985Farin, , 1986, one obtains the following three-step procedure for computing the control points of micro-triangles (see Fig. 2, right) from the input data (9), applied to each T ∈ . The triple (τ 0 , τ 1 , τ 2 ) provides the barycentric coordinates of the split point Z in T (U 0 , U 1 , U 2 ).
Step 1. Set and similarly for the remaining control points along the edges of T . Then compute and similarly for I 11 and I 21 .
Step 2. Compute C 0 , C 1 , and C 2 by one of the constructions detailed below.
The first step ensures interpolation of the input data. For given split points, the second step fixes the only degree of freedom in the construction, i.e., it computes the central control point of each micro-triangle. The construction is finalised in the third step, which computes all the remaining control points to assure global C 1 continuity.
We now move on to particular constructions and show how to perform Step 2 for each of them in turn. Although triangles with boundary edges require special treatment in most constructions mentioned below, we address boundaries later in Section 3.8.

CT o : the Clough-Tocher construction
The original construction of Clough and Tocher (1965) fixes the extra degree of freedom per edge in by requiring that the normal derivative at each edge be linear. Thus, this construction yields splines from the reduced space Ŝ 1 3 , see (10), by choosing v ε ⊥ ε. We will refer to this construction as CT o , where the 'o' stands for 'orthogonal'.
A formula for computing C 2 based on a general direction v ε , derived from the linear derivative requirement of Ŝ 1 3 , is given in Section 2.3 of Speleers (2010): with (λ 0 , λ 1 , 0) the barycentric coordinates of the projection of Z onto the edge ε = U 0 U 1 in the direction of v ε .
In the special case of v ε ⊥ ε, λ 0 and λ 1 are computed by orthogonal projection as Clough and Tocher (1965) mention several variants of this construction and note that the original 3-split into microtriangles was suggested by T.K. Hsieh.
The advantage of this original construction is that it is completely local. More precisely, all control points for the three micro-patches per input triangle T are computed using only the input data at the vertices of T . This makes it very simple to implement. On the other hand, the extra degree of freedom per edge is not used in an optimal way. Additionally, using the perpendicular means that the construction is not invariant with respect to affine reparameterisations, not even to non-uniform scaling of the parameter plane. While full invariance in the above sense may seem unnecessary, non-uniform rescaling of parameter space is a common operation in CAD kernels and it is therefore desirable to have this type of invariance. This is elaborated further in Section 5.1.

CT i : an invariant modification of CT o
The idea developed here was mentioned by Farin (1985) and later also by Schumaker and Speleers (2010) and Speleers (2010). Instead of using the perpendicular direction, i.e., v ε ⊥ ε, we set v ε to be parallel to the line connecting Z and Z in the neighbouring triangle T .
In this case, λ 0 and λ 1 in (15) need to be set to where × denotes the planar cross-product of two vectors. This formula is easily derived by intersecting the line ZZ with the edge U 0 U 1 . This construction is invariant to affine reparameterisations, hence we denote it by CT i , where the 'i' stands for 'invariant'. On the other hand, it requires the positions of split points from adjacent macro-triangles and thus it is a less local construction than CT o . Farin (1985) describes a construction that has quadratic precision and minimises the C 2 discontinuity between microtriangles along edges of macro-triangles. We abbreviate this construction by Fa.

Fa: Farin's construction
First, set Q of the macro-patch corresponding to T to guarantee quadratic precision (see Fig. 1, right): All the other control points of the macro-patch are computed from the input data at vertices of T via (11). The macro-patch is then subdivided into three micro-patches according to the split point Z , e.g. using the de-Casteljau algorithm (Böhm et al., 1984). This gives initial positions of all micro-triangle control points and similarly for the remaining points. Note that S is actually not needed in the construction because it is later replaced in Step 3. By formulating a constrained minimisation problem using (7), one obtains, e.g. by the standard Lagrange multipliers method, where The barycentric coordinates involved in (20)-(21) are defined by Z = τ 0 Z + τ 1 U 0 + τ 2 U 1 and Z =τ 0Z +τ 1 U 1 +τ 2 U 0 . Note that this procedure is done separately for each non-boundary edge using information from edge-adjacent triangles. The other control points, C 0 and C 1 , are computed similarly.

FO: Foley-Opitz construction
The method of Foley and Opitz (1992) constructs a hybrid cubic patch with additional control points that allow for crossboundary derivative control. Mann (1998) described how that construction can be modified to fit into the Clough-Tocher framework. We refer to this construction by FO.
First, repeat the following procedure three times, once for each edge of T : set Q in the macro-patch corresponding to T to guarantee cubic precision across the edge U 0 U 1 : Similarly to Farin's method (Section 3.3), the macro-patch is then subdivided into micro-patches via (19), but only C 2 is actually required. The other two applications result in C 0 and C 1 . As shown in Mann (1998), this construction has cubic precision in the sense that if the positions and gradients at the 6 vertices shown in Fig. 3, left, are sampled from a common cubic, the three micro-triangles corresponding to T describe patches of the same cubic. Additionally, it is, by construction, invariant to affine reparameterisations.

Ka: Kashyap's construction
Farin's method described in Section 3.3 minimises the C 2 discontinuity between micro-triangles along edges of macrotriangles. To achieve cubic precision, Farin's minimisation can be used to minimise the C 2 discontinuity of macro-patches, which are in turn used to obtain C 0 , C 1 , and C 2 . As in the case of FO, this procedure is applied three times.
First, fit a macro-patch over T that minimises the C 2 discontinuity across the edge given by U 0 U 1 . This is done using the same equations (18)-(21) as in Farin's method, with the following modifications: and the barycentric coordinates involved are defined by U 3 = τ 0 U 2 + τ 1 U 0 + τ 2 U 1 and U 2 =τ 0 U 3 +τ 1 U 1 +τ 2 U 0 . Note that they are not, in general, the same as the coordinates used in Farin's method in Section 3.3. Repeating this procedure for the other two edges of T gives the remaining control points.
This construction, which we denote by Ka, first appeared in Kashyap (1996) and is based on the 18-point interpolant introduced there. It was later rediscovered in the presented form in Mann (1998Mann ( , 1999. It has cubic precision in the same sense as FO of Section 3.4; see Fig. 3, left.

MG: mid-edge gradients
All the methods discussed so far make use of the data (9) defined on vertices of only. However, if we want to convert a given surface into a Clough-Tocher spline, we can fix the extra degree of freedom per edge in by reading extra data off the given surface. Given the role of the control points computed in Step 2, it is natural to try to use them to fit the gradients at mid-points of edges in . However, as only one degree of freedom is available per edge, the mid-edge gradient needs to be projected in a certain direction v ε not parallel to ε.
There are two obvious options available. One is to set v ε ⊥ ε, leading to a construction denoted by MG o . This was investigated e.g. in Schmidt et al. (2001), Bastian and Schmidt (2001). Similarly to the original Clough-Tocher construction, this ensures complete locality, but does not provide invariance to affine reparameterisations. The other option is to set v ε = ZZ , cf. Section 3.2, leading to a construction MG i that is invariant to affine reparameterisations. Turning back to (5) with d = 3, the gradient of the patch at (U 0 + U 1 )/2 is given by with A(U 0 , U 1 , U 2 ) again equal to the area of T . The central control point P 1,1,1 is then easily computed from the linear where g ε is the mid-edge gradient sampled from the input surface. This covers both options mentioned above. P 1,1,1 , computed three times, once for each edge of T , yields the control points C 0 , C 1 , and C 2 in Step 2. Input data for cubic reproduction over a macro-triangle are shown in Fig. 3, right. A more general treatment of cross-boundary derivatives for triangular patches can be found in Foley et al. (1993).

Iterative methods
We now turn our attention to the iterative method described in Farin and Kashyap (1992) and Kashyap (1996), which is based on Farin's construction in Section 3.3. Although using Farin's method in Step 2 minimises the C 2 discontinuity between adjacent patches, Step 3 subsequently changes the positions of some of the control points involved in the C 2 conditions (7) to achieve global C 1 continuity. Thus, Step 2 can be performed again, then Step 3, and so on.
This gives an iterative process. It weakens the locality of the construction, but typically produces results of better visual quality; see Section 5. We note that this iterative procedure, one step of which is denoted by IM, can be applied to any of the Clough-Tocher variants described above.
A second procedure, presented in Section 4.2.2 of Kashyap (1996), minimises the C 2 discontinuities within macropatches. The result is a single cubic patch over each macro-triangle, but with only C 0 continuity between adjacent macro-patches. While it was suggested in Kashyap (1996) that this procedure can be used in combination with the former iterative method, we did not observe any significant improvement in any of our experiments based on this combination. The adjustment of control points when followed by a step of IM is so small that one cannot observe any visual difference (not even using reflection lines or curvature plots). We thus do not include this second iterative procedure in our study in Section 5.

Boundary triangles
So far, we have not considered triangles whose edges are boundary edges of . Note that only triangles with at least one edge on the boundary of are considered as boundary triangles. From the set of constructions detailed above, only CT o and MG o , due to their complete locality per triangle, need no modification near boundaries. All the other constructions cannot be used over boundary triangles in the form they are defined.
If T is a boundary triangle then its neighbour across its boundary edge ε is not available. In that case, we first consider two options based on (15). One is to use (16), i.e., the perpendicular direction, to compute C 2 . The other is to use the edge midpoint (λ 0 = λ 1 = 1/2), which does not preclude invariance to affine reparameterisations.
If mid-edge gradients are known on boundary edges (i.e., not necessarily on all edges), one can use (25) with either v ε ⊥ ε or v ε in the direction given by Z , the split point of T , and the boundary edge midpoint. Again, only the second choice may lead to invariance to affine reparameterisations. While it is likely that the degree of freedom per boundary edge can be put to a better use in some of the constructions, we do not investigate this any further.

Delaunay triangulation
Before we can compare the above variants of the original Clough-Tocher construction, we need to generate suitable input data. We therefore turn to the problem of generating triangulations with the necessary gradient information from trimmed surfaces appearing in CAD models. While there are many suitable triangulation algorithms (Piegl and Richard, 1995;Piegl and Tiller, 1998;Kumar et al., 2001), the examples in this paper used the CADfix product (TranscenData, 2015) to import CAD models and generate the required data.
CADfix performs the triangulation in the parameter space of each surface, using an incremental constrained Delaunay triangulation (CDT) algorithm similar to the approach described by Kallmann et al. (2004). For NURBS surfaces, the parameter space is naturally given by the rectangle in R 2 that is the tensor product of the spline domains in each parametric direction. For other surface types, a parameterisation is constructed that ensures that the edges bounding each face project to closed loops in parameter space. The input to the triangulation algorithm also includes a maximum permissible error between the surface and its meshed representation, and a flag that indicates whether there are constraints on the aspect ratio of the constructed triangles; see Fig. 4.
Given this input, the triangulation algorithm starts by polygonising each boundary edge, beginning with a small number of samples and adaptively adding samples in any spans where the current representation exceeds the permissible error. The representation of each edge constructed in this stage is a C 1 cubic spline with derivatives that are also sampled from the edge. Each sample of a boundary loop corresponds to a parameter-space node in R 2 , and every pair of adjacent samples has a constrained segment added in the CDT algorithm. The output of this stage is therefore a triangulation of the face which satisfies the error bound on the edges but may be a very poor representation of the internal surface geometry. Note that this edge-first triangulation strategy guarantees that adjacent surfaces share the same boundary vertices and thus makes it possible to produce watertight models.
The second stage of the algorithm iteratively refines the triangulation by examining the error at triangle centroids and triangle edge midpoints, until all triangles have an error that lies below the given bound. For all the examples in this paper, Table 1 Property comparison. See Section 5.1 for a detailed explanation of the meanings of the rows. In each row, best ranked constructions are highlighted in grey. All criteria apply to all three split point options, with the exception of invariance to affine reparameterisations in the second row. In the third row, data required from edge-adjacent triangles are listed; the index ranges in {4, 5, 6}; see Fig. 3 11/0 11/0 11/0 11/0 11/0 1 1 /3 1 1 /3 error was estimated with respect to the CT o construction (with Bary), and the same triangulation and input data are fixed for any comparisons between constructions. Note that this error estimate, while being relatively cheap to compute, does not provide an upper bound. An upper bound on the error could be obtained, if needed, by using local convex hulls provided by the Bézier control nets of Clough-Tocher micro-patches; cf. (5).
If required, the constraints on triangle aspect ratio are also satisfied in this second stage by inserting additional samples that are purely to split triangles with a high aspect ratio; see Fig. 4, right. Unfortunately, at the end of this process, there can still be a mismatch between the C 1 cubic splines constructed for each edge and the boundaries of the Clough-Tocher spline surface constructed for each face. An example is shown in Fig. 18. This is a necessary consequence of the geometry of the trimming curves, the image of which is generally not a cubic spline in the spline surface; generically, only straight edges in parameter space map to cubic curves on Clough-Tocher splines.
We can repair this mismatch rather inelegantly, by simply moving the boundary control points constructed for the Clough-Tocher spline to the positions given by the edge cubic spline control points. Following the description given in Section 3, this involves a different value for T 01 and T 10 (for example) in calculating Step 1. This modification restores C 0 continuity between adjacent Clough-Tocher splines, at the cost of introducing a jump in derivative at the edges of any triangle touching the boundary. In practice, we find that the modifications required by this correction are typically in-plane with the surface, so the spline is usually close to satisfying G 1 continuity conditions even though it is no longer C 1 . This is shown in Fig. 18.

Comparison and results
We now compare and contrast the constructions introduced in Section 3. Our assessment criteria include invariance to affine reparameterisations, polynomial reproduction, visual quality (reflection lines, Gaussian and mean curvature), and approximation error. Table 1 presents a quantitative comparison of the seven Clough-Tocher variants. The best ranked entries for each considered criterion (row) are shaded. None of the constructions performs best across all criteria and only FO and Ka perform best in three out of four of them. Note, however, that different applications may put emphasis on different criteria.

Quantitative comparison
We now address each criterion in more detail.
Polynomial reproduction. Reproduction of polynomials of a certain degree is an important property for nearly any spline construction. Higher polynomial reproduction leads to smaller approximation errors and faster convergence rates under refinement (Lai and Schumaker, 1998). Additionally, cubic precision, instead of quadratic, generally leads to fairer surfaces as we show in the examples below.
Invariance. As already mentioned above, shape invariance with respect to affine reparameterisations guarantees that a rescaling (or any other affine transform) of parameter space does not change the surface. More precisely, let A be a regular 2 × 2 matrix and T a translation vector. Then all vertices U i of are mapped to AU i + T and all gradients ∇f i are transformed to Ā ∇f i , where Ā is the inverse transpose of A (this follows from the chain rule). All positions f i in 3D remain unchanged. Then, a spline s(u) is invariant to affine reparameterisations if s(U ) =s(AU + T ) for all values of U ∈ , where s uses the transformed data. We investigate this property further in the example shown in Fig. 9. Clearly, using Inc2 as split points precludes this invariance. We thus include only Bary and Inc3 in Table 1.
A related concept is that of affine covariance (often also called invariance) in 3D. When an affine transform is applied to the input data in 3D, the resulting spline is the same affine image of the one produced by the original input data. It follows from the definitions of the Clough-Tocher constructions that unless incentres from 3D are used, all variants are affinely covariant in the above sense. This property is not as important as the former invariance in our context; our focus is on converting existing models, not their affine images. And although the input NURBS models are affinely covariant in 3D, the triangulation algorithm used to generate our input data is not.  Farin's cubic given in (26). From left to right: Function values and gradients were sampled at a regular grid of 7 × 7 vertices over the unit square in parameter space; only the highlighted region is reconstructed to avoid boundary influence. A 3D view of the cubic with the triangular grid superimposed. Edges of individual micro-triangles (Ka with Bary). The cubic without the outer-most layer. Locality. Complete locality to one triangle simplifies implementation and boundary triangles do not require special treatment. On the other hand, constructions which make use of data from adjacent triangles lead to fairer surfaces, as we demonstrate in our examples below.
Only CT o and MG o possess complete locality. More precisely, the construction of the macro-patch over T is based on data attached to T only. CT i and MG i have the same locality with respect to position and gradients in 3D, but require the split points of edge-adjacent triangles. Finally, all the other constructions need 3D data from adjacent triangles. The data required from edge-adjacent triangles are listed in Table 1.
Storage. Note that the 11 floats per vertex in Table 1 correspond to 3 floats for 3D position plus 6 for gradients plus 2 for parameter values.
Storing 3 extra floats per edge can hardly be considered a serious disadvantage. However, if full gradients, not just their projections in a particular direction, are stored, 11/3 becomes 11/6. An interesting compromise, as already mentioned in Section 3.8, is to store only boundary mid-edge gradients (as opposed to all mid-edge gradients).
Having explored theoretical aspects of Clough-Tocher spline variants, we now proceed to their visual assessment and error comparison on several examples, ranging from simple test cases to full trimmed NURBS models.
The errors presented throughout the paper were, unless specified otherwise, computed as point to mesh distances over dense samples of both the input function (or model) and the Clough-Tocher approximation using the Computational Geometry Algorithms Library (CGAL, 2015). The errors are reported relative to the models' bounding box diagonals.

Academic examples
Our first example is Farin's cubic function (Farin, 1985, Section 6) c which we used to test polynomial precision of Clough-Tocher spline variants. The general setup is described in Fig. 5. Reflection lines, and Gaussian and mean curvature plots are shown in Fig. 6. The results are in accord with the expected respective quadratic and cubic precision shown in Table 1. While barycentres were used as split points in Fig. 6, we emphasise that the split point does not influence the precision of a particular construction. Additionally, running the iterative algorithm IM on the three constructions with only quadratic precision (first three columns) will force them to converge to the input cubic.
Applying IM to any of the cubic precision constructions will not modify the reproduced cubic (which already minimises C 2 discontinuities).
Our second example is one of the Franke bivariate test functions (Franke, 1979): Maximum errors are summarised in Table 2, computed as function value differences. As all error plots over the whole spline look very similar across all constructions, we show only one of them, for Ka with barycentres, in Fig. 7. Constructions that need special treatment over boundary triangles were modified by using edge perpendiculars as described in Section 3.8. Using edge midpoints gave nearly identical results.
Shaded renderings and mean curvature plots on the whole domain, and reflection lines on a subpatch are shown in Fig. 8; all seven constructions are included. The split point has negligible influence on the resulting shape and thus results using Bary only are shown. The influence of split points is discussed further in Section 5.3. Gaussian curvature plots are not shown as they did not provide sufficient visual differences.
Observe that while both MG o and MG i give the smallest maximum errors (Table 2), they lead to visually poor results (Fig. 8). On the other hand, fairest shapes were achieved by Fa, FO, and Ka. These results can be further improved by applying IM several times.
As our next example, we consider a simple sine wave and data sampled from it using a sparse triangulation. As can be observed in Fig. 9, best results with CT o are achieved when the 2D and 3D triangulation 'agree' as much as possible (i.e., when the parameterisation is close to isometric). The example shows that while CT o can produce worse results than its invariant modification CT i in some extreme situations, proper scaling in parameter space leads to better shapes using the original Clough-Tocher construction CT o . We obtained similar results for MG o when compared with its invariant counterpart MG i .
One approach to quantifying fairness of Clough-Tocher splines is to consider the maximum magnitude of the differences of the right-hand and left-hand sides in (7), i.e., the maximum magnitude of this C 2 discontinuity measure across edges between micro-and macro-patches. This measure is included in Fig. 9. Note how the C 2 discontinuity values closely correspond to the visual comparison using reflection lines and the significant improvement after affine reparameterisation.
Our next test case is the sphere. We use the rational quartic Bézier triangle representation (Farin et al., 1987) of one of its octants and a tessellation of its equilateral parametric triangle into 4 l macro-triangles, where l denotes the tessellation  [−27.48, 39.47]. In this case, grey-scale (black to white) was used because it better reveals curvature discontinuities. Bottom row: Reflection lines on the close-up part shown in Fig. 7. All results are based on Bary. The rightmost plot shows Ka after 10 steps of the iterative procedure IM.

Fig. 9.
A simple wave example. The triangulation in 3D (top right) remains fixed, but we consider two parameterisations, 1 and 2 , which are related by a non-uniform scaling in parameter space. As expected, CT i gives identical results on both parameterisations due to its invariance to affine reparameterisations. In contrast, the shape and fairness of CT o are heavily influenced by the parameterisation; compare CT o ( 1 ) and CT o ( 2 ). Top views with reflection lines are shown. All three examples were generated using Bary. The four numbers under each reflection line plot are C 2 discontinuities: mean and maximum C 2 discontinuity at macro-edges and the same two measures at micro-edges, in that order.
level. Convergence results for Gaussian and mean curvature, and approximation error are reported in Fig. 10 for levels l ∈ [1, . . . , 5], i.e., over regular triangulations with 4 to 1024 triangles. The case with l = 0, i.e., with only one macro-triangle, is not included, as a single triangle gives no information on the fairness between macro-patches.
A visual comparison of Clough-Tocher approximations of an octant of the unit sphere for l = 2, in terms of their fairness and approximation error, is shown in Fig. 11. FO produced indistinguishable results from those generated by Ka. KaG denotes the variant of Ka where mid-edge gradients (Section 3.8) are used at boundary edges (both projection options gave indistinguishable results); the improvement over Ka is significant along boundary edges and it compares favourably even to MG o , which relies on mid-edge gradients at all edges. We observed this behaviour across all tessellation levels; see Fig. 10. Note that the error plots (on logarithmic scale) exhibit the expected convergence rates, i.e., cubic for CT o , FO, and Ka, which possess quadratic precision (in the case of Ka, this is due to boundary effects), and quartic for MG o and KaG, which have cubic precision over the entire triangulation.
Fa, while superior to CT o , produced worse results than FO, Ka, and KaG, both visually and error-wise, and is therefore not included. Only barycentres were considered in all constructions because other split points gave (nearly) the same results due to symmetry. For the same reason, CT i and MG i are not shown as they are indistinguishable from their included counterparts CT o and MG o , respectively.

CAD models
Our first CAD model example is the car front wing model shown in Fig. 12. We used only one of its trimmed NURBS patches in our tests. We consider two triangulations, C and Q , both of which are shown in parameter space in Fig. 4. The corresponding spline approximations with reflection lines are shown in Fig. 13.
While the triangulation C may seem inferior because it contains skinny triangles, it leads to superior results in terms of lateral artifacts when compared to Q . These artifacts are caused by the well-known 'dinosaur back' effect, which arises in nearly all spline approximations when control meshes are not aligned with features on a model (Farin et al., 2002;Augsdörfer et al., 2011). Additionally, once a triangulation has been fixed, the macro-edges of all macro-patches in 3D are uniquely determined and shared by all Clough-Tocher variants. On the other hand, skinny triangles may give rise to foldovers in the resulting spline; see Fig. 14, top row. However, in our tests it was possible to avoid foldovers by using Inc2   . 15. C 2 discontinuity graphs for the car front wing model. Solid lines represent C 2 discontinuities across macro-edges and dashed lines correspond to C 2 discontinuities across micro-edges. The two graphs on the left are for C , the two on the right for Q . The horizontal axes in all four graphs state the number of steps of IM applied. FO with mid-edge gradients at boundary edges produced nearly indistinguishable results from those shown for KaG.
or Inc3. Additionally, Inc2 or Inc3 help to improve the quality of the results over Bary in macro-patches corresponding to skinny triangles (Fig. 14, bottom row).
We gathered C 2 discontinuity data of the constructions and also after several applications of IM. The results are summarised in Fig. 15, including mean and maximal C 2 discontinuities. Observe that while the iterative procedure IM reduces C 2 discontinuities across macro-edges, it is not necessarily the case for micro-edges. Additionally, the maximal C 2 discontinuities across micro-edges are typically greater than those across macro-edges. Errors for all Clough-Tocher spline variants considered in this paper are reported in Table 3 for C and Q .  As we have seen in Section 5.2, the split point position has little influence on the resulting shape when the triangulations are close to regular. On the other hand, as shown in Fig. 14, this influence is much higher for highly irregular meshes containing skinny triangles. We explore this in more detail in Fig. 16. The top row shows triangles in 2D and 3D which are of similar shape. All cases can be expected to lead to good splits of the macro-patch into micro-patches. In contrast, the bottom row shows two triangles that are of very different shapes. While the red split point for Inc2 is sensibly placed in 2D, its blue image in 3D lies too close to one of the macro-edges, giving rise to a skinny micro-patch and potential artifacts. Furthermore, Bary may give rise to poorly-shaped micro-patches caused by uneven distribution of angles between inner micro-patch edges at the image of the split point on the spline; see Fig. 14, top.
In summary, Inc3 improves fairness of Clough-Tocher splines in some extreme cases, but Bary is generally the best (and also simplest) choice. Moreover, recall that Inc2 precludes invariance to affine reparameterisations (Section 5.1). Our second CAD model example is the toy car model shown in Fig. 17. The model is topologically closed but not watertight due to trimming. The triangle layout used for conversion to Clough-Tocher splines is visualised in Fig. 18, left. The middle image shows the result when each input NURBS patch is treated separately. This leads to gaps (as discussed in Section 4). These gaps can be avoided by repositioning some of the boundary vertices in Step 1, giving a watertight model (Fig. 18, right). Observe that the adjustment reduces the continuity of the spline in the immediate neighbourhood of the B-rep edges from C 1 to C 0 , but the magnitude of the discontinuity remains relatively low. While not ideal, it is a considerable improvement over the version with gaps in Fig. 18, middle, and provides an analysis-ready model. The effect of the adjustment leading to watertight models composed of several Clough-Tocher splines is further evaluated in Fig. 19.  Fig. 18. The car model from Fig. 17 converted into a collection of 42 Clough-Tocher splines using Ka with Bary. Left: The layout of macro-patches. Middle: A conversion result without adjusting control points corresponding to input B-rep edges; note the gaps. Right: A watertight conversion result after adjusting control points to match common cubic splines which approximate B-rep edges; see Section 4. Fig. 19. The car model converted to a collection of 42 splines. The adjustment achieves watertightness at the expense of reducing the continuity across some internal edges from C 1 to C 0 ; this corresponds to discontinuities in reflection lines. To maintain invariance to affine reparameterisations for all the methods shown (all based on Inc3), boundary triangles used the midpoint modification described in Section 3.8 where necessary.

Summary
We now present a brief summary of our comparative study. Our main observations are: • CT o , CT i , and Fa produced worst results in the majority of our experiments, mostly due to the fact that they possess only quadratic reproduction.
• With mid-edge gradients available, MG o and MG i resulted in the lowest approximation error, since more information from the converted surface is used compared to constructions that do not rely on these extra gradients. On the other hand, while Ka and FO often give nearly indistinguishable results, they also performed best in most cases in terms of fairness, with Ka slightly better when the adjustment of control points at boundary edges is used.
• The construction KaG which stores mid-edge gradients only at boundary edges is therefore a good compromise between approximation error and fairness.
• Inc3 split points lead to fairer surfaces than Bary in extreme cases such as the skinny triangles shown in Fig. 14. We do not recommend Inc2 split points, as they performed worse than Inc3 in most cases.
• The iterative minimisation of C 2 discontinuity (IM) does improve fairness, but typically only across the macro-edges that appear in the original triangulation, not across micro-edges.
• The parameterisation and triangulation that are used have a significant impact on the resulting shape. While triangles with a low aspect ratio seem desirable for avoiding folded surfaces (Fig. 14), a large number of triangles or poor alignment to model features can lead to much worse surfaces (see Fig. 13). This makes it difficult to design an optimal triangulation strategy.

Conclusion and future work
We have performed a comparative study of seven C 1 cubic Clough-Tocher spline constructions, combined with three split point types. An iterative method aimed at improving C 2 discontinuities, which can be applied to any spline variant, was also included.
We compared the constructions on visual fairness, invariance to affine reparameterisations, polynomial precision and approximation error. The constructions were tested on both synthetic data and on triangulated CAD models. Our results show that no construction is optimal across all assessment criteria and that surface quality depends heavily on the triangulation and parameterisation that are used. Nevertheless, if only one construction was to be chosen, we would recommend KaG, the construction of Kashyap (1996) augmented by boundary mid-edge gradients.
Avenues for future research include improving the stitching process of several Clough-Tocher splines which are separately parameterised. We achieved C 0 continuity by adjusting some of the boundary control points, however, it should be possible to join several splines with at least G 1 continuity across their interfaces without compromising their internal C 1 continuity.
Our tests showed strong boundary effects for the majority of constructions. We believe that boundary influence should be investigated in more detail and freedoms at boundary edges put to better use, ideally in conjunction with the boundary adjustment mentioned above.
The triangulation and parameterisation that are used strongly influence the fairness of the resulting approximations. It would be interesting to investigate this influence in more detail. This is also linked to the use of Clough-Tocher splines in isogeometric analysis (Jaxon and Qian, 2014).