Use tangent-tangent cross product when checking Newton convergence

[dhermes]

During curve-curve intersection, before exiting (also, Fortran code) from Newton iteration, the actual cross product of the normalized tangent directions should be checked. From the commit message of 270a9a0dad4fefbcaea499131bde2b60c455e842:

... the normalized tangent vectors (i.e. T(s) = B(s) / ||B(s)||) were near parallel, but not parallel. (For comparison, the actually tangent test cases have cross product T1(s) x T2(t) on the order of 10^{-6} or smaller.)

[dhermes]

In [1]: from tests.functional import utils

In [2]: import numpy as np

In [3]: from bezier import _curve_helpers

In [4]: CURVES, INTERSECTIONS = utils.curve_intersections_info()

In [5]: def normalized_tangents(s, id1, t, id2):
   ...:     curve1 = CURVES[id1].curve
   ...:     tangent1 = _curve_helpers.evaluate_hodograph(s, curve1._nodes)
   ...:     n_tangent1 = tangent1 / np.linalg.norm(tangent1, ord=2)
   ...: 
   ...:     curve2 = CURVES[id2].curve
   ...:     tangent2 = _curve_helpers.evaluate_hodograph(t, curve2._nodes)
   ...:     n_tangent2 = tangent2 / np.linalg.norm(tangent2, ord=2)
   ...: 
   ...:     cross_product = (
   ...:         n_tangent1[0, 0] * n_tangent2[1, 0] -
   ...:         n_tangent1[1, 0] * n_tangent2[0, 0])
   ...:     return cross_product
   ...: 

In [6]: normalized_tangents(0.5, '1', 0.5, '6')
Out[6]: 0.0

In [7]: normalized_tangents(2.0 / 3.0, '14', 1.0 / 3.0, '15')
Out[7]: 0.0

In [8]: normalized_tangents(1.0, '1', 0.0, '18')
Out[8]: 0.0

In [9]: normalized_tangents(0.5, '10', 3.0 / 10.0, '23')
Out[9]: 0.0

In [10]: normalized_tangents(0.5, '28', 0.5, '29')
Out[10]: 0.0

In [11]: normalized_tangents(1.0 / 3.0, '38', 0.5, '39')
Out[11]: 0.0

In [12]: normalized_tangents(0.5, '1', 0.5, '55')
Out[12]: 0.0

In [13]: normalized_tangents(0.5, '56', 0.5, '57')
Out[13]: 0.0

In [14]: normalized_tangents(0.0, '58', 0.0, '59')
Out[14]: 0.0

In [15]: normalized_tangents(0.25, '60', 0.0, '59')
Out[15]: 0.0

In [16]: normalized_tangents(1.0 / 3.0, '61', 1.0 / 3.0, '62')
Out[16]: 0.0

In [17]: normalized_tangents(
    ...:     float.fromhex('0x1.6225dd18931fap-3'), '63',
    ...:     float.fromhex('0x1.1aa1c972f3091p-1'), '64')
    ...:
Out[17]: 0.0

In [18]: normalized_tangents(0.375, '65', 0.75, '66')
Out[18]: 0.0

[dhermes]

In [1]: from tests.functional import utils

In [2]: import numpy as np

In [3]: from bezier import _curve_helpers

In [4]: from bezier import _intersection_helpers

In [5]: CURVES, INTERSECTIONS = utils.curve_intersections_info()

In [6]: import bezier

In [7]: bezier._HAS_SPEEDUP
Out[7]: False

In [8]: original_newton_iterate = _intersection_helpers.newton_iterate

In [9]: def custom_newton_iterate(evaluate_fn, s, t):
   ...:     converged, new_s, new_t = original_newton_iterate(evaluate_fn, s, t)
   ...:     if isinstance(evaluate_fn, _intersection_helpers.NewtonSimpleRoot):
   ...:         tangent1 = _curve_helpers.evaluate_multi(
   ...:             evaluate_fn.first_deriv1, np.array([new_s]))
   ...:         tangent2 = _curve_helpers.evaluate_multi(
   ...:             evaluate_fn.first_deriv2, np.array([new_t]))
   ...:         n_tangent1 = tangent1 / np.linalg.norm(tangent1, ord=2)
   ...:         n_tangent2 = tangent2 / np.linalg.norm(tangent2, ord=2)
   ...:         cross_product = (
   ...:             n_tangent1[0, 0] * n_tangent2[1, 0] -
   ...:             n_tangent1[1, 0] * n_tangent2[0, 0])
   ...:         print('    Converged: {}'.format(converged))
   ...:         print('            s: {!r}'.format(new_s))
   ...:         print('            t: {!r}'.format(new_t))
   ...:         print('Cross-product: {!r}'.format(cross_product))
   ...: 
   ...:     return converged, new_s, new_t
   ...: 

In [10]: _intersection_helpers.newton_iterate = custom_newton_iterate

In [11]: def intersect(id1, id2):
    ...:     curve1 = CURVES[id1].curve
    ...:     curve2 = CURVES[id2].curve
    ...:     try:
    ...:         curve1.intersect(curve2)
    ...:         return None
    ...:     except Exception as exc:
    ...:         return exc
    ...:

In [12]: intersect('1', '6')

In [13]: intersect('14', '15')
    Converged: False
            s: 0.66668701171875
            t: 0.33331298828125
Cross-product: 0.0

In [14]: intersect('1', '18')

In [15]: intersect('10', '23')
    Converged: True
            s: 0.5
            t: 0.3
Cross-product: 0.0
    Converged: True
            s: 0.5
            t: 0.3
Cross-product: 0.0

In [16]: intersect('28', '29')
    Converged: True
            s: 0.5
            t: 0.5
Cross-product: 0.0
    Converged: True
            s: 0.5
            t: 0.5
Cross-product: 0.0

In [17]: intersect('38', '39')
    Converged: False
            s: 0.3333343505867286
            t: 0.5000015258800931
Cross-product: -3.051760185957283e-06
    Converged: False
            s: 0.33333163791257236
            t: 0.4999974568688584
Cross-product: 5.086262283066492e-06

In [18]: intersect('1', '55')
    Converged: True
            s: 0.5
            t: 0.5
Cross-product: 0.0
    Converged: True
            s: 0.5
            t: 0.5
Cross-product: 0.0

In [19]: intersect('56', '57')
    Converged: False
            s: 0.4999814162744777
            t: 0.49996283254887297
Cross-product: -1.3811830307868217e-09
Out[19]: NotImplementedError('Unsupported multiplicity...')

In [20]: intersect('58', '59')
    Converged: False
            s: 0.99993896484375
            t: 0.99993896484375
Cross-product: 0.0
    Converged: False
            s: 0.9999923706181448
            t: 0.9999949137454298
Cross-product: -2.5431596226693855e-06

In [21]: intersect('60', '59')
    Converged: False
            s: 0.25006103515625
            t: 0.99993896484375
Cross-product: 0.0
    Converged: False
            s: 0.25000762939686166
            t: 0.9999949137354255
Cross-product: 2.543164624961945e-06

In [22]: intersect('61', '62')
Out[22]: NotImplementedError('The number of candidate intersections...')

In [23]: intersect('63', '64')
    Converged: False
            s: 0.1729234820690634
            t: 0.552014725012572
Cross-product: 5.152259668683712e-06
    Converged: True
            s: 0.9725947698747117
            t: 0.008262384061543731
Cross-product: -0.8787137527215808
    Converged: True
            s: 0.6666077520083393
            t: 0.8877064436020611
Cross-product: 0.9601657637176324

In [24]: intersect('65', '66')
    Converged: True
            s: 0.375
            t: 0.75
Cross-product: 0.0
    Converged: True
            s: 0.375
            t: 0.75
Cross-product: 0.0
    Converged: True
            s: 0.8653254924944822
            t: 0.33184901501103564
Cross-product: -0.48982055304770056
    Converged: True
            s: 0.5721745075055179
            t: 0.9181509849889643
Cross-product: 0.7527863231309911

In [25]: intersect('73', '74')  # Not tangent
    Converged: True
            s: 0.3150272455455984
            t: 0.010333395340453272
Cross-product: -0.018317404395668535

In [26]: intersect('75', '76')  # Not tangent
    Converged: True
            s: 0.09298759307023234
            t: 0.07544268626657058
Cross-product: -0.026956884210448454

In [27]: intersect('77', '78')  # Not tangent
    Converged: True
            s: 0.3308706505135163
            t: 0.43940853337568364
Cross-product: -0.0011845877680295245

[dhermes]

To be extra thorough, here is what happens for the 61-62 case:

In [28]: from bezier._geometric_intersection import from_linearized

In [29]: from bezier._geometric_intersection import intersect_one_round

In [30]: from bezier._geometric_intersection import Linearization

In [31]: from bezier._geometric_intersection import prune_candidates

In [32]: from bezier._geometric_intersection import SubdividedCurve

In [33]: nodes_first = CURVES['61'].curve._nodes

In [34]: nodes_second = CURVES['62'].curve._nodes

In [35]: curve_first = SubdividedCurve(nodes_first, nodes_first)

In [36]: curve_second = SubdividedCurve(nodes_second, nodes_second)

In [37]: candidate1 = Linearization.from_shape(curve_first)

In [38]: candidate2 = Linearization.from_shape(curve_second)

In [39]: candidates = [(candidate1, candidate2)]

In [40]: intersections = []

In [41]: for _ in range(15):
    ...:     candidates = intersect_one_round(candidates, intersections)
    ...:     if len(candidates) > 64:
    ...:         candidates = prune_candidates(candidates)
    ...:

In [42]: len(candidates)
Out[42]: 74

In [43]: for first, second in candidates:
    ...:     try:
    ...:         from_linearized(first, second, intersections)
    ...:     except:
    ...:         pass
    ...:
    Converged: False
            s: 0.3332737195588042
            t: 0.3332141124889551
Cross-product: -4.201689426006311e-08
    Converged: False
            s: 0.3332776773175139
            t: 0.333222027146069
Cross-product: -3.6625078908607495e-08
    Converged: False
            s: 0.33328176534545406
            t: 0.33323020237511886
Cross-product: -3.144387284370002e-08

...

    Converged: False
            s: 0.3333980011314736
            t: 0.33346267682807973
Cross-product: -4.951612876173098e-08
    Converged: False
            s: 0.3334018088613824
            t: 0.3334702932453512
Cross-product: -5.552090088904466e-08
    Converged: False
            s: 0.3334054819127319
            t: 0.33347764032331795
Cross-product: -6.163874818895603e-08