Multi-sided B-spline surfaces over curved, multi-connected domains

We propose a new surface representation, the Generalized B-spline (GBS) patch, that combines ribbon interpolants given in B-spline form. A GBS patch can connect to tensor-product B-spline surfaces with arbitrary G m continuity. It supports ribbons not only along the perimeter loop, but also around holes in the interior of the patches. This is a follow-up paper of a recent publication (Várady et al., 2020) that described multi-sided Bézier surfaces over curved multi-sided domains. While the fundamental concept is retained, several new details have been elaborated. The weighting functions are modiﬁed to be products of B-spline and Bernstein basis functions, multiplied by rational terms. A new local parameterization method is introduced using harmonic functions, that handles periodic hole loops, as well. Interior shape control is adapted to the B-spline representation of the ribbons. Several examples illustrate the capabilities of the proposed scheme.


Introduction
Complex free-form objects in CAGD are predominantly composed of four-sided tensor-product surfaces. However, difficulties may arise, when smoothly connected multi-sided patches need to be created. In particular, problems emerge in design situations where the boundary curves are concave with high curvature variation, and explicit boundary and cross-derivative constraints need to be satisfied for hole loops in the interior of a surface. To overcome these problems we propose a new multi-sided surface representation, the Generalized B-spline (GBS) surface, that combines ribbon interpolants given in B-spline form and connects to tensor-product B-splines with arbitrary G m continuity. This is a follow-up paper of a recent publication (Várady et al., 2020) that described multi-sided Bézier surfaces over curved domains. In that paper, we have argued that the use of curved domains is practically indispensable in complex situations. We combined Bézier ribbons and produced multi-connected C ∞ surfaces, using a parameterization based on barycentric coordinates. We have also pointed out that the Bézier boundaries may need to be split into segments (i) in order to produce tighter control structures, and (ii) to avoid unfavorable parameterizations when a ribbon may excessively penetrate into the patch (see Section 6 in the previous paper).
For the above piecewise representation we applied an overlapping parameterization (as it will be explained later in Section 4), but we also saw the prospective advantages of using B-splines of arbitrary degrees and knot vectors, with a natural local parameterization. This has led to the current scheme, where we combine open (clamped) and closed (periodic) B-spline ribbons. Open ribbons connect two adjacent convex corners; closed ribbons define hole loops in the interior. The new parameterization builds on harmonic functions, and supports periodic B-splines. The patch equation is modified, and new blending functions are used, which are products of B-spline and Bernstein basis functions, multiplied by rational terms.
Curvenet-based design, where complex, multi-sided surfaces are created by interpolation, provides ample evidence for the benefits of this representation. We show three simple examples here. The first model is bounded by a winding, concave B-spline ribbon, and another representing a planar profile on an extruded surface, see Figure 1a. It would be hard to produce such a constrained patch by trimming; in contrast, the GBS patch yields a nice shape. The base surface of the second model ( Figure 1b) has several concave B-spline boundaries. A periodic B-spline ribbon with vertical cross-derivatives is added in the middle, and the base patch is smoothly connected. Traditional CAD would require trimming and lofting operations to reach this configuration, while the GBS patch produces a nice, multi-connected shape without auxiliary steps. The third model is a half bottle ( Figure 1c); in the middle there is a label area that is defined on a developable surface. We prescribe an interior ribbon that matches this surface and produce a single GBS patch by creating smooth blends to the ribbons of the bottle. The paper is structured as follows. In Section 2 we extend our previous review on prior work, then in Section 3 we describe the new surface equation. In Section 4 we briefly review the curved domain creation algorithm, then focus on the new harmonic parameterization. In Section 5 we briefly review methods for interior shape control and in Section 6 we analyze our work by further test examples.

Previous work
There are various techniques to define general topology surface models, including the composition of trimmed quadrilateral surfaces, polyhedral design and curvenet-based design. Irrespective of the chosen modeling method, our primary focus here is the representation of multi-sided surfaces. These can be categorized as (i) smoothly connected macropatches composed of quadrilaterals, (ii) transfinite interpolation and (iii) control point based representations. A recent publication by Peters (2019) exhaustively surveys the variety of macropatch approaches, while a good summary of transfinite schemes can be found in Várady et al. (2011).
In the following, we will restrict our discussion to control point based methods. Concerning those that generalize the Bernstein-Bézier representation and thus reproduce Bézier boundaries, we have already reviewed several approaches in Várady et al. (2016Várady et al. ( , 2017. Recent advances include G 2 continuity conditions for toric patches (Sun and Zhu, 2018), removing the twist compatibility assumptions of GB patches (Hettinga and Kosinka, 2018), and extending the GB patch to concave, and then curved domains (Salvi and Várady, 2018;Várady et al., 2020). According to our best knowledge, the idea of creating curved domains for control point based surfaces to mimic the 3D boundary configuration first appeared in our 2020 paper.
There are also control point based, multi-sided patches that are bounded by B-spline curves. The approach in Pla-Garcia et al. (2006) is somewhat similar to ours, as they map B-spline basis functions to a special domain structure based on a symmetric planar curve network with fixed knot vectors. A recent publication by Hettinga and Kosinka (2020) describes a hole filling approach, created at the irregular points of a Catmull-Clark subdivision. This represents a solution to an important problem, using modified (cubic and quintic) blending functions in the cross direction, extending those of the convex GB patch (Várady et al., 2016). The B-spline boundaries always have two segments; the ribbons are directly deduced from the control polyhedron and ensure G 2 continuity with the surrounding surface elements.
The GBS patch proposed in this paper lifts most (implicit or explicit) assumptions of previous representations. It does not require uniform knots or symmetric arrangements, and combines B-spline ribbons with arbitrary knot vectors, degrees, and geometric continuity. In contrast with patches with convex domains, it is capable of handling complex concave boundaries and hole loops with high curvature variations.

Surface equations
The GBS patch is the normalized sum of n ribbon interpolants prescribing positional and cross-derivative constraints. The interpolants are tensor product surfaces, defined as B-splines along the boundary curves, and Bernstein polynomials in the cross direction. As was mentioned in the Introduction, ribbons can be one of two kinds: open ribbons bounded by convex corners, which are given as B-splines with interpolating end conditions; or closed ribbons for interior hole loops, which are given as periodic B-splines. Let the (non-uniform) B-spline of the interpolant at side i have degree d i and knot vector ξ i = [ξ 0 , ξ 1 , . . . ξ k i ]. The interpolants are defined by (p i + 1) × (m i + 1) control points, defining boundary positions and cross-derivatives for G m i -continuous 1 connections: where s i , h i ∈ [0, 1] are parameters along the boundary, and in the cross direction, respectively. These are local to the interpolant for side i over a curved, possibly multiply connected planar region Ω ⊂ R 2 . P i j,k refers to the j-th control point in the k-th row, which gets multiplied by a product of the corresponding B-spline basis function N d i ,ξ i j and the (Note that half of the control points of the tensor product surface, with k > m i , is not used, similarly to the half-Bézier interpolants of Salvi and Várady, 2018.) The ribbon for side i needs to vanish along all other sides; to this end we multiply the basis functions by rational weights µ i j,k , that are constant 1 for closed ribbons, while for open ribbons they depend on a subset of local h-coordinates h i = {h i−1 , h i , h i+1 } (using cyclic indexing, with 1 following n) as follows: Similarly to CD Bézier patches, these rational weights ensure that the effect of each ribbon is localized to its own side, and vanishes along its neighbors. (When p i − m i ≤ j ≤ m i , which occurs e.g. in the case of a two-segment boundary with G 2 continuity, µ i j,k takes on the product of the related terms, as in Hettinga and Kosinka, 2020.) The patch is then defined over the domain Ω, with global parameters (u, v), as a sum over the ribbons: normalized by the sum of the basis functions We see that GBS patches are very similar to CD Bézier patches -the key differences are in the use of B-spline basis functions along the boundary, and also in the local parameterization, which will be discussed in 4.2. Note that GBS patches interpolate the ribbons along their boundaries, and this property is naturally preserved when a new knot is inserted. The question arises, whether this will change the interior of the patch. It is easy to prove that if a knot gets inserted in an internal knot interval, i.e., not in the first and not in the last one, the patch remains algebraically identical. At the same time, a new knot in one of the end intervals will change the local parameterization as described in 4.2, thus it will also change the shape of the interior to a minimal extent.

Domain and parameterization
In this section we describe our method for curved domain generation and parameterization.

Curved domain generation
The curved domains of GBS patches are generated with the method described in Várady et al. (2020). In essence, we compute a development of the boundary curves into the tangent planes of the interpolants (Sternberg, 2012), meaning we flatten the curves while preserving their shape within the surface (i.e., their geodesic curvature). We first discretize the boundary curves into polylines, and compute the tangent planes of the B-spline interpolants at the sample points. The geodesic curvature at a point is approximated by simply projecting its neighborhood into the tangent plane and measuring the resulting planar angle. The polyline is then developed into the plane by placing the segments isometrically, in sequence, with lengths (l 1 , . . . ,l N ) and angles (θ 1 , . . . ,θ N ).
The resulting planar curve (p 0 ,p 1 , . . .p N ) might not form a closed loop, which is corrected via a two-step procedure. First the angles are rescaled to be compatible with forming a planar polygon: after which the displacement error between the first and last points is distributed along the curve: In the presence of interior boundary loops, we first compute the development of the outermost perimeter loop and evaluate its corresponding GBS surface. The planar developments of interior loops are then placed and oriented with the help of their projections onto the perimeter surface. We refer to Várady et al. (2020) for more details. Note that interpolants might be oriented in such a manner that interior boundaries become geodesic curves and develop into straight lines before the loop is closed, which would result in highly distorted domains. In such cases, we project curves onto their smoothed osculating planes, computed with the method of Huang and Ju (2016).

Local parameterization
Given a curved domain Ω ⊂ R 2 , our next goal is to define local side parameters (s i , h i ) : Ω → [0, 1] 2 . Recall that each side is composed of several segments, corresponding to knot spans of the B-spline ribbons. Naïvely, the local h i -parameter lines of ribbon i would be spaced linearly along the entirety of sides (i − 1) and (i + 1), and similarly the s i parameters would be constant along the adjacent sides -see Figure 2a. However, as was observed in Várady et al., 2020 (Section 6), stretching a ribbon with Bézier cross-derivatives along the entire length of adjacent sides might be inappropriate, as the effect of even relatively small ribbons would extend over the entire surface, leading to shape problems. Such issues were avoided by splitting longer boundary curves into shorter (Bézier) segments, which restricted the effect of adjacent ribbons. With GBS patches, similar shape problems might arise, unless the local parameterizations are adapted to contain the influence of the ribbons. Thus, for GBS patches the h i -parameter lines increase linearly only along the directly adjacent segments of sides (i − 1) and (i + 1), and the s i -coordinates are similarly constrained, as shown in Figure 2d. Another difference compared to (Várady et al., 2020) is that in that work, local parameters were computed for each Bézier segment independently, using generalized barycentric coordinates, resulting in overlapping parameterizations, see Figure 2b. In contrast, for GBS patches, the h-coordinate is common for each segment, while the s-coordinate lines of different segments are localized within the knot lines, as shown in Figure 2e.
We choose to define both s i and h i as harmonic functions, satisfying the minimal necessary set of boundary conditions. 2 Unlike the method based on generalized barycentric coordinates, this approach can be generalized easily to the case of periodic interior loops as well. The boundary conditions are summarized as follows: • The coordinate function s i (for side i) should be constant 0 along the neighboring segment of side i − 1 and 1 for side i + 1, while along side i it should interpolate between 0 and 1, in accordance with the knot vector. 3 Note that the value of s i is unspecified along the remaining parts of the domain boundary. • The coordinate function h i should be constant 0 along side i, interpolate between 0 and 1 along the adjacent segments of sides i − 1 and i + 1; and it should take on the value 1 everywhere else along the boundary.
We look for harmonic functions, which are constrained minimizers of the Dirichlet energy: where f = s i or h i , and b f encodes the aforementioned Dirichlet boundary conditions for the constrained subset of the boundary D f ⊂ ∂Ω. Harmonic functions are preferred, as they are guaranteed to be C ∞ -continuous in the interior, and take their minimal and maximal values on the boundary (Hélein and Wood, 2008). The continuous problem (7) can be discretized using piecewise-linear functions over a triangle mesh M = (V, E, F ) to arrive at a convex quadratic energy of the vertex values where w i j are the cotangent weights (Pinkall and Polthier, 1993). The solution of this problem then satisfies the discrete Laplace equation where L is the cotangent Laplacian matrix. When f = s i , the solution of (9) will also satisfy Neumann natural boundary conditions along unconstrained segments, making the s-coordinate lines perpendicular to the boundary (Stein et al., 2018).

5
Note, that although we compute the coordinates over a tessellation, solutions of the Laplace equation can also be approximated without discretizing the domain interior (see e.g. Sawhney and Crane (2020) and the references therein) -thus, in principle, our patches can be evaluated at arbitrary points of the domain.
Periodic interior loops. Multiply connected domains might have interior loops with periodic parameterization, for which the previously described method must be modified. Computing the h-coordinate is relatively straightforward: the boundary conditions are constant 0 along the current inner loop and 1 everywhere else. For the s-coordinate, we need to compute a function interpolating discontinuous boundary conditions along the loop.
Let us compute a piecewise-linear polyline C connecting the s = 0 point on the periodic loop with the outer perimeter loop, which will serve as a virtual cut in the domain. In our implementation, we trace the gradient field of the h-coordinate, but other choices are possible as well. We will require that the function values jump by 1 when crossing the cut. We modify the discretized optimization problem (8) accordingly as follows: where is the required jump while crossing the edge e i j . The solution of this problem satisfies the discrete Poisson equation where the right-hand side ρ i = − e i j w i j ρ i j encodes the jump conditions along the cut. The solution of this linear system might take values outside [0, 1], in which case the function values are simply shifted by ±1. It has been shown in the context of similar problems in seamless mesh parameterization (Aigerman and Lipman, 2015;Bright et al., 2017) and quad meshing (Tong et al., 2006), that the solution of (12) depends only on the homotopy type of the cut within the domain. Compare the overlapping multi-loop parameterization of Várady et al. (2020) shown in Figure 2c, with our new periodic scheme shown in Figure 2f.

Editing the interior
In our previous paper, we have raised the issue of adding interior control points to the surface formula, when the internal curvature distribution did not entirely match the designer's expectation and further adjustments are needed. For curved domain Bézier patches the relative power of the ribbons decreases when high-degree Bézier boundaries are used, however with B-splines, complexity does not affect the degree, and this phenomenon is less frequent. If we wish to modify the shape of the interior, we can add further control points Q i , i = 1, . . . , q, multiplied by appropriate blending functions B Int i (u, v). "Appropriate" means that the interior control points should not affect the boundary constraints, just the transitions between the ribbons. Both the basic equation (Eq. 3) and the sum of the blending functions (Eq. 4) need to be modified: Here we show two interior shape editing constructions we found useful; others can be defined in a similar fashion. In the first case, assume that the ribbons are defined by two rows of control points and we wish to retain G 1 continuity. We add ribbon control points {P Rib i, j }, corresponding to Q i above, associated with the i-th side. The basis functions 6 (a) Without ribbon control points (b) With ribbon control points on 3 sides corresponding to B Int i are defined as B-spline basis functions along the boundary, and quartic Bernstein functions in the cross-direction: Possible values of the index j range between 2 and p i − 2. The effect of these points will vanish on the adjacent sides due to the B-spline basis term, and on the distant sides due to the quartic Bernstein functions. A simple default is to set P Rib i, j = P i 1, j + (P i 1, j − P i 0, j ). This construction can be applied for arbitrary sides, and generalized for ribbons with more than two rows, retaining higher degree geometric continuity. A GBS patch -with and without ribbon control points -is shown in Figure 3. In our second example we introduce a central control point P C , where the distance parameters h i are used in the blending function. For example, in the G 1 case guarantees that the effect of P C vanishes on all sides. Here ζ is a user-defined constant, controlling the strength of the control point. We can make the following observations: 1. Higher degree continuity can be achieved by increasing the exponent in Eq. 16. 7 2. Our practical experience shows that for interior editing it may be preferable to interpolate internal surface points. For example, it is easy to compute a central control point by repositioning a selected surface point at (u 0 , v 0 ) into a given positionŜ 0 , by solving the equation A GBS patch after lifting a surface point to modify the fullness of the patch is shown in Figure 4.

Discussion and case studies
In this section we discuss some interesting features of GBS patches, and show related test examples. We do not deal with the editing operations, i.e., how the B-spline ribbons were actually defined, and our forthcoming analysis mainly relates to the capabilities of GBS patches using an already specified input. 6.1. Convex vs. curved domains In our previous paper (Várady et al., 2020) we have shown several examples concerning the difficulties of modeling complex multi-sided patches with concave boundaries. Classical approaches often fail, and the patches may wrinkle or the shape may fold under itself. As an illustration, we have picked the Charrot-Gregory patch (Charrot and Gregory, 1984) to be compared with the GBS patch; see curvature maps in Figure 5.

A complex vertex blend
This is an interesting example, where two setback vertex blends (Várady and Hoffmann, 1998) need to be combined. The blending parameters are set in such a way that the short edge in the middle degenerates, and the two six-sided vertex blends are represented as a single eight-sided patch with strong concave boundary curves. The edge configuration and the GBS patch with B-spline ribbons and curvature map are shown in Figure 6.

Polyhedral design
GBS patches are well-suited to define a collection of patches, defined indirectly by a control polyhedron and an associated free-form curve network. An example is shown in Figure 7.

G 2 connection
While curvenet-based models typically utilize G 1 ribbons, it may occur that the adjacent patches need to be connected with G 2 continuity. In these cases three-layer ribbons are used; once these satisfy G 2 continuity, the related multi-sided patches will also connect by G 2 . An example with an 8-sided and a 4-sided patch, meeting in along a strongly curved boundary, is shown in Figure 8.

Multi-sided patches at extraordinary vertices
One important utilization of multi-sided patches is to fill holes at the extraordinary vertices of quadrilateral meshes, for example, in Catmull-Clark subdivision. The quality of these patches is particularly important, and various solutions have been published to ensure high surface qualities (see constrained quadrilateral patches in the work of Peters (2019), and the recent multi-sided B-spline scheme by Hettinga and Kosinka, 2020). The key idea here is to borrow positional and cross-derivative information from the surrounding regular patches and insert a multi-sided patch that 9 (a) The beam6 model (b) The monk7 model Figure 9: GBS patches as G 2 caps in subdivision surfaces with mean curvature map and isophote lines near the center. has compatible ribbons. GBS patches with three-layer ribbons also provide straightforward G 2 solutions. We have picked examples from the test suite of Peters (2020) and examined the curvature map and the distribution of the isophotes around the extraordinary points. The surfaces beam6 and monk7 are shown in Figure 9.

Surfaces with hole loops
We show two test parts that illustrate how multi-loop modeling works. The first is a sheet metal part -first without the interior hole ribbon, then enforcing the ribbon constraint, see Figure 10. The second example is a plastic bottle (Figure 11), showing the original surface and the modified one with the interior loop. The cross-derivatives for the side ribbons and the hole loop must be orthogonal to the symmetry plane of the bottle in order to smoothly join the two halves.

Sheet metal parts
Two surface models defined by concept car sketches are shown here. The first model is a surface over the hood; the patch with contouring and curvature map is shown in Figure 12. The second surface represents the back part of a sports car; the patch and its curvature map is shown in Figure 13. In the latter case, the control structure contains two ribbons with three layers for additional shape control.

Special cases
This model has three surfaces ( Figure 14). The largest is a four-sided GBS patch with a long, complex boundary, being connected to an extruded surface, terminated by a two-sided GBS patch, which naturally fits into the concept of our representation (Figure 14b) . Observe, that at the common vertex of the blue and yellow ribbons of the large patch, the ribbons are incompatible. The entire patch interior is nicely defined by our representation; however, the patch is singular at the common vertex. This incompatibility was necessary to provide a smooth connection to the extruded surface, as contouring shows.

Conclusion and future work
We have proposed a new multi-sided surface representation that interpolates complex convex-concave B-spline boundary curves with associated cross-derivatives and periodic B-splines in the interior of the patches. The basic elements of the formulation include curved domain generation, harmonic parameterization for multi-connected domains and a special combination of B-spline and Bernstein basis functions with rational terms. We have shown some interesting shapes being represented by single GBS patches. We believe that these can hardly be produced by means of former multi-sided schemes based on convex polygonal domains.
At the same time, plenty of problems remained for future research. One particular issue is the weight deficiency in the interior of the patches. We have shown that adding new layers and ribbon control points can help, but for perfect shape control, a uniform structure of interior control points would be required that supplements boundary data. Another problem is that harmonic parameterization produces somewhat unevenly distributed h parameter lines for long, strongly curved boundaries. Resolving this would ensure that the cross-derivatives affect the interior always to the same extent. Finally, we remark that designing appropriate B-spline ribbons may become difficult in certain situations, thus interactive CAD operations that support the designers' work would facilitate the wider use of GBS patches.