Using surface_intersction.f90 from Python.

In the process, started the module `_surface_intersection` to
break off pieces used during intersection (the process is
incomplete).

docs/algorithm-helpers.rst

@@ -26,10 +26,10 @@ This is to help with the exposition of the computation and
 .. autofunction:: bezier._curve_helpers._newton_refine
 .. autoclass:: bezier._surface_helpers._IntersectionClassification
    :members:
-.. autofunction:: bezier._surface_helpers.newton_refine
 .. autofunction:: bezier._surface_helpers.classify_intersection
 .. autofunction:: bezier._surface_helpers._jacobian_det
 .. autofunction:: bezier._surface_helpers.two_by_two_det
+.. autofunction:: bezier._surface_intersection._newton_refine
 .. autofunction:: bezier._algebraic_intersection.bezier_roots
 .. autofunction:: bezier._algebraic_intersection.lu_companion
 .. autofunction:: bezier._algebraic_intersection.bezier_value_check

docs/native-libraries.rst

@@ -134,6 +134,7 @@ The C headers for ``libbezier`` will be included in the installed package
        curve_intersection.h
        helpers.h
        surface.h
+       surface_intersection.h
      bezier.h
 
 .. testcleanup:: show-headers, show-lib, show-dll, show-pxd, os-x-dylibs
@@ -162,6 +163,7 @@ Cython declaration files are provided:
      _curve_intersection.pxd
      _helpers.pxd
      _surface.pxd
+     _surface_intersection.pxd
 
 For example, ``cimport bezier._curve`` will provide all the functions
 in ``bezier/curve.h``.

setup_helpers.py

@@ -130,6 +130,18 @@
     'curve',
     'curve_intersection',
 )
+FORTRAN_MODULES['surface_intersection'] = (
+    'types',
+    'helpers',
+    os.path.join(QUADPACK_DIR, 'd1mach'),
+    os.path.join(QUADPACK_DIR, 'dqelg'),
+    os.path.join(QUADPACK_DIR, 'dqpsrt'),
+    os.path.join(QUADPACK_DIR, 'dqk21'),
+    os.path.join(QUADPACK_DIR, 'dqagse'),
+    'curve',
+    'surface',
+    'surface_intersection',
+)
 FORTRAN_SOURCE_FILENAME = os.path.join('src', 'bezier', '{}.f90')
 OBJECT_FILENAME = os.path.join('src', 'bezier', '{}.o')
 SPEEDUP_FILENAME = os.path.join('src', 'bezier', '_{}_speedup.c')

src/bezier/__init__.py

@@ -52,6 +52,11 @@
     _HAS_CURVE_INTERSECTION_SPEEDUP = True
 except ImportError:  # pragma: NO COVER
     _HAS_CURVE_INTERSECTION_SPEEDUP = False
+try:
+    import bezier._surface_intersection_speedup  # noqa: F401
+    _HAS_SURFACE_INTERSECTION_SPEEDUP = True
+except ImportError:  # pragma: NO COVER
+    _HAS_SURFACE_INTERSECTION_SPEEDUP = False
 
 
 # NOTE: The ``__version__`` and ``__author__`` are hard-coded here, rather

src/bezier/_surface_helpers.py

@@ -39,6 +39,7 @@
 from bezier import _curve_helpers
 from bezier import _helpers
 from bezier import _intersection_helpers
+from bezier import _surface_intersection
 from bezier import curved_polygon
 try:
     from bezier import _surface_speedup
@@ -966,172 +967,6 @@ def _jacobian_det(nodes, degree, st_vals):
             bs_bt_vals[:, 1] * bs_bt_vals[:, 2])
 
 
-def newton_refine_solve(jac_both, x_val, surf_x, y_val, surf_y):
-    r"""Helper for :func:`newton_refine`.
-
-    We have a system:
-
-    .. code-block:: rest
-
-       [A C][ds] = [E]
-       [B D][dt]   [F]
-
-    This is not a typo, ``A->B->C->D`` matches the data in ``jac_both``.
-    We solve directly rather than using a linear algebra utility:
-
-    .. code-block:: rest
-
-       ds = (D E - C F) / (A D - B C)
-       dt = (A F - B E) / (A D - B C)
-
-    Args:
-        jac_both (numpy.ndarray): A 1x4 matrix of entries in a Jacobian.
-        x_val (float): An ``x``-value we are trying to reach.
-        surf_x (float): The actual ``x``-value we are currently at.
-        y_val (float): An ``y``-value we are trying to reach.
-        surf_y (float): The actual ``x``-value we are currently at.
-
-    Returns:
-        Tuple[float, float]: The pair of values the solve the
-         linear system.
-    """
-    a_val, b_val, c_val, d_val = jac_both[0, :]
-    #       and
-    e_val = x_val - surf_x
-    f_val = y_val - surf_y
-    # Now solve:
-    denom = a_val * d_val - b_val * c_val
-    delta_s = (d_val * e_val - c_val * f_val) / denom
-    delta_t = (a_val * f_val - b_val * e_val) / denom
-    return delta_s, delta_t
-
-
-def newton_refine(nodes, degree, x_val, y_val, s, t):
-    r"""Refine a solution to :math:`B(s, t) = p` using Newton's method.
-
-    Computes updates via
-
-    .. math::
-
-       \left[\begin{array}{c}
-           0 \\ 0 \end{array}\right] \approx
-           \left(B\left(s_{\ast}, t_{\ast}\right) -
-           \left[\begin{array}{c} x \\ y \end{array}\right]\right) +
-           \left[\begin{array}{c c}
-               B_s\left(s_{\ast}, t_{\ast}\right) &
-               B_t\left(s_{\ast}, t_{\ast}\right) \end{array}\right]
-           \left[\begin{array}{c}
-               \Delta s \\ \Delta t \end{array}\right]
-
-    For example, (with weights
-    :math:`\lambda_1 = 1 - s - t, \lambda_2 = s, \lambda_3 = t`)
-    consider the surface
-
-    .. math::
-
-       B(s, t) =
-           \left[\begin{array}{c} 0 \\ 0 \end{array}\right] \lambda_1^2 +
-           \left[\begin{array}{c} 1 \\ 0 \end{array}\right]
-               2 \lambda_1 \lambda_2 +
-           \left[\begin{array}{c} 2 \\ 0 \end{array}\right] \lambda_2^2 +
-           \left[\begin{array}{c} 2 \\ 1 \end{array}\right]
-               2 \lambda_1 \lambda_3 +
-           \left[\begin{array}{c} 2 \\ 2 \end{array}\right]
-               2 \lambda_2 \lambda_1 +
-           \left[\begin{array}{c} 0 \\ 2 \end{array}\right] \lambda_3^2
-
-    and the point
-    :math:`B\left(\frac{1}{4}, \frac{1}{2}\right) =
-    \frac{1}{4} \left[\begin{array}{c} 5 \\ 5 \end{array}\right]`.
-
-    Starting from the **wrong** point
-    :math:`s = \frac{1}{2}, t = \frac{1}{4}`, we have
-
-    .. math::
-
-       \begin{align*}
-       \left[\begin{array}{c} x \\ y \end{array}\right] -
-          B\left(\frac{1}{2}, \frac{1}{4}\right) &= \frac{1}{4}
-          \left[\begin{array}{c} -1 \\ 2 \end{array}\right] \\
-       DB\left(\frac{1}{2}, \frac{1}{4}\right) &= \frac{1}{2}
-           \left[\begin{array}{c c} 3 & 2 \\ 1 & 6 \end{array}\right] \\
-       \Longrightarrow \left[\begin{array}{c}
-           \Delta s \\ \Delta t \end{array}\right] &= \frac{1}{32}
-           \left[\begin{array}{c} -10 \\ 7 \end{array}\right]
-       \end{align*}
-
-    .. image:: images/newton_refine_surface.png
-       :align: center
-
-    .. testsetup:: newton-refine-surface
-
-       import numpy as np
-       import bezier
-       from bezier._surface_helpers import newton_refine
-
-    .. doctest:: newton-refine-surface
-
-       >>> nodes = np.asfortranarray([
-       ...     [0.0, 0.0],
-       ...     [1.0, 0.0],
-       ...     [2.0, 0.0],
-       ...     [2.0, 1.0],
-       ...     [2.0, 2.0],
-       ...     [0.0, 2.0],
-       ... ])
-       >>> surface = bezier.Surface(nodes, degree=2)
-       >>> surface.is_valid
-       True
-       >>> (x_val, y_val), = surface.evaluate_cartesian(0.25, 0.5)
-       >>> x_val, y_val
-       (1.25, 1.25)
-       >>> s, t = 0.5, 0.25
-       >>> new_s, new_t = newton_refine(nodes, 2, x_val, y_val, s, t)
-       >>> 32 * (new_s - s)
-       -10.0
-       >>> 32 * (new_t - t)
-       7.0
-
-    .. testcleanup:: newton-refine-surface
-
-       import make_images
-       make_images.newton_refine_surface(
-           surface, x_val, y_val, s, t, new_s, new_t)
-
-    Args:
-        nodes (numpy.ndarray): Array of nodes in a surface.
-        degree (int): The degree of the surface.
-        x_val (float): The :math:`x`-coordinate of a point
-            on the surface.
-        y_val (float): The :math:`y`-coordinate of a point
-            on the surface.
-        s (float): Approximate :math:`s`-value to be refined.
-        t (float): Approximate :math:`t`-value to be refined.
-
-    Returns:
-        Tuple[float, float]: The refined :math:`s` and :math:`t` values.
-    """
-    lambda1 = 1.0 - s - t
-    (surf_x, surf_y), = evaluate_barycentric(
-        nodes, degree, lambda1, s, t)
-    if surf_x == x_val and surf_y == y_val:
-        # No refinement is needed.
-        return s, t
-
-    # NOTE: This function assumes ``dimension==2`` (i.e. since ``x, y``).
-    jac_nodes = jacobian_both(nodes, degree, 2)
-
-    # The degree of the jacobian is one less.
-    jac_both = evaluate_barycentric(
-        jac_nodes, degree - 1, lambda1, s, t)
-
-    # The first column of the jacobian matrix is B_s (i.e. the
-    # left-most values in ``jac_both``).
-    delta_s, delta_t = newton_refine_solve(
-        jac_both, x_val, surf_x, y_val, surf_y)
-    return s + delta_s, t + delta_t
-
-
 def update_locate_candidates(
         candidate, next_candidates, x_val, y_val, degree):
     """Update list of candidate surfaces during geometric search for a point.
@@ -1231,13 +1066,13 @@ def locate_point(nodes, degree, x_val, y_val):
     # We take the average of all centroids from the candidates
     # that may contain the point.
     s_approx, t_approx = mean_centroid(candidates)
-    s, t = newton_refine(
+    s, t = _surface_intersection.newton_refine(
         nodes, degree, x_val, y_val, s_approx, t_approx)
 
     actual = evaluate_barycentric(nodes, degree, 1.0 - s - t, s, t)
     expected = np.asfortranarray([[x_val, y_val]])
     if not _helpers.vector_close(actual, expected, eps=_LOCATE_EPS):
-        s, t = newton_refine(
+        s, t = _surface_intersection.newton_refine(
             nodes, degree, x_val, y_val, s, t)
     return s, t
 

src/bezier/_surface_intersection.pxd

@@ -0,0 +1,20 @@
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     https://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""Cython wrapper for ``surface_intersection.f90``."""
+
+
+cdef extern from "bezier/surface_intersection.h":
+    void newton_refine_surface(
+        int *num_nodes, double *nodes, int *degree,
+        double *x_val, double *y_val, double *s, double *t,
+        double *updated_s, double *updated_t)