[BUG] Geometric Curve-Curve Intersection breaks down when linearized segments are vastly different in size

[KIKI007]

Hi,

I found a bug when using this library.

nodes1 = np.asfortranarray([
    [2.8888563469481592, -1.5080712474161317],
    [10.700182708846356, -8.8105073796645534]])
curve1 = bezier.Curve(nodes1, degree=1)
nodes2 = np.asfortranarray([
    [0, 0.5, 1],
    [0, 0.20000000298023224, 0]])
curve2 = bezier.Curve(nodes2, degree=2)
intersections = curve1.intersect(curve2)

The code return empty list [] but the actual intersection should be as following:

Do I miss to set up some parameters?

Best
Ziqi Wang

[dhermes]

It's definitely giving the wrong answer, I'm not sure if it's a true "bug" or if it's a problem with robustness against round-off and all real inputs. I have computed the exact answer as follows:

Create the bezier.Curve objects

In [1]: import bezier

In [2]: import numpy as np

In [3]: import sympy

In [4]: s, x, y = sympy.symbols("s, x, y")

In [5]: nodes1 = np.asfortranarray(
   ...:     [
   ...:         [2.8888563469481592, -1.5080712474161317],
   ...:         [10.700182708846356, -8.8105073796645534],
   ...:     ]
   ...: )

In [6]: nodes2 = np.asfortranarray([[0.0, 0.5, 1.0], [0.0, 0.20000000298023224, 0.0]])

In [7]: curve1 = bezier.Curve.from_nodes(nodes1)

In [8]: curve2 = bezier.Curve.from_nodes(nodes2)

Convert to Exact Symbolic Representations B(s) and Implicit f(x, y) for each

In [9]: x1_sym, y1_sym = curve1.to_symbolic()

In [10]: x2_sym, y2_sym = curve2.to_symbolic()

In [11]: implicit1 = curve1.implicitize()

In [12]: implicit2 = curve2.implicitize()

Plug in f_2(x_1(s), y_1(s)) = 0 and f_1(x_2(t), y_2(t)) = 0 to solve for the s and t roots

In [13]: root_s, = [root for root in sympy.roots(implicit2.subs({x: x1_sym, y: y1_sym})) if root.is_real and 0 <= root <= 1]

In [14]: root_t, = [root for root in sympy.roots(implicit1.subs({x: x2_sym, y: y2_sym})) if root.is_real and 0 <= root <= 1]

In [15]: root_s
Out[15]: 593656466891835232324098750426407741369/328933483085511928302742135948263082452 - 68719476736*sqrt(2278791103713948174937452054559133115761832488923789846)/82233370771377982075685533987065770613

In [16]: float(root_s)
Out[16]: 0.5433012701892219

In [17]: root_t
Out[17]: -335324269632744659912633/66444492187637950902414 + sqrt(2278791103713948174937452054559133115761832488923789846)/272156640000565046896287744

In [18]: float(root_t)
Out[18]: 0.5

Verify that the only roots on [0, 1] match

In [19]: b1_s = sympy.Matrix([x1_sym, y1_sym]).subs({s: root_s})

In [20]: b2_t = sympy.Matrix([x2_sym, y2_sym]).subs({s: root_t})

In [21]: delta = b1_s - b2_t

In [22]: delta.simplify()

In [23]: delta
Out[23]:
Matrix([
[0],
[0]])

Make Sure This Works in Floating Point Still

In [24]: curve1.evaluate(float(root_s))
Out[24]:
array([[0.5],
       [0.1]])

In [25]: curve2.evaluate(float(root_t))
Out[25]:
array([[0.5],
       [0.1]])

[dhermes]

I've learned a bit more; after 12 refinements the candidate segments from both nodes1 and nodes2 are linearized. There are two pairs of candidates:

curve1: [0, 1]
curve2: [2046/4096, 2047/4096]
----------------------------------------
curve1: [0, 1]
curve2: [2047/4096, 2048/4096]

I think the primary issue is that the size of B_1([0, 1]) is so much bigger than the size of B_2([2046/4096, 2047/4096]) or B_2([2047/4096, 2048/4096]).

This may necessitate a change in the way geometric intersection is done. The current approach (I'm working on documenting it here https://bezier.readthedocs.io/en/latest/algorithms/geometric-curve-intersection.html) is to stop subdividing once one of the curves becomes linear. However, the massive difference in scale (4096 difference in the parameter range) may be causing issues here.

[dhermes]

Based on https://doi.org/10.1016/j.cagd.2019.101791 / https://github.com/dhermes/condition-number-bezier-curve-intersection, I calculated the 2-norm condition number of the intersection to be 6.591740295307331. In other words, the intersection is VERY well-conditioned, and the algorithm is to blame here.

[dhermes]

The issue happens here when trying to intersect B_1([0, 1]) with B_2([2047/4096, 2048/4096]):

In [25]: first, second
Out[25]:
(<bezier.hazmat.geometric_intersection.Linearization at 0x7fce42071880>,
 <bezier.hazmat.geometric_intersection.Linearization at 0x7fce421cc640>)

In [26]: first.start_node
Out[26]: array([ 2.88885635, 10.70018271])

In [27]: first.end_node
Out[27]: array([-1.50807125, -8.81050738])

In [28]: second.start_node
Out[28]: array([0.49975586, 0.09999998])

In [29]: second.end_node
Out[29]: array([0.5, 0.1])

In [30]: s, t, success = segment_intersection(
    ...:     first.start_node, first.end_node, second.start_node, second.end_node
    ...: )

In [31]: s, t, success
Out[31]: (0.5433012701892219, 1.0000000000009721, True)

In [32]: _py_helpers.in_interval(s, 0.0, 1.0)
Out[32]: True

In [33]: _py_helpers.in_interval(t, 0.0, 1.0)
Out[33]: False

[KIKI007]

Thank you very much for investigating my problem.

If it is a robustness problem, how should we handle this case? To give a small random perturbation to the endpoints of the line?

[dhermes]

If it is a robustness problem, how should we handle this case? To give a small random perturbation to the endpoints of the line?

That could work. Another thing that could work would be to manually do the work of bisection by replacing the line with multiple sub-segments of it via Curve.subdivide().

[dhermes]

Just a single subdivision does the trick (though this is still not something you should have to resort to):

In [1]: import bezier                                                                                                                                                                        

In [2]: import numpy as np                                                                                                                                                                   

In [3]: nodes1 = np.asfortranarray( 
   ...:     [ 
   ...:         [2.8888563469481592, -1.5080712474161317], 
   ...:         [10.700182708846356, -8.8105073796645534], 
   ...:     ] 
   ...: )                                                                                                                                                                                    

In [4]: nodes2 = np.asfortranarray([[0.0, 0.5, 1.0], [0.0, 0.20000000298023224, 0.0]])                                                                                                       

In [5]: curve1 = bezier.Curve.from_nodes(nodes1)                                                                                                                                             

In [6]: curve2 = bezier.Curve.from_nodes(nodes2)                                                                                                                                             

In [7]: curve1.intersect(curve2)                                                                                                                                                             
Out[7]: array([], shape=(2, 0), dtype=float64)

In [8]: curve1a, curve1b = curve1.subdivide()                                                                                                                                                

In [9]: curve1a.intersect(curve2)                                                                                                                                                            
Out[9]: array([], shape=(2, 0), dtype=float64)

In [10]: intersections_b = curve1b.intersect(curve2)                                                                                                                                         

In [11]: intersections_b                                                                                                                                                                     
Out[11]: 
array([[0.08660254],
       [0.5       ]])

In [12]: s_val = 0.5 + 0.5 * intersections_b[0, 0]                                                                                                                                           

In [13]: t_val = intersections_b[1, 0]                                                                                                                                                       

In [14]: curve1.evaluate(s_val)                                                                                                                                                              
Out[14]: 
array([[0.5],
       [0.1]])

In [15]: curve2.evaluate(t_val)                                                                                                                                                              
Out[15]: 
array([[0.5],
       [0.1]])

In [16]: s_val                                                                                                                                                                               
Out[16]: 0.5433012701892219

Also note that 0.5433012701892219 is exactly equal to the value we got when using SymPy to compute the root in exact arithmetic, so that's nice at least. (This just reflects the fact that the root condition number is O(1).)