def classify_intersection(intersection, edge_nodes1, edge_nodes2):
    r"""Determine which curve is on the "inside of the intersection".

    .. note::

       This is a helper used only by :meth:`.Surface.intersect`.

    This is intended to be a helper for forming a :class:`.CurvedPolygon`
    from the edge intersections of two :class:`.Surface`-s. In order
    to move from one intersection to another (or to the end of an edge),
    the interior edge must be determined at the point of intersection.

    The "typical" case is on the interior of both edges:

    .. image:: ../images/classify_intersection1.png
       :align: center

    .. testsetup:: classify-intersection1, classify-intersection2,
                   classify-intersection3, classify-intersection4,
                   classify-intersection5, classify-intersection6,
                   classify-intersection7, classify-intersection8,
                   classify-intersection9

       import numpy as np
       import bezier
       from bezier import _curve_helpers
       from bezier._intersection_helpers import Intersection
       from bezier._surface_helpers import classify_intersection

       def hodograph(curve, s):
           return _curve_helpers.evaluate_hodograph(
               s, curve._nodes)

       def curvature(curve, s):
           nodes = curve._nodes
           tangent = _curve_helpers.evaluate_hodograph(
               s, nodes)
           return _curve_helpers.get_curvature(
               nodes, tangent, s)

    .. doctest:: classify-intersection1
       :options: +NORMALIZE_WHITESPACE

       >>> nodes1 = np.asfortranarray([
       ...     [1.0, 1.75, 2.0],
       ...     [0.0, 0.25, 1.0],
       ... ])
       >>> curve1 = bezier.Curve(nodes1, degree=2)
       >>> nodes2 = np.asfortranarray([
       ...     [0.0, 1.6875, 2.0],
       ...     [0.0, 0.0625, 0.5],
       ... ])
       >>> curve2 = bezier.Curve(nodes2, degree=2)
       >>> s, t = 0.25, 0.5
       >>> curve1.evaluate(s) == curve2.evaluate(t)
       array([[ True],
              [ True]])
       >>> tangent1 = hodograph(curve1, s)
       >>> tangent1
       array([[1.25],
              [0.75]])
       >>> tangent2 = hodograph(curve2, t)
       >>> tangent2
       array([[2. ],
              [0.5]])
       >>> intersection = Intersection(0, s, 0, t)
       >>> edge_nodes1 = (nodes1, None, None)
       >>> edge_nodes2 = (nodes2, None, None)
       >>> classify_intersection(intersection, edge_nodes1, edge_nodes2)
       <IntersectionClassification.FIRST: 0>

    .. testcleanup:: classify-intersection1

       import make_images
       make_images.classify_intersection1(
           s, curve1, tangent1, curve2, tangent2)

    We determine the interior (i.e. left) one by using the `right-hand rule`_:
    by embedding the tangent vectors in :math:`\mathbf{R}^3`, we
    compute

    .. _right-hand rule: https://en.wikipedia.org/wiki/Right-hand_rule

    .. math::

       \left[\begin{array}{c}
           x_1'(s) \\ y_1'(s) \\ 0 \end{array}\right] \times
       \left[\begin{array}{c}
           x_2'(t) \\ y_2'(t) \\ 0 \end{array}\right] =
       \left[\begin{array}{c}
           0 \\ 0 \\ x_1'(s) y_2'(t) - x_2'(t) y_1'(s) \end{array}\right].

    If the cross product quantity
    :math:`B_1'(s) \times B_2'(t) = x_1'(s) y_2'(t) - x_2'(t) y_1'(s)`
    is positive, then the first curve is "outside" / "to the right", i.e.
    the second curve is interior. If the cross product is negative, the
    first curve is interior.

    When :math:`B_1'(s) \times B_2'(t) = 0`, the tangent
    vectors are parallel, i.e. the intersection is a point of tangency:

    .. image:: ../images/classify_intersection2.png
       :align: center

    .. doctest:: classify-intersection2
       :options: +NORMALIZE_WHITESPACE

       >>> nodes1 = np.asfortranarray([
       ...     [1.0, 1.5, 2.0],
       ...     [0.0, 1.0, 0.0],
       ... ])
       >>> curve1 = bezier.Curve(nodes1, degree=2)
       >>> nodes2 = np.asfortranarray([
       ...     [0.0, 1.5, 3.0],
       ...     [0.0, 1.0, 0.0],
       ... ])
       >>> curve2 = bezier.Curve(nodes2, degree=2)
       >>> s, t = 0.5, 0.5
       >>> curve1.evaluate(s) == curve2.evaluate(t)
       array([[ True],
              [ True]])
       >>> intersection = Intersection(0, s, 0, t)
       >>> edge_nodes1 = (nodes1, None, None)
       >>> edge_nodes2 = (nodes2, None, None)
       >>> classify_intersection(intersection, edge_nodes1, edge_nodes2)
       <IntersectionClassification.TANGENT_SECOND: 4>

    .. testcleanup:: classify-intersection2

       import make_images
       make_images.classify_intersection2(s, curve1, curve2)

    Depending on the direction of the parameterizations, the interior
    curve may change, but we can use the (signed) `curvature`_ of each
    curve at that point to determine which is on the interior:

    .. _curvature: https://en.wikipedia.org/wiki/Curvature

    .. image:: ../images/classify_intersection3.png
       :align: center

    .. doctest:: classify-intersection3
       :options: +NORMALIZE_WHITESPACE

       >>> nodes1 = np.asfortranarray([
       ...     [2.0, 1.5, 1.0],
       ...     [0.0, 1.0, 0.0],
       ... ])
       >>> curve1 = bezier.Curve(nodes1, degree=2)
       >>> nodes2 = np.asfortranarray([
       ...     [3.0, 1.5, 0.0],
       ...     [0.0, 1.0, 0.0],
       ... ])
       >>> curve2 = bezier.Curve(nodes2, degree=2)
       >>> s, t = 0.5, 0.5
       >>> curve1.evaluate(s) == curve2.evaluate(t)
       array([[ True],
              [ True]])
       >>> intersection = Intersection(0, s, 0, t)
       >>> edge_nodes1 = (nodes1, None, None)
       >>> edge_nodes2 = (nodes2, None, None)
       >>> classify_intersection(intersection, edge_nodes1, edge_nodes2)
       <IntersectionClassification.TANGENT_FIRST: 3>

    .. testcleanup:: classify-intersection3

       import make_images
       make_images.classify_intersection3(s, curve1, curve2)

    When the curves are moving in opposite directions at a point
    of tangency, there is no side to choose. Either the point of tangency
    is not part of any :class:`.CurvedPolygon` intersection

    .. image:: ../images/classify_intersection4.png
       :align: center

    .. doctest:: classify-intersection4
       :options: +NORMALIZE_WHITESPACE

       >>> nodes1 = np.asfortranarray([
       ...     [2.0, 1.5, 1.0],
       ...     [0.0, 1.0, 0.0],
       ... ])
       >>> curve1 = bezier.Curve(nodes1, degree=2)
       >>> nodes2 = np.asfortranarray([
       ...     [0.0, 1.5, 3.0],
       ...     [0.0, 1.0, 0.0],
       ... ])
       >>> curve2 = bezier.Curve(nodes2, degree=2)
       >>> s, t = 0.5, 0.5
       >>> curve1.evaluate(s) == curve2.evaluate(t)
       array([[ True],
              [ True]])
       >>> intersection = Intersection(0, s, 0, t)
       >>> edge_nodes1 = (nodes1, None, None)
       >>> edge_nodes2 = (nodes2, None, None)
       >>> classify_intersection(intersection, edge_nodes1, edge_nodes2)
       <IntersectionClassification.OPPOSED: 2>

    .. testcleanup:: classify-intersection4

       import make_images
       make_images.classify_intersection4(s, curve1, curve2)

    or the point of tangency is a "degenerate" part of two
    :class:`.CurvedPolygon` intersections. It is "degenerate"
    because from one direction, the point should be classified as
    :attr:`~.IntersectionClassification.FIRST` and from another as
    :attr:`~.IntersectionClassification.SECOND`.

    .. image:: ../images/classify_intersection5.png
       :align: center

    .. doctest:: classify-intersection5
       :options: +NORMALIZE_WHITESPACE

       >>> nodes1 = np.asfortranarray([
       ...     [1.0, 1.5, 2.0],
       ...     [0.0, 1.0, 0.0],
       ... ])
       >>> curve1 = bezier.Curve(nodes1, degree=2)
       >>> nodes2 = np.asfortranarray([
       ...     [3.0, 1.5, 0.0],
       ...     [0.0, 1.0, 0.0],
       ... ])
       >>> curve2 = bezier.Curve(nodes2, degree=2)
       >>> s, t = 0.5, 0.5
       >>> curve1.evaluate(s) == curve2.evaluate(t)
       array([[ True],
              [ True]])
       >>> intersection = Intersection(0, s, 0, t)
       >>> edge_nodes1 = (nodes1, None, None)
       >>> edge_nodes2 = (nodes2, None, None)
       >>> classify_intersection(intersection, edge_nodes1, edge_nodes2)
       <IntersectionClassification.TANGENT_BOTH: 6>

    .. testcleanup:: classify-intersection5

       import make_images
       make_images.classify_intersection5(s, curve1, curve2)

    The :attr:`~.IntersectionClassification.TANGENT_BOTH` classification
    can also occur if the curves are "kissing" but share a zero width
    interior at the point of tangency:

    .. image:: ../images/classify_intersection9.png
       :align: center

    .. doctest:: classify-intersection9
       :options: +NORMALIZE_WHITESPACE

       >>> nodes1 = np.asfortranarray([
       ...     [0.0, 20.0, 40.0],
       ...     [0.0, 40.0,  0.0],
       ... ])
       >>> curve1 = bezier.Curve(nodes1, degree=2)
       >>> nodes2 = np.asfortranarray([
       ...     [40.0, 20.0,  0.0],
       ...     [40.0,  0.0, 40.0],
       ... ])
       >>> curve2 = bezier.Curve(nodes2, degree=2)
       >>> s, t = 0.5, 0.5
       >>> curve1.evaluate(s) == curve2.evaluate(t)
       array([[ True],
              [ True]])
       >>> intersection = Intersection(0, s, 0, t)
       >>> edge_nodes1 = (nodes1, None, None)
       >>> edge_nodes2 = (nodes2, None, None)
       >>> classify_intersection(intersection, edge_nodes1, edge_nodes2)
       <IntersectionClassification.TANGENT_BOTH: 6>

    .. testcleanup:: classify-intersection9

       import make_images
       make_images.classify_intersection9(s, curve1, curve2)

    However, if the `curvature`_ of each curve is identical, we
    don't try to distinguish further:

    .. image:: ../images/classify_intersection6.png
       :align: center

    .. doctest:: classify-intersection6
       :options: +NORMALIZE_WHITESPACE

       >>> nodes1 = np.asfortranarray([
       ...     [-0.125 , -0.125 , 0.375 ],
       ...     [ 0.0625, -0.0625, 0.0625],
       ... ])
       >>> curve1 = bezier.Curve(nodes1, degree=2)
       >>> nodes2 = np.asfortranarray([
       ...     [-0.25, -0.25, 0.75],
       ...     [ 0.25, -0.25, 0.25],
       ... ])
       >>> curve2 = bezier.Curve(nodes2, degree=2)
       >>> s, t = 0.5, 0.5
       >>> curve1.evaluate(s) == curve2.evaluate(t)
       array([[ True],
              [ True]])
       >>> hodograph(curve1, s)
       array([[0.5],
              [0. ]])
       >>> hodograph(curve2, t)
       array([[1.],
              [0.]])
       >>> curvature(curve1, s)
       2.0
       >>> curvature(curve2, t)
       2.0
       >>> intersection = Intersection(0, s, 0, t)
       >>> edge_nodes1 = (nodes1, None, None)
       >>> edge_nodes2 = (nodes2, None, None)
       >>> classify_intersection(intersection, edge_nodes1, edge_nodes2)
       Traceback (most recent call last):
         ...
       NotImplementedError: Tangent curves have same curvature.

    .. testcleanup:: classify-intersection6

       import make_images
       make_images.classify_intersection6(s, curve1, curve2)

    In addition to points of tangency, intersections that happen at
    the end of an edge need special handling:

    .. image:: ../images/classify_intersection7.png
       :align: center

    .. doctest:: classify-intersection7
       :options: +NORMALIZE_WHITESPACE

       >>> nodes1a = np.asfortranarray([
       ...     [0.0, 4.5, 9.0 ],
       ...     [0.0, 0.0, 2.25],
       ... ])
       >>> curve1a = bezier.Curve(nodes1a, degree=2)
       >>> nodes2 = np.asfortranarray([
       ...     [11.25, 9.0, 2.75],
       ...     [ 0.0 , 4.5, 1.0 ],
       ... ])
       >>> curve2 = bezier.Curve(nodes2, degree=2)
       >>> s, t = 1.0, 0.375
       >>> curve1a.evaluate(s) == curve2.evaluate(t)
       array([[ True],
              [ True]])
       >>> intersection = Intersection(0, s, 0, t)
       >>> edge_nodes1 = (nodes1a, None, None)
       >>> edge_nodes2 = (nodes2, None, None)
       >>> classify_intersection(intersection, edge_nodes1, edge_nodes2)
       Traceback (most recent call last):
         ...
       ValueError: ('Intersection occurs at the end of an edge',
                    's', 1.0, 't', 0.375)
       >>>
       >>> nodes1b = np.asfortranarray([
       ...     [9.0, 4.5, 0.0],
       ...     [2.25, 2.375, 2.5],
       ... ])
       >>> curve1b = bezier.Curve(nodes1b, degree=2)
       >>> curve1b.evaluate(0.0) == curve2.evaluate(t)
       array([[ True],
              [ True]])
       >>> intersection = Intersection(1, 0.0, 0, t)
       >>> edge_nodes1 = (nodes1a, nodes1b, None)
       >>> classify_intersection(intersection, edge_nodes1, edge_nodes2)
       <IntersectionClassification.FIRST: 0>

    .. testcleanup:: classify-intersection7

       import make_images
       make_images.classify_intersection7(s, curve1a, curve1b, curve2)

    As above, some intersections at the end of an edge are part of
    an actual intersection. However, some surfaces may just "kiss" at a
    corner intersection:

    .. image:: ../images/classify_intersection8.png
       :align: center

    .. doctest:: classify-intersection8
       :options: +NORMALIZE_WHITESPACE

       >>> nodes1 = np.asfortranarray([
       ...     [0.25, 0.0, 0.0, 0.625, 0.5  , 1.0 ],
       ...     [1.0 , 0.5, 0.0, 0.875, 0.375, 0.75],
       ... ])
       >>> surface1 = bezier.Surface(nodes1, degree=2)
       >>> nodes2 = np.asfortranarray([
       ...     [0.0625, -0.25, -1.0, -0.5  , -1.0, -1.0],
       ...     [0.5   ,  1.0 ,  1.0,  0.125,  0.5,  0.0],
       ... ])
       >>> surface2 = bezier.Surface(nodes2, degree=2)
       >>> curve1, _, _ = surface1.edges
       >>> edge_nodes1 = [curve.nodes for curve in surface1.edges]
       >>> curve2, _, _ = surface2.edges
       >>> edge_nodes2 = [curve.nodes for curve in surface2.edges]
       >>> s, t = 0.5, 0.0
       >>> curve1.evaluate(s) == curve2.evaluate(t)
       array([[ True],
              [ True]])
       >>> intersection = Intersection(0, s, 0, t)
       >>> classify_intersection(intersection, edge_nodes1, edge_nodes2)
       <IntersectionClassification.IGNORED_CORNER: 5>

    .. testcleanup:: classify-intersection8

       import make_images
       make_images.classify_intersection8(
           s, curve1, surface1, curve2, surface2)

    .. note::

       This assumes the intersection occurs in :math:`\mathbf{R}^2`
       but doesn't check this.

    .. note::

       This function doesn't allow wiggle room / round-off when checking
       endpoints, nor when checking if the cross product is near zero,
       nor when curvatures are compared. However, the most "correct"
       version of this function likely should allow for some round off.

    Args:
        intersection (.Intersection): An intersection object.
        edge_nodes1 (Tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray]): The
            nodes of the three edges of the first surface being intersected.
        edge_nodes2 (Tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray]): The
            nodes of the three edges of the second surface being intersected.

    Returns:
        IntersectionClassification: The "inside" curve type, based on
        the classification enum.

    Raises:
        ValueError: If the intersection occurs at the end of either
            curve involved. This is because we want to classify which
            curve to **move forward** on, and we can't move past the
            end of a segment.
    """
    if intersection.s == 1.0 or intersection.t == 1.0:
        raise ValueError(
            "Intersection occurs at the end of an edge",
            "s",
            intersection.s,
            "t",
            intersection.t,
        )

    nodes1 = edge_nodes1[intersection.index_first]
    tangent1 = _curve_helpers.evaluate_hodograph(intersection.s, nodes1)
    nodes2 = edge_nodes2[intersection.index_second]
    tangent2 = _curve_helpers.evaluate_hodograph(intersection.t, nodes2)
    if ignored_corner(
        intersection, tangent1, tangent2, edge_nodes1, edge_nodes2
    ):
        return CLASSIFICATION_T.IGNORED_CORNER

    # Take the cross product of tangent vectors to determine which one
    # is more "inside" / "to the left".
    cross_prod = _helpers.cross_product(
        tangent1.ravel(order="F"), tangent2.ravel(order="F")
    )
    if cross_prod < -ALMOST_TANGENT:
        return CLASSIFICATION_T.FIRST

    elif cross_prod > ALMOST_TANGENT:
        return CLASSIFICATION_T.SECOND

    else:
        # NOTE: A more robust approach would take ||tangent1|| and ||tangent2||
        #       into account when comparing (tangent1 x tangent2) to the
        #       "almost zero" threshold. We (for now) avoid doing this because
        #       normalizing the tangent vectors has a "cost" of ~6 flops each
        #       and that cost would happen for **every** single intersection.
        return classify_tangent_intersection(
            intersection, nodes1, tangent1, nodes2, tangent2
        )