**4. Algorithms for (Rational) Bernstein Polynomials**

This section contains algorithms and procedures for Bernstein polynomials that make use of the properties presented in Section 2. These functions include: evaluating bounds, using the convex hull property to quickly find conservative bounds; evaluating extrema, through an iterative procedure that computes a solution within a desired tolerance; minimum spatial distance, applying a similar iterative procedure to find the minimum spatial distance between two Bernstein polynomials; and collision detection, which quickly determines whether a collision may exist or does not exist.

#### *4.1. Evaluating Bounds*

Property A1 allows one to quickly establish conservative bounds on the Bernstein polynomial. For example, given the 2D Bernstein polynomial, see Equation (5), with coefficients given by

$$\mathbf{P}\_5 = \begin{bmatrix} 0 & 1 & 2 & 3 & 4 & 5 \\ 5 & 0 & 2 & 5 & 7 & 5 \end{bmatrix} \text{ \textquotedblleft}$$

lower and upper bounds *C*min and *C*max satisfying *C*min ≤ *C*(*t*) ≤ *C*max, ∀*t* ∈ [*t*0, *tf* ] can be derived using Equation (A1). Figure 18 exhibits the Bernstein polynomial (solid blue line) given the above coefficients (orange dots connected with dashes). The most conservative estimate of the minimum and maximum *Y* values of the Bernstein polynomial is given by the coefficients with the lowest and highest *Y* values, respectively. The lower bound is 0

and the upper bound is 7. As mentioned in Property A6 and Equation (A8), the Bernstein coefficients converge to the curve as the polynomial is degree elevated. This fact can be used to derive tighter bounds. A degree elevation of 20 results in a lower bound of 1.93 and an upper bound of 5.89. This is a closer estimate of the actual minimum and maximum, 2.26 and 5.70, respectively (see red dotted lines and Section 4.2). Figure 18 also illustrates degree elevations of 5, 10, and 15. Since the degree elevation matrix, see Equation (A7), is independent of the Bernstein coefficients, a database of elevation matrices can be computed ahead of time to produce tight estimates of the bounds at a low computational cost.

**Figure 18.** Bounds for Bernstein polynomials. The solid blue line is the Bernstein polynomial, the dashed lines connect the coefficients of each different order, the black dotted line represents the convex hull of the 5th degree Bernstein polynomial, and the red dotted lines represent the actual extrema.

#### *4.2. Evaluating Extrema*

The extrema of a Bernstein polynomial are calculated using an iterative procedure similar to the one proposed in [44]. This is done by recursively splitting the curve and using the Convex Hull (Property A1) to obtain an estimate within some desired tolerance. Algorithm 1 outlines the process for determining the maximum of a Bernstein polynomial and can easily be modified to determine the minimum of a Bernstein polynomial.

The inputs to Algorithm 1 are the Bernstein polynomial's coefficients, P = {**P***n*}, **P***<sup>n</sup>* = [*P*0,*n*, ... , *Pn*,*n*], an arbitrarily large *negative* global lower bound, *α*, and a desired tolerance, . Note that in order to compute a reliable maximum, *α* ≤ min (*P*). In practice, *α* is set to the lowest possible value that can be reliably represented in the computer.

Line 1 finds the maximum of the two endpoints of the Bernstein polynomial, where *n* is the degree of the polynomial. This makes use of the End Points (Property A2). Next, we determine the upper bound by simply finding the maximum of P. The *if* statement on line 3 determines whether the global lower bound should be replaced with the current lower bound. The next *if* statement, line 6, will prune the current set of Bernstein coefficients. This is valid because *α* always provides a lower bound on the global maximum. If the upper bound of any subset is below *α*, then we know that it is impossible for any point on that subset to be the global maximum. The final *if* statement, line 9, determines whether the difference between the upper and lower bounds is within the desired tolerance and returns the global minimum bound *α* if the tolerance is met.

The *else* statement, starting on line 11, splits the curve and then recursively calls Algorithm 1 again. The function split() utilizes the de Casteljau algorithm (Property A5). One of two different splitting points, *tdiv*, can be employed. The first option simply splits the curve in half, i.e., *tdiv* <sup>=</sup> *<sup>t</sup>*<sup>0</sup> <sup>+</sup> *tf* <sup>−</sup>*t*<sup>0</sup> <sup>2</sup> . The second option uses the index of the largest valued coefficient, *iub* <sup>=</sup> argmax(P), to determine the splitting point, i.e., *tdiv* <sup>=</sup> *<sup>t</sup>*<sup>0</sup> + (*tf* <sup>−</sup> *<sup>t</sup>*0)*iub n* .

Algorithm 1 (and its converse) is employed to find the minimum and maximum of the 5th degree Bernstein polynomial depicted in Figure 18 (red lines). The execution time to compute the minimum is 320 μs on a Lenovo ThinkPad laptop using an Intel Core i7-8550U, 1.80 GHz CPU. The implementation can be found in [40].

**Input:**P, *α*, **<sup>1</sup>** *Plb* = max{P[0],P[*n*]} **<sup>2</sup>** *Pub* = max(P) **<sup>3</sup> if** *Plb* > *α* **then <sup>4</sup>** *α* = *Plb* **<sup>5</sup> end <sup>6</sup> if** *α* > P*ub* **then <sup>7</sup> return** *α* **<sup>8</sup> end <sup>9</sup> if** *Pub* − *Plb* < **then <sup>10</sup> return** *α* **<sup>11</sup> else <sup>12</sup>** <sup>P</sup> *<sup>A</sup>*,P*<sup>B</sup>* <sup>=</sup> split(P) **<sup>13</sup>** *<sup>α</sup><sup>A</sup>* <sup>=</sup> Algorithm 1(<sup>P</sup> *<sup>A</sup>*, *<sup>α</sup>*, ) **<sup>14</sup>** *<sup>α</sup><sup>B</sup>* <sup>=</sup> Algorithm 1(P*B*, *<sup>α</sup>*, ) **<sup>15</sup>** *α* = max(*αA*, *αB*) **<sup>16</sup> end <sup>17</sup> return** *α*

#### **Algorithm 1:** Evaluating Maximum Value of a Bernstein Polynomial

#### *4.3. Minimum Spatial Distance*

The minimum spatial distance between two Bernstein polynomials can be computed using the method outlined in [44]. This is done by exploiting the Convex Hull property (Property A1), the End Point Values property (Property A2), the de Casteljau Algorithm (Property A5), and the Gilbert-Johnson-Keerthi (GJK) algorithm [45]. The latter is widely used in computer graphics to compute the minimum distance between convex shapes.

The algorithm for the minimum spatial distance between two Bernstein polynomials is shown in Algorithm 2. The first two inputs to the function are the sets of Bernstein coefficients, P = {**P***m*} and Q = {**Q***n*}, which define the two Bernstein polynomials in question. The last two inputs are the global upper bound on the minimum distance, *α*, and a desired tolerance .

The upper\_bound() function on line 1 finds the furthest distance between the end points of the two Bernstein polynomials, i.e., lower\_bound(P, Q) = max{P[0] − Q[0],P[0] −Q[*n*],P[*m*] − Q[0],P[*m*] − Q[*n*]} where *m* and *n* are the degrees of the polynomials represented by P and Q, respectively. This is a valid upper bound on the minimum distance between the two polynomials due to End Point Values (Property A2).

The lower\_bound() function on line 2 finds the lower bound on the distance between the two polynomials by using the GJK algorithm. This is a valid lower bound because of Property A1, Convex Hull. The next three *if* statements on lines 3, 6, and 9 are very similar to those seen in Algorithm 1. Line 3 updates the global upper bound *α* if the current upper bound is smaller. Line 6 prunes the current iteration since it is impossible the current lower bound, *lower*, to be the minimum distance if it is larger than the global upper bound. Line 9 returns *α* if the desired tolerance is met.

The lines within the *else* statement split the Bernstein polynomials defined by P and Q and recursively call Algorithm 2. Like in Algorithm 1, the first option for splitting would be to simply split at the halfway point. The second option for splitting the curves is outlined in [44] and uses the location at which the minimum distance occurs to choose the splitting point. Figure 19a illustrates the minimum distance between several different Bernstein polynomials. The code to generate this plot can be found in [40]. The execution time to compute the minimum spatial distance is 3.29 ms on a Lenovo ThinkPad laptop using an Intel Core i7-8550U, 1.80 GHz CPU.

**Input:**P, Q, *α*, **<sup>1</sup>** *upper* = upper\_bound(P, Q) **<sup>2</sup>** *lower* = lower\_bound(P, Q) **<sup>3</sup> if** *upper* < *α* **then <sup>4</sup>** *α* = *upper* **<sup>5</sup> end <sup>6</sup> if** *α* < *lower* **then <sup>7</sup> return** *α* **<sup>8</sup> end <sup>9</sup> if** *upper* − *lower* < **then <sup>10</sup> return** *α* **<sup>11</sup> else <sup>12</sup>** <sup>P</sup> *<sup>A</sup>*,P*<sup>B</sup>* <sup>=</sup> split(P) **<sup>13</sup>** <sup>Q</sup>*A*, <sup>Q</sup>*<sup>B</sup>* <sup>=</sup> split(Q) **<sup>14</sup>** *<sup>α</sup>* <sup>=</sup> min(*α*, Algorithm <sup>2</sup> (<sup>P</sup> *<sup>A</sup>*, <sup>Q</sup>*A*, *<sup>α</sup>*)) **<sup>15</sup>** *<sup>α</sup>* <sup>=</sup> min(*α*, Algorithm <sup>2</sup> (<sup>P</sup> *<sup>A</sup>*, <sup>Q</sup>*B*, *<sup>α</sup>*)) **<sup>16</sup>** *<sup>α</sup>* <sup>=</sup> min(*α*, Algorithm <sup>2</sup> (P*B*, <sup>Q</sup>*A*, *<sup>α</sup>*)) **<sup>17</sup>** *<sup>α</sup>* <sup>=</sup> min(*α*, Algorithm <sup>2</sup> (P*B*, <sup>Q</sup>*B*, *<sup>α</sup>*)) **<sup>18</sup> end <sup>19</sup> return** *α*

**Algorithm 2:** Minimum Distance Between Two Bernstein Polynomials

**Remark 1.** *Note that Algorithm 2 can also be employed to compute the minimum distance between a Bernstein polynomial and a point or a convex shape. This is shown in Figure 19b.*

**Figure 19.** (**a**) Minimum distance between curves. (**b**) Minimum distance between a curve and a polygon. All distances are measured to the blue curve. A red curve or polygon indicates that a collision exists.
