dhermes / bezier

Helper for Bézier Curves, Triangles, and Higher Order Objects
Apache License 2.0
255 stars 35 forks source link

[QUESTION] how to find the maximum y or or minimum y point on a Bezier Curve of degree 3. #235

Closed piscvau closed 4 years ago

piscvau commented 4 years ago

Thiese are questions from maths incompetent user of the package, so please be tolerant!...

For a degree 3 bezier Curve, is there an easy way to determine the point on the curve which has the minimum y coordinate or alternatively the maximum y coordinate.

Also is it possible to split the bezier curve at a chosen point defined by its coordinates. (spun off as #236)

dhermes commented 4 years ago

So a few remarks here:

  1. These questions are likely better as StackOverflow questions than issues. (I'd even get some of that coveted StackOverflow reputation out of answering them.) Though admittedly I'm not checking StackOverflow for bezier questions.
  2. There are 2 questions here; it's best to keep an issue focused on a single question (so you should file two issues next time). I have split off #236 for your second question.

I'll follow this up with another comment with some ideas on how to find the extremal coordinates.

dhermes commented 4 years ago

Consider a curve

B(s) = [x(s)]
       [y(s)]

Finding the extremal values of x(s) for s in [0, 1] corresponds to comparing all of x(0), x(s*), ..., x(1) where s* can be any of the critical points of x(s), i.e. solutions to x'(s) = 0 inside [0, 1].

Consider the following quadratic curve

In [1]: import bezier

In [2]: import numpy as np

In [3]: import seaborn

In [4]: seaborn.set()

In [5]: nodes = np.asfortranarray([
   ...:     [0.0, 3.0, 1.0],
   ...:     [0.0, 1.0, 3.0],
   ...: ])

In [6]: curve = bezier.Curve.from_nodes(nodes)

In [7]: ax = curve.plot(256)

In [8]: ax.axis("scaled")
Out[8]: (-0.09, 1.89, -0.15000000000000002, 3.15)

In [9]: ax.figure.savefig("issue-235.png", dpi=300, bbox_inches="tight")

issue-235

We can convert it to the hodograph curve by taking the first finite difference and multiplying by the degree. In order to keep things simple, we can also focus on x'(s) and eliminate y'(s) from the hodograph by just setting the y-coordinates to the function y(s) = s.

In [10]: hodograph_nodes = curve.degree * (nodes[:, 1:] - nodes[:, :-1])

In [11]: hodograph_nodes[1, :] = np.arange(curve.degree) / (curve.degree - 1)

In [12]: x_prime = bezier.Curve.from_nodes(hodograph_nodes)

To find places where x'(s) = 0 we artificially create a line where x(s) = 0 and the y-values are on the same line as we used above.

In [13]: all_zero = np.asfortranarray([
    ...:     [0.0, 0.0],
    ...:     [0.0, 1.0],
    ...: ])

In [14]: all_zero_curve = bezier.Curve.from_nodes(all_zero)

In [15]: x_prime.intersect(all_zero_curve)
Out[15]:
array([[0.6],
       [0.6]])

This means there is exactly one critical point at s* = 0.6. Now to find the extremal x-values we can plug in the endpoints and the critical points:

In [16]: curve.evaluate_multi(np.asfortranarray([0.0, 0.6, 1.0]))
Out[16]:
array([[0.  , 1.8 , 1.  ],
       [0.  , 1.56, 3.  ]])

There is of course another option which is to evaluate the curve at a lot of points and then find the extrema:

In [17]: s_values = np.linspace(0.0, 1.0, 2**14 + 1)

In [18]: evaluated = curve.evaluate_multi(s_values)

In [19]: np.argmin(evaluated[0, :])
Out[19]: 0

In [20]: np.argmax(evaluated[0, :])
Out[20]: 9830

In [21]: evaluated[:, (0, 9830)]
Out[21]:
array([[0.        , 1.8       ],
       [0.        , 1.55992188]])

In [22]: s_values[(0, 9830),]
Out[22]: array([0.        , 0.59997559])

but clearly this will be less precice and will likely require more computation to find the extrema.


I'm pre-emptively closing this issue @piscvau because I think it sufficiently answers the question. Feel free to request that I re-open, I'm happy to continue the discussion.

dhermes commented 4 years ago

I suppose yet another approach would be to find critical points in exact arithmetic, but this mostly defeats the purpose of this library:

In [23]: import sympy

In [24]: x_sym, _ = curve.to_symbolic()

In [25]: s = sympy.Symbol("s")

In [26]: x_sym.diff(s)
Out[26]: 6 - 10*s

In [27]: sympy.solve([x_sym.diff(s)], [s])
Out[27]: {s: 3/5}
piscvau commented 3 years ago

Thanks a lot for this detailed answer which allowed me to write the function I needed. I am very grateful. Here is the code I ended up with:

import bezier
import numpy as np
from matplotlib import pyplot
import seaborn
seaborn.set()

nodes = np.asfortranarray([
    [0.0, 1.0, 3.0, 4.0, 6.0],
    [0.0, 3.0, 1.0, 1.0, 1.0],
])
curve = bezier.Curve.from_nodes(nodes)
ax = curve.plot(256)
ax.axis("scaled")

hodograph_nodes = curve.degree * (nodes[:, 1:] - nodes[:, :-1])
hodograph_nodes[0, :] = np.arange(curve.degree) / (curve.degree - 1)
y_prime = bezier.Curve.from_nodes(hodograph_nodes)
all_zero = np.asfortranarray([
    hodograph_nodes[0, :],
    np.zeros(hodograph_nodes[0].size)])
all_zero_curve = bezier.Curve.from_nodes(all_zero)
value = y_prime.intersect(all_zero_curve)
point = curve.evaluate(value[0][0])

ax.plot(point[0, :], point[1, :], marker='o', linestyle='None', color='black')

pyplot.show()