Using a full Newton's method **always** in from_linearized().

This will **mostly** do the same amount of work: one Newton step. The change is to account for tangent intersections that are "accidentally well-behaved". For example, the curve-curve intersection functional test cases 43 and 44 have linearized segments that intersect inside `[0, 1] x [0, 1]` even though they "shouldn't" because of a tangency. (An alternate approach here would be to check when the convex hulls are tangent and do a check similar to `tangent_bbox_intersection()`.)
dhermes · Mar 5, 2018 · d06430f · d06430f
1 parent c998445
commit d06430f
Show file tree

Hide file tree

Showing 5 changed files with 41 additions and 54 deletions.
diff --git a/docs/curve-curve-intersection.rst b/docs/curve-curve-intersection.rst
@@ -288,7 +288,7 @@ numbers, we can compute the intersection to machine precision:
    ... ]) / 66.0
    >>> max_err = np.max(np.abs(intersections - expected_ints))
    >>> binary_exponent(max_err)
-   -55
+   -54
    >>> s_vals = np.asfortranarray(intersections[0, :])
    >>> points = curve1.evaluate_multi(s_vals)
    >>> expected_pts = np.asfortranarray([
@@ -367,7 +367,7 @@ numbers, we can compute the intersection to machine precision:
    ... ]) / 16.0
    >>> max_err = np.max(np.abs(points - expected_pts))
    >>> binary_exponent(max_err)
-   -51
+   -53
 
 .. image:: images/curves21_and_22.png
    :align: center
@@ -436,7 +436,7 @@ larger.
    ... ]) / 48.0
    >>> max_err = np.max(np.abs(intersections - expected_ints))
    >>> binary_exponent(max_err)
-   -53
+   -52
    >>> s_vals = np.asfortranarray(intersections[0, :])
    >>> points = curve1.evaluate_multi(s_vals)
    >>> expected_pts = np.asfortranarray([
@@ -475,14 +475,14 @@ larger.
    ... ])
    >>> max_err = np.max(np.abs(intersections - expected_ints))
    >>> binary_exponent(max_err)
-   -51
+   -52
    >>> expected_pts = np.asfortranarray([
    ...     [s_val1, s_val2],
    ...     [ 0.375, 0.375 ],
    ... ])
    >>> max_err = np.max(np.abs(points - expected_pts))
    >>> binary_exponent(max_err)
-   -51
+   -52
 
 .. image:: images/curves8_and_27.png
    :align: center
@@ -585,7 +585,7 @@ Detecting Self-Intersections
    ... ]) / 3.0
    >>> max_err = np.max(np.abs(intersections - expected_ints))
    >>> binary_exponent(max_err)
-   -54
+   -53
    >>> s_vals = np.asfortranarray(intersections[0, :])
    >>> left.evaluate_multi(s_vals)
    array([[-0.09375 , -0.25  ],

diff --git a/src/bezier/_geometric_intersection.py b/src/bezier/_geometric_intersection.py
@@ -42,9 +42,8 @@
     _speedup = None
 
 
-# Set the threshold for exponent at half the bits available,
-# this way one round of Newton's method can finish the job
-# by squaring the error.
+# Set the threshold for exponent at half the bits available, this way one round
+# of Newton's method can (usually) finish the job by squaring the error.
 _ERROR_VAL = 0.5**26
 _MAX_INTERSECT_SUBDIVISIONS = 20
 _MAX_CANDIDATES = 64
@@ -758,22 +757,22 @@ def from_linearized(first, second, intersections):
     # pylint: disable=too-many-return-statements
     s, t, success = segment_intersection(
         first.start_node, first.end_node, second.start_node, second.end_node)
-    do_full_newton = False
+    bad_parameters = False
     if success:
         if not (_helpers.in_interval(s, 0.0, 1.0) and
                 _helpers.in_interval(t, 0.0, 1.0)):
-            do_full_newton = True
+            bad_parameters = True
     else:
         if first.error == 0.0 and second.error == 0.0:
             raise ValueError(_UNHANDLED_LINES)
 
-        # Just fall back to a full Newton iteration starting in the middle of
+        # Just fall back to a Newton iteration starting in the middle of
         # the given intervals.
-        do_full_newton = True
+        bad_parameters = True
         s = 0.5
         t = 0.5
 
-    if do_full_newton:
+    if bad_parameters:
         # In the unlikely case that we have parallel segments or segments
         # that intersect outside of [0, 1] x [0, 1], we can still exit
         # if the convex hulls don't intersect.
@@ -783,16 +782,9 @@ def from_linearized(first, second, intersections):
     # Now, promote ``s`` and ``t`` onto the original curves.
     orig_s = (1 - s) * first.curve.start + s * first.curve.end
     orig_t = (1 - t) * second.curve.start + t * second.curve.end
-    if do_full_newton:
-        refined_s, refined_t = _intersection_helpers.full_newton(
-            orig_s, first.curve.original_nodes,
-            orig_t, second.curve.original_nodes)
-    else:
-        # Perform one step of Newton iteration to refine the computed
-        # values of s and t.
-        refined_s, refined_t = _intersection_helpers.newton_refine(
-            orig_s, first.curve.original_nodes,
-            orig_t, second.curve.original_nodes)
+    refined_s, refined_t = _intersection_helpers.full_newton(
+        orig_s, first.curve.original_nodes,
+        orig_t, second.curve.original_nodes)
 
     refined_s, success = _helpers.wiggle_interval(refined_s)
     if not success:
@@ -853,9 +845,10 @@ def add_intersection(s, t, intersections):
     norm_candidate = np.linalg.norm([candidate_s, candidate_t], ord=2)
     for existing_s, existing_t in intersections:
         # NOTE: |(1 - s1) - (1 - s2)| = |s1 - s2| in exact arithmetic, so
-        #       we don't bother comparing to ``candidate_s`` / ``candidate_t``.
-        #       Due to round-off, these may be slightly different, but only
-        #       up to machine precision.
+        #       we just compute ``s1 - s2`` rather than using
+        #       ``candidate_s`` / ``candidate_t``. Due to round-off, these
+        #       differences may be slightly different, but only up to machine
+        #       precision.
         delta_s = s - existing_s
         delta_t = t - existing_t
         norm_update = np.linalg.norm([delta_s, delta_t], ord=2)

diff --git a/src/bezier/curve_intersection.f90 b/src/bezier/curve_intersection.f90
@@ -804,7 +804,7 @@ subroutine from_linearized( &
     integer(c_int), intent(out) :: status
     ! Variables outside of signature.
     real(c_double) :: s, t
-    logical(c_bool) :: success, do_full_newton
+    logical(c_bool) :: success, bad_parameters
 
     status = Status_SUCCESS
     does_intersect = .FALSE.  ! Default value.
@@ -814,25 +814,25 @@ subroutine from_linearized( &
          curve2%nodes(:, 1), curve2%nodes(:, num_nodes2), &
          s, t, success)
 
-    do_full_newton = .FALSE.
+    bad_parameters = .FALSE.
     if (success) then
        if (.NOT. ( &
             in_interval(s, 0.0_dp, 1.0_dp) .AND. &
             in_interval(t, 0.0_dp, 1.0_dp))) then
-          do_full_newton = .TRUE.
+          bad_parameters = .TRUE.
        end if
     else
        ! NOTE: If both curves are lines, the intersection should have
        !       been computed already via ``check_lines()``.
 
        ! Just fall back to a full Newton iteration starting in the middle of
        ! the given intervals.
-       do_full_newton = .TRUE.
+       bad_parameters = .TRUE.
        s = 0.5_dp
        t = 0.5_dp
     end if
 
-    if (do_full_newton) then
+    if (bad_parameters) then
        call convex_hull_collide( &
             num_nodes1, curve1%nodes, POLYGON1, &
             num_nodes2, curve2%nodes, POLYGON2, success)
@@ -848,19 +848,10 @@ subroutine from_linearized( &
     s = (1.0_dp - s) * curve1%start + s * curve1%end_  ! orig_s
     t = (1.0_dp - t) * curve2%start + t * curve2%end_  ! orig_t
 
-    if (do_full_newton) then
-       call full_newton( &
-            s, num_nodes1, root_nodes1, &
-            t, num_nodes2, root_nodes2, &
-            refined_s, refined_t, status)
-    else
-       ! Perform one step of Newton iteration to refine the computed
-       ! values of s and t.
-       call newton_refine_intersect( &
-            s, num_nodes1, root_nodes1, &
-            t, num_nodes2, root_nodes2, &
-            refined_s, refined_t, status)
-    end if
+    call full_newton( &
+         s, num_nodes1, root_nodes1, &
+         t, num_nodes2, root_nodes2, &
+         refined_s, refined_t, status)
 
     if (status /= Status_SUCCESS) then
        return
@@ -1100,9 +1091,10 @@ subroutine add_intersection( &
     ! determined by ``ulps_away``).
     do index_ = 1, num_intersections
        ! NOTE: |(1 - s1) - (1 - s2)| = |s1 - s2| in exact arithmetic, so
-       !       we don't bother comparing to ``(1 - s)`` / ``(1 - t)``.
-       !       Due to round-off, these difference may be slightly different,
-       !       but only up to machine precision.
+       !       we just compute ``s1 - s2`` rather than using
+       !       ``candidate_s`` / ``candidate_t``. Due to round-off, these
+       !       differences may be slightly different, but only up to machine
+       !       precision.
        workspace(1) = s - intersections(1, index_)
        workspace(2) = t - intersections(2, index_)
        if (norm2(workspace) < NEWTON_ERROR_RATIO * norm_candidate) then

diff --git a/tests/fortran/test_surface_intersection.f90 b/tests/fortran/test_surface_intersection.f90
@@ -1876,16 +1876,16 @@ subroutine test_surfaces_intersect(success)
          10, cubic1, 3, 10, cubic2, 3, &
          segment_ends, segments, &
          num_intersected, contained, status)
-    start = 0.60937510406326156_dp
-    end_ = 0.64417397312262581_dp
+    start = 0.60937510406326056_dp
+    end_ = 0.64417397312262270_dp
     case_success = ( &
          num_intersected == 1 .AND. &
          allocated(segment_ends) .AND. &
          all(segment_ends == [3]) .AND. &
          allocated(segments) .AND. &
          segment_check(segments(1), start, end_, 4) .AND. &
-         segment_check(segments(2), 0.97192953004411253_dp, 1.0_dp, 2) .AND. &
-         segment_check(segments(3), 0.0_dp, 0.029255079571202415_dp, 3) .AND. &
+         segment_check(segments(2), 0.97192953004411586_dp, 1.0_dp, 2) .AND. &
+         segment_check(segments(3), 0.0_dp, 0.029255079571205871_dp, 3) .AND. &
          contained == SurfaceContained_NEITHER .AND. &
          status == Status_SUCCESS)
     call print_status(name, case_id, case_success, success)

diff --git a/tests/functional/test_curve_curve.py b/tests/functional/test_curve_curve.py
@@ -73,7 +73,7 @@
         },
         22: {
             (0, 0): 8,  # Established on CentOS 5 (i686 Docker image)
-            (1, 0): 18,  # Established on CentOS 5 (i686 Docker image)
+            (1, 0): 21,  # Established on CentOS 5 (i686 Docker image)
         },
         23: {
             (0, 0): 14,  # Established on Ubuntu 16.04
@@ -157,6 +157,8 @@
         31: {'success': True},
         41: {'success': True},
         42: {'bad_multiplicity': True},
+        43: {'success': True},
+        44: {'success': True},
         45: {'too_many': 74},
     },
     ALGEBRAIC: {},
@@ -171,7 +173,7 @@
     ALGEBRAIC: {},
 }
 INCORRECT_COUNT = {
-    GEOMETRIC: (43, 44),
+    GEOMETRIC: (),
     ALGEBRAIC: (),
 }
 if base_utils.IS_MAC_OS_X or base_utils.IS_PYPY: