Typical curve with G1 constraints for curve completion

This paper presents a novel algorithm for planar G1 interpolation using typical curves with monotonic curvature. The G1 interpolation problem is converted into a system of nonlinear equations and sufficient conditions are provided to check whether there is a solution. The proposed algorithm was applied to a curve completion task. The main advantages of the proposed method are its simple construction, compatibility with NURBS, and monotonic curvature.


Introduction
G1 interpolation is an essential problem in many applications such as path planning and curve completion. The goal is to find a transition curve that matches the positions and associated unit tangent vectors at two given endpoints [1]. Curve completion aims to find a pleasing contour to fill in an object boundary that is partially occluded [2,3]. This process differs significantly for humans and computers. For example, the visual system of humans can automatically complete missing parts of curves ( Fig. 1), which is called visual curve completion [4][5][6]. For computers, it is necessary to identify a fairing curve matching "boundary conditions" among potential transition curves, which is called curve completion.
In computer-aided design (CAD) and computeraided geometric design, the measurement of the fairness of curves typically depends on the curvature distribution. Farin and Sapidis [7] assumed that a fairing curve should have relatively few monotonic curvature variation segments. By virtue of linear curvature variation and the ability to minimize changes in the total curvature, the Euler spiral is considered to be the most pleasing curve for shape completion [2], so it has been widely used in various applications such as path planning and curve completion [2,[8][9][10][11]. However, the Euler spiral is defined in the form of a transcendental function that is computationally intensive and inefficient. Furthermore, the Euler spiral is not compatible with the NURBS methods, which are the standard methods in existing CAD software [12].
To overcome these shortcomings, some methods using polynomial curves or arc splines to approximate Euler curves have been proposed [3,[13][14][15][16]. Another solution is to use a polynomial curve as a design object directly. In highway design, a transition curve is typically defined between a line and a circular curve or between two circular curves. Walton et al. [17][18][19] used planar cubic Bézier spirals to construct transition curves with G2 continuity for different boundary constraints and the results of refs. [17,19] were extended in refs. [20,21]. Traditional methods for designing Bézier curves with monotonic curvature focused on adjusting the positions of control points until Mineur proposed the concept of a typical curve [22], which indicates a special subset of planar Bézier curves with control edges that maintain a specific geometric constraint. Sánchez-Reyes [23] pointed out that these constrained typical curves belong to the family of offset rational sinusoidal spirals and can also be expressed as a subset of rational Bézier curves called p-Bézier curves [24]. Recently, Wang et al. [25] proposed a geometric proof of the necessary and sufficient conditions for the curvature monotonicity of typical curves. In 2006, Farin [26] extended typical curves to 3D class-A Bézier curves, whose control edges were generated from a class-A matrix and initial control edge. However, the conditions of the class-A matrix have been proven to be incomplete [27] and several counterexamples have been reported [28]. Yoshida et al. [29] proposed an interactive method to generate general class-A Bézier curves by perturbing the elements of a planar typical class-A matrix. Furthermore, typical curves converge to logarithmic spiral segments as the degree is elevated and the drawable regions of typical curves for boundary constraints under different degrees can be obtained.
In recent years, geometric continuity has also been discussed based on trigonometric Bézier curves with shape parameters [30][31][32][33]. Generalized trigonometric Bézier curves are flexible and can achieve G2 continuity. However, these methods are more complex than polynomials and have not provided a way to obtain a monotonic variation of curvature. This study analyzed the sufficient conditions for the G1 interpolation problem using typical curves and developed a novel algorithm to find a fairing solution, which is applied to curve completion. For a given set of constraints, multiple solutions may exist, so a suitable criterion must be defined to find the optimal solution automatically. Furthermore, a designer can modify data manually to obtain a fair solution to fit different situations. The main contributions of this work can be summarized as follows. The sufficient conditions for G1 interpolation based on typical curves are provided and curve completion using Bézier curves with monotonic curvature is realized.
The remainder of this paper is organized as follows. Typical curves Section introduces typical curves comprehensively. Methods Section discusses the G1 interpolation problem, provides sufficient conditions for a typical curve solution, and presents a novel algorithm to solve G1 boundary constraints. In Results and discussion Section, some examples are presented to demonstrate the practicability and superiority of the proposed algorithm. The method is then applied in the application of curve completion. Finally, this paper is concluded in Conclusions Section.

Typical curves
A planar Bézier curve of degree k can be expressed as where b i are the two-dimensional control points and B i, be the control edge. Then, a typical curve is obtained with control edges satisfying where s is a positive scale factor and R θ is a secondorder rotation matrix with a rotation angle θ∈½− π 2 ; π 2 . A k-degree typical curve exhibits monotonic curvature variation if and only if When θ∈½0; π 2 , it is assumed that the control edge rotates counterclockwise ( Fig. 2(a)) and the relative curvature is assumed to be positive ( Fig. 2(b)). In this case, if the scale factor s > 1, then the curve P(t) has a monotonically decreasing curvature. If 0 < s < 1, then P(t) has a monotonically increasing curvature. When θ∈½− π 2 ; 0, it is assumed that the control edge rotates clockwise ( Fig. 2(c)) and the relative curvature is negative ( Fig.   2(d)), so the monotonicity of curvature is opposite to θ∈½ 0; π 2 . As a result of the construction process of control edge vectors, the generated curve can only be 'C' shaped rather than 'S' shaped.

Methods
This section presents the G1 interpolation algorithm based on typical curves. First, the G1 interpolation problem based on typical curves is introduced briefly. The given endpoint-orientation pairs are then divided into two cases, and a novel algorithm for constructing a typical curve under the given boundary constraints is proposed.

G1 Hermite interpolation problem
In many interpolation problems, not only positions, but also the corresponding derivative values or higher-order derivatives at endpoints are required to be equal, which is called Hermite interpolation. If the interpolation curve matches the given unit tangent vectors at the two endpoints, it is a two-point G1 Hermite interpolation. Additionally, if the given curvatures at the two endpoints are also matched, it is called two-point G2 Hermite interpolation [1]. Although G2 interpolation has better continuity, it also introduces stronger constraints for the interpolation curve and the high-order continuity at endpoints is sometimes difficult to guarantee. In some applications, G1 interpolation is preferable to G2 interpolation when considering cost effectiveness, such as curve completion or path planning using Euler spirals [1]. Consider two different points P A = (x A , y A ) T and P B = (x B , y B ) T with the corresponding unit tangent vectors T A = (cosθ A , sinθ A ) T and T B = (cosθ B , sinθ B ) T , where the three vectors T A , T B , and P A P B = P B − P A are not all parallel. G1 Hermite interpolation is used to find a fairing curve that joins P A and P B on the condition that the tangent vectors at the endpoints match T A and T B , respectively.
Walton et al. [1,34] proposed an improved Euler spiral algorithm for shape completion. However, the Euler spiral is not compatible with NURBS. Therefore, the paper proposes a more intuitive algorithm for the G1 interpolation problem based on a typical curve, which is easier to construct.

G1 interpolation problem based on a typical curve Proposition 1
Consider the positions of two endpoints P A = (x A , y A ) T and P B = (x B , y B ) T , and their associated unit tangent vectors T A = (cosθ A , sinθ A ) T and T B = (cosθ B , sinθ B ) T , where θ A and θ B are orientation angles with multiple values. When θ A or θ B are increased or reduced by 2m ⋅ π(m ∈ Z), the unit tangent vector remains  unchanged. The sufficient condition for the existence of a k-degree typical curve matching the boundary constraints is that the following equation has real number solutions satisfying s ⋅ cos θ ≥ 1 (s ≥ 1) or s ≤ cos θ (0 < s < 1): where ‖V 0 ‖ and s are free variables with positive values and θ = (θ B − θ A )/(k − 1).

Proof
When the degree k is fixed, θ = (θ B − θ A )/(k − 1) is the rotation angle of the typical curve. If there are real number solutions such that s ⋅ cos θ ≥ 1 (s ≥ 1) or s ≤ cos θ (0 < s < 1), which are the sufficient and necessary conditions for the monotonic curvature of a typical curve, the curvature of the Bézier curve constructed using this method must be monotonous. Now we are going to analyze the sufficient conditions for the existence of solutions for Eq. (4). Considering the degree k as a free variable in Eq. (4), there are three free variables k ∈ Z + , s > 0, and Because Bézier curves are invariant under affine transformations, one can place the first endpoint P A = (x A , y A ) T at the origin and the second endpoint P B = (x B , y B ) T at the positive half of the x axis using a pretransformation, which leads to x' A = 0, y' A = 0, x' B > 0 and y' B = 0 (Fig. 3). Then, Eq. (5) is converted into The superscript ' indicates a new variable obtained from the pre-transformation. The pre-transformation can be used to construct a new local coordinate system. Equation (6) is a polynomial of degree k and one can infer the number of positive roots from the Descartes rule of signs [35]. The number of sign changes in the coefficients is equal to the number of sign changes in sin[θ' A + (i − 1) ⋅ θ], i = 1, …, k. Therefore, a sufficient condition to guarantee that Eq. (6) has at least one positive solution is that sin[θ' A + (i − 1) ⋅ θ] changes signs an odd number of times, meaning that after pre-transformation, the sign of the sine angle for all control edges changes an odd number of times.

Angle sign change condition
Consider P A as the origin and P A P B as the positive direction of the x axis. In this new local coordinate system, the function sinα (θ' A ≤ α ≤ θ' B ) only changes its sign an odd number of times.

Proposition 2
If the ASC condition is satisfied, then one can always find an appropriate degree k such that the solutions of Eq. (4) or (6) satisfy Eq. (3).

Proof
Because the affine transformation does not change the solution of the curve, Eqs. (4,6) have the same solution and the ASC condition guarantees a positive solution for Eq. (4). When the degree k gradually increases in Eq. (6) and the rotation angle θ ¼ θ B −θ A k−1 → 0 cosθ → 1, there is an s k > 0 such that s k ⋅ cos θ ≥ 1 (s k ≥ 1) or s k ≤ cos θ (0 < s k < 1), meaning the solutions of Eq. (4) or (6) satisfy Eq. (3).

Corollary 1
For Eq. (4) or (6), if there is a typical curve satisfying the constraints when k = m ≥ 2, then there are also typical curve solutions when k ≥ m + 1. This means that the G1 interpolation problem has multiple solutions consisting of typical curves and the limit of the solutions increases θ → 0 with an increase in the degree k.
The ASC condition is a sufficient condition for the positive solution of Eq. (6). Even if the ASC condition is not satisfied, Eq. (4) or (6) may still have a positive root. In contrast, the ASC condition cannot guarantee solutions satisfying Eq. (3) under a fixed degree.

Two cases of G1 constraints
For convenience of expression, homogeneous coordinates will be used to represent points and vectors for the remainder of this paper. For a 2D point (x, y) T , its homogeneous coordinate can be written as (ωx, ωy, ω) T , where ω ≠ 0. This study used ω = 1. For a 2D vector (x, y) T , the homogeneous coordinate can only be expressed as (x, y, 0) T .
For an arbitrary P A = (x A , y A , 1) T , P B = (x B , y B , 1) T , T A = (cosθ A , sinθ A , 0) T , and T B = (cosθ B , sinθ B , 0) T , one can divide the G1 constraints into two cases according to the relationship between the relative positions of the three vectors T A = (cosθ A , sinθ A , 0) T , T B = (cosθ B , sinθ B , 0) T , and P A P B = (x B − x A , y B − y A , 0) T . This analysis does not consider the degenerate case in which all three vectors are parallel. Let −π ≤ φ A ≤ π denote the angle from vector T A to vector P A P B (Fig. 4). The sign of φ A indicates the direction and φ A > 0 indicates that T A is on the right side of P A P B . Similarly, −π ≤ φ B ≤ π denotes the angle from vector P A P B to vector T B .
Considering the signs of these two angles, one can divide all potential constraints into two cases as follows: Case I φ A ⋅ φ B > 0, meaning T A and T B are located on opposite sides of P A P B (Fig. 5).
Case II φ A ⋅ φ B < 0, meaning that T A and T B are located on the same side of P A P B (Fig. 6).
Because θ A and θ B can be increased or decreased by 2m ⋅ π (m ∈ Z) while keeping the unit tangent vectors T A and T B unchanged, the rotation angle θ = (θ B − θ A )/(k − 1) can be either positive or negative, and the labeling of the given endpoint-orientation pairs can be swapped. Based on these features, there are different constructions for typical curves under the given conditions. This paper presents a practical method for choosing a suitable curve based on the position relationship between the unit tangent vectors T A T B and the vector P A P B .

Proposed algorithm
Given two endpoint-orientation pairs, the G1 interpolation algorithm is used to find the optimal solution to Eq. (4), where the degree k ≥ 2 is an integral and the rotation angle θ may have multiple values. s > 0 and ‖V 0 ‖ > 0 are the variables to be determined. The value of k can be increased from two to a set maximum value k max . For a fixed k, one must determine the value of θ = (θ B − θ A )/(k − 1), which is equivalent to determining Δθ = θ B − θ A . The proposed algorithm provides a rule for determining Δθ automatically. For case I, let Δθ = φ A + φ B . For case II, let In this manner, the algorithm can find the typical curve with the lowest degree and smallest total angle between T A and T B .
Once k and Δθ are determined, one can use an optimization process to obtain s and ‖V 0 ‖ such that becomes zero. If s > 0, ‖V 0 ‖ > 0, and Eq. (3) is satisfied, then a typical curve segment can be constructed using these values. Otherwise, the steps of the iteration should be repeated by increasing the degree k. The pseudo code for the proposed algorithm is presented in Algorithm 1.
The algorithm first calculates an angle Δθ, which is determined by the position relationship between, T A , T B , and P A P B . Then, Eq. (7) is minimized to a value of zero iteratively by the optimization process. This step was implemented using the interior point method, which can find the minimum value of a constrained nonlinear multivariable function. For the obtained s > 0 and ‖V 0 ‖ > 0, the algorithm checks whether Eq. (3) is satisfied. Once the solutions satisfy Eq. (3), iteration can be stopped and a typical curve can be generated from the obtained data.
The algorithm proposed above avoids outputting multiple solutions by choosing an appropriate rotation angle θ = (θ B − θ A )/(k − 1), but this does not mean that there is only one typical curve solution under the given constraints. The main principle of the proposed algorithm is to select the minimum rotation angle Δθ = θ B − θ A under the given conditions and then choose the lowest degree k to satisfy the constraints. Additional solutions can be obtained by increasing the value of the degree k, modifying the angle value Δθ = θ B − θ A while leaving the unit tangents T A and T B unchanged, or even by exchanging the labels of P A and P B or T A and T B . These methods are discussed in detail in Results and discussion Section.

Results and discussion
This section presents some numerical experiments based on the proposed algorithm and discusses the existence of multiple solutions.

Examples of two cases of constraints
This section presents different results for two cases of constraints.

Example 1
Consider data P A = (0, 0, 1) T , θ A = 0, T A = (1, 0, 0) T , P B = (50, 50, 1) T , θ B ¼ 2π 3 , and Fig. 7). It is easy to verify that this is an instance of case I that satisfies the ASC condition. In the proposed algorithm, Δθ = θ B − θ A is chosen as Δθ ¼ 2π 3 and a solution is obtained when the degree k = 5 (Fig. 7). Figure 7(a) presents the boundary constraints and resulting curve, where P A and T A are marked in blue, and P B and T B are marked in green. Figure 7(b) presents an increasing curvature plot. The data obtained from the typical curve are presented in Table 1.

Discussion of multiple solutions
Under the given constraints, there are three main reasons for generating multiple solutions: an increase in the  degree k, different chosen values of Δθ, and the exchange of labels of two endpoint-orientation pairs. Examples of these cases are discussed below.
Increasing the degree k According to Corollary 1, multiple solutions with different degrees exist.  (Table 3).

Different values of Δθ
The proposed algorithm provides a criterion for the determination Δθ. This is a practical rule. In fact, one can add 2m ⋅ π(m ∈ Z) to the tangent vector angle while keeping the tangent direction unchanged. Therefore, one can obtain new solutions that satisfy G1 interpolation by adding 2m ⋅ π(m ∈ Z) to the Δθ selected by the previously proposed algorithm. Figure 10 presents two typical curves with different Δθ, both meeting the given condi-  (Table 4).

Exchange of two endpoint-orientation pairs
The construction process of typical curves depends on the sequence of endpoints. If one exchanges the labels of endpoint-orientation pairs, one can obtain different solutions that meet the given conditions. As shown in Fig. 11, the red and carmine curves are two different typical curves that match the same G1 constraint. The difference between the two results is that the red curve takes the blue endpoint-orientation as the starting point and the green endpoint-orientation as the ending point. This , and Δθ ¼ − 7π 6 . The specific data for these two typical curves are presented in Table 5.

Curve completion
Curve completion is a geometric continuation of the boundaries of objects that are temporarily interrupted by occlusion [2]. Curve completion requires a visually pleasing shape and is similar to G1 interpolation, but more expansive. The principles for finding an optimal completion curve include minimizing the total curvature square [36] and minimizing the total curvature variation [2]. In this section, these rules are simplified to realize curve completion by achieving monotonic curvature. Two sets of experimental results are presented to demonstrate the application of the proposed method to curve completion and discuss the differences between the proposed approach and other interpolation methods. Figure 12(a) presents a banana that is partially covered by a piece of tape. Figure 12(b) presents the completed curve drawn by the proposed method and Fig. 12(c) presents a local zoomed-in view of the completed curve. Figure 12(d) presents the increasing curvature. Figure 13(a) presents two overlapping leaves, one of which is partially occluded by the other. The proposed method completes the occluded boundary using a typical curve, as shown in Fig. 13(b). The corresponding monotonic curvature is presented in Fig. 13(c).  For the curve completion problem, many studies and methods have been discussed, including methods based on Euler spirals [2,8,10], Euler arc splines [3], and cubic B-spline curves [37], but these methods have their own disadvantages. The Euler spiral is expressed by transcendent equations, which may lead to low computational efficiency. Another shortcoming of the Euler spiral is incompatibility with NURBS. A typical curve is a polynomial curve in the Bézier form and the G1 interpolation problem can be transformed into a simple system of nonlinear equations.
An Euler arc spline [3], which consists of several arcs with the same arc length, is considered to be an extension of the Euler spiral. Although it can be represented by NURBS, the continuity at the internal connecting points is only G1 continuous, rather than continuous curvature. In contrast, the curvature of a typical curve is continuous and monotonous over the entire curve.
For the cubic B-spline curve with monotonic curvature mentioned in ref. [37], with an increase in the number of control vectors, the B-spline curve gradually converges to a straight line as the parameter trends toward one. This feature makes the bending curvature inadequate and leads to poor control flexibility. Furthermore, the algorithm for curve completion using a monotonically cubic B-spline curve is a brutal algorithm and the cases of boundary conditions have not been fully elucidated. In the proposed method, the bending curvature can be controlled well based on the rotation angle θ or degree k, and the construction of a typical curve is simpler than the cubic B-spline curve constructed in ref. [37]. Additionally, the proposed algorithm can handle more constrained cases, which makes it more widely applicable.

Conclusions
This paper proposed a novel algorithm for G1 interpolation based on typical curves. G1 interpolation was converted into a system of nonlinear equations and a sufficient condition was provided to determine whether there is a typical curve solution with monotonic curvature in advance.