Taking a look at a numba port of earcut

[Huite]

Here's a trial version that uses a numpy structured array with integers instead of pointers.

Some notes:

• ObjectPool is currently a misnomer, since it's not a pool but just a big block of memory.
• Numba does not seem to like functions that return a single structured array element (np.void). The functions here return the integer index instead.
• There's likely still some porting issues since the answers aren't exactly the same. For the big water example, the last 296 triangles or something?
• From earlier experience, numba does not optimize recursion that well. This uses a stack to push and pop.
• Performance seems a lot worse: 70 ms vs 1 ms (!). This is somewhat surprising to me, almost makes me wonder whether it's allocating somewhere unexpected...

from typing import NamedTuple

import numpy as np
import numba as nb


IntDType = np.int64
FloatDType = np.float64


NodeDType = np.dtype(
    [
        ("i", IntDType),
        ("x", FloatDType),
        ("y", FloatDType),
        ("prev", IntDType),
        ("next", IntDType),
        ("z", IntDType),
        ("prev_z", IntDType),
        ("next_z", IntDType),
        ("steiner", bool),
        ("address", IntDType),
    ]
)
Node = np.void

ITEMSIZE = NodeDType.itemsize
NULL = np.void((-1, np.nan, np.nan, -1, -1, -1, -1, -1, True, -1), dtype=NodeDType)
ListPointer = int


class ObjectPool(NamedTuple):
    allocated: np.ndarray
    size: int
    blocksize: int


@nb.njit(inline="always")
def construct(nodes: ObjectPool, index: int, i: int, x: float, y: float):
    node = nodes.allocated[index]
    node["i"] = i
    node["x"] = x
    node["y"] = y
    node["address"] = index
    node["steiner"] = False
    return index + 1


@nb.njit(inline="always")
def isnull(p) -> bool:
    return p["i"] < 0


@nb.njit(inline="always")
def isvalid(p) -> bool:
    return p["i"] > -1


@nb.njit
def earcut(vertices, hole_indices):
    n_vertex = len(vertices)
    has_holes = len(hole_indices) > 0
    if has_holes:
        outer_length = hole_indices[0]
    else:
        outer_length = len(vertices)

    # create_linked_list for the exterior and eliminate_holes will never exceed
    # the initial memory requirements, as every node is featured only once.
    # Hence, the object pool will not be reallocated.
    size = int(1.5 * n_vertex)
    pool = ObjectPool(
        np.full(size * ITEMSIZE, -1, dtype=np.uint8).view(NodeDType),
        size,
        size,
    )
    index = 0
    outer_node_address, index = create_linked_list(
        pool, index, vertices, 0, outer_length, True
    )

    if outer_node_address == -1:
        return np.empty((0, 3), dtype=IntDType)

    outer_node = pool.allocated[outer_node_address]
    if has_holes:
        outer_node_address, index = eliminate_holes(
            pool, index, vertices, hole_indices, outer_node
        )

    xmin = 0.0
    ymin = 0.0
    size = 0.0
    # if the shape is not too simple, we'll use z-order curve hash later; calculate polygon bbox
    # if n_vertex > 80:
    xmin = vertices[:, 0].min()
    xmax = vertices[:, 0].max()
    ymin = vertices[:, 1].min()
    ymax = vertices[:, 1].max()
    size = max(xmax - xmin, ymax - ymin)
    # else:
    #    xmin = 0.0
    #    ymin = 0.0
    #    inv_size = 0.0

    triangles = earcut_linked(
        pool, index, outer_node_address, n_vertex, xmin, ymin, size
    )
    return triangles.copy()


@nb.njit
def create_linked_list(pool, index, vertices, start, end, clockwise):
    if clockwise == (signed_area(vertices, start, end) > 0):
        step = 1
    else:
        start, end = end - 1, start - 1
        step = -1

    last_address = -1
    for i in range(start, end, step):
        x, y = vertices[i]
        last_address, index = insert_node(pool, index, i, x, y, last_address)

    if last_address != -1:
        last = pool.allocated[last_address]
        _next = pool.allocated[last["next"]]
        if equals(last, _next):
            remove_node(pool, last)
            last = _next

    return last["address"], index


# eliminate colinear or duplicate points
@nb.njit(inline="always")
def filter_points(pool, start_address, end_address) -> ListPointer:
    if end_address == -1:
        end_address = start_address

    start = pool.allocated[start_address]
    end = pool.allocated[end_address]

    p = start
    again = True

    while again or (p["address"] != end["address"]):
        again = False

        p_next = pool.allocated[p["next"]]
        p_prev = pool.allocated[p["prev"]]
        if not p["steiner"] and (equals(p, p_next) or area(p_prev, p, p_next) == 0):
            remove_node(pool, p)
            p = end = p_prev
            if p["address"] == p_next["address"]:
                return -1
            again = True

        else:
            p = p_next

    return end["address"]


@nb.njit(inline="always")
def push(array, value, size):
    array[size] = value
    return size + 1


@nb.njit(inline="always")
def pop(array, size):
    return array[size - 1], size - 1


@nb.njit(inline="always")
def set_triangle(triangles, triangle_count, a, b, c):
    triangles[triangle_count, 0] = a["i"]
    triangles[triangle_count, 1] = b["i"]
    triangles[triangle_count, 2] = c["i"]
    return triangle_count + 1


@nb.njit
def earcut_linked(pool, index, ear_address, n_vertex, xmin, ymin, size):
    # main ear slicing loop which triangulates a polygon (given as a linked _list)
    triangles = np.full((n_vertex * 2, 3), -1, dtype=IntDType)
    triangle_count = 0

    stack = np.full(n_vertex, -1, dtype=IntDType)
    stack[0] = ear_address
    stack_size = 1
    pass_stack = np.full(n_vertex, 0, dtype=np.uint8)
    pass_stack[0] = 0
    pass_size = 1
    # stack = [ear_address]
    # pass_stack = [0]

    while stack_size > 0:
        # while len(stack) > 0:
        address, stack_size = pop(stack, stack_size)
        _pass, pass_size = pop(pass_stack, pass_size)
        # address = stack.pop()
        # _pass = pass_stack.pop()
        if address == -1:
            continue

        ear = pool.allocated[address]
        # interlink polygon nodes in z-order
        if (_pass == 0) and (size != 0.0):
            index_curve(pool, ear, xmin, ymin, size)

        # iterate through ears, slicing them one by one
        stop = ear
        while ear["prev"] != ear["next"]:
            _prev = pool.allocated[ear["prev"]]
            _next = pool.allocated[ear["next"]]

            if (
                is_ear_hashed(pool, ear, xmin, ymin, size)
                if (size != 0)
                else is_ear(pool, ear)
            ):
                # cut off the triangle
                triangle_count = set_triangle(
                    triangles, triangle_count, _prev, ear, _next
                )
                remove_node(pool, ear)

                # skipping the next vertice leads to less sliver triangles
                ear = stop = pool.allocated[_next["next"]]
                continue

            ear = _next
            # if we looped through the whole remaining polygon and can't find any more ears
            if ear["address"] == stop["address"]:
                # try filtering points and slicing again
                if _pass == 0:
                    end_address = filter_points(pool, ear["address"], -1)
                    stack_size = push(stack, end_address, stack_size)
                    pass_size = push(pass_stack, 1, pass_size)
                    # stack.append(end_address)
                    # pass_stack.append(1)

                    # if this didn't work, try curing all small self-intersections locally
                elif _pass == 1:
                    ear_address, triangle_count = cure_local_intersections(
                        pool, ear, triangles, triangle_count
                    )
                    stack_size = push(stack, ear_address, stack_size)
                    pass_size = push(pass_stack, 2, pass_size)
                    # stack.append(ear_address)
                    # pass_stack.append(2)

                    # as a last resort, try splitting the remaining polygon into two
                elif _pass == 2:
                    # This function may add entries to the stack
                    pool, index, stack_size, pass_size = split_earcut(
                        pool, index, stack, stack_size, pass_stack, pass_size, ear
                    )

                break

    return triangles[:triangle_count]


@nb.njit(inline="always")
def split_earcut(pool, index, stack, stack_size, pass_stack, pass_size, start):
    do = True
    a = start

    while do or (a["address"] != start["address"]):
        do = False
        _next = pool.allocated[a["next"]]
        b = pool.allocated[_next["next"]]

        while b["address"] != a["prev"]:
            if a["i"] != b["i"] and is_valid_diagonal(pool, a, b):
                # split the polygon in two by the diagonal
                c_address, pool, index = split_polygon(pool, index, a, b)
                c = pool.allocated[c_address]

                # filter colinear points around the cuts
                a_address = filter_points(pool, a["address"], a["next"])
                c_address = filter_points(pool, c["address"], c["next"])

                # stack.append(c_address)
                # pass_stack.append(0)
                # stack.append(a_address)
                # pass_stack.append(0)

                # run earcut on each half
                stack_size = push(stack, c_address, stack_size)
                pass_size = push(pass_stack, 0, pass_size)
                stack_size = push(stack, a_address, stack_size)
                pass_size = push(pass_stack, 0, pass_size)
                return pool, index, stack_size, pass_size
            else:
                b = pool.allocated[b["next"]]
        a = pool.allocated[a["next"]]
    return pool, index, stack_size, pass_size


# check whether a polygon node forms a valid ear with adjacent nodes
@nb.njit(inline="always")
def is_ear(pool, ear) -> bool:
    a = pool.allocated[ear["prev"]]
    b = ear
    c = pool.allocated[ear["next"]]

    if area(a, b, c) >= 0:
        return False  # reflex, can't be an ear

    # now make sure we don't have other points inside the potential ear
    _next = pool.allocated[ear["next"]]
    p = pool.allocated[_next["next"]]
    p_prev = pool.allocated[p["prev"]]
    p_next = pool.allocated[p["next"]]

    while p["address"] != ear["prev"]:
        if (
            point_in_triangle(
                a["x"], a["y"], b["x"], b["y"], c["x"], c["y"], p["x"], p["y"]
            )
            and area(p_prev, p, p_next) >= 0
        ):
            return False
        p = p_next

    return True


@nb.njit(inline="always")
def is_ear_hashed(pool, ear, xmin, ymin, size):
    a = pool.allocated[ear["prev"]]
    b = ear
    c = pool.allocated[ear["next"]]

    if area(a, b, c) >= 0:
        return False  # reflex, can't be an ear

    # triangle bbox; min & max are calculated like this for speed
    min_tx = (
        (a["x"] if a["x"] < c["x"] else c["x"])
        if a["x"] < b["x"]
        else (b["x"] if b["x"] < c["x"] else c["x"])
    )
    min_ty = (
        (a["y"] if a["y"] < c["y"] else c["y"])
        if a["y"] < b["y"]
        else (b["y"] if b["y"] < c["y"] else c["y"])
    )
    max_tx = (
        (a["x"] if a["x"] > c["x"] else c["x"])
        if a["x"] > b["x"]
        else (b["x"] if b["x"] > c["x"] else c["x"])
    )
    max_ty = (
        (a["y"] if a["y"] > c["y"] else c["y"])
        if a["y"] > b["y"]
        else (b["y"] if b["y"] > c["y"] else c["y"])
    )

    # z-order range for the current triangle bbox;
    zmin = z_order(min_tx, min_ty, xmin, ymin, size)
    zmax = z_order(max_tx, max_ty, xmin, ymin, size)

    # first look for points inside the triangle in increasing z-order
    p_address = ear["next_z"]
    while p_address != -1:
        p = pool.allocated[p_address]
        if p["z"] > zmax:
            break

        p_prev = pool.allocated[p["prev"]]
        p_next = pool.allocated[p["next"]]

        if (
            p["address"] != ear["prev"]
            and p["address"] != ear["next"]
            and point_in_triangle(
                a["x"], a["y"], b["x"], b["y"], c["x"], c["y"], p["x"], p["y"]
            )
            and area(p_prev, p, p_next) >= 0
        ):
            return False
        p_address = p["next_z"]

    # then look for points in decreasing z-order
    p_address = ear["prev_z"]

    while p_address != -1:
        p = pool.allocated[p_address]
        if p["z"] < zmin:
            break

        ear_prev = pool.allocated[ear["prev"]]
        ear_next = pool.allocated[ear["next"]]
        p_prev = pool.allocated[p["prev"]]
        p_next = pool.allocated[p["next"]]

        if (
            p["address"] != ear_prev["address"]
            and p["address"] != ear_next["address"]
            and point_in_triangle(
                a["x"], a["y"], b["x"], b["y"], c["x"], c["y"], p["x"], p["y"]
            )
            and area(p_prev, p, p_next) >= 0
        ):
            return False
        p_address = p["prev_z"]

    return True


# go through all polygon nodes and cure small local self-intersections
@nb.jit(inline="always")
def cure_local_intersections(pool, start, triangles, triangle_count):
    do = True
    p = start

    while do or p["address"] != start["address"]:
        do = False

        p_next = pool.allocated[p["next"]]
        a = pool.allocated[p["prev"]]
        b = pool.allocated[p_next["next"]]

        if (
            not equals(a, b)
            and intersects(a, p, p_next, b)
            and locally_inside(pool, a, b)
            and locally_inside(pool, b, a)
        ):
            triangle_count = set_triangle(triangles, triangle_count, a, p, b)

            # remove two nodes involved
            remove_node(pool, p)
            remove_node(pool, p_next)

            p = start = b

        p = p_next

    address = filter_points(pool, p["address"], -1)
    return address, triangle_count


# link every hole into the outer loop, producing a single-ring polygon without holes
@nb.njit
def eliminate_holes(pool, index, vertices, hole_indices, outer_node):
    n_hole = len(hole_indices)
    hole_addresses = np.zeros(n_hole, IntDType)
    hole_x = np.zeros(n_hole, FloatDType)

    for i, start in enumerate(hole_indices):
        start = hole_indices[i]

        if i < (n_hole - 1):
            end = hole_indices[i + 1]
        else:
            end = len(vertices)
        _list_address, index = create_linked_list(
            pool, index, vertices, start, end, False
        )
        _list = pool.allocated[_list_address]

        _next = pool.allocated[_list["next"]]
        if _list["address"] == _next["address"]:
            _list["steiner"] = True

        address, x = get_leftmost(pool, _list)
        hole_addresses[i] = address
        hole_x[i] = x

    # process holes from left to right
    order = np.argsort(hole_x)
    for i in order:
        address = hole_addresses[i]
        hole = pool.allocated[address]
        pool, index = eliminate_hole(pool, index, hole, outer_node)
        _next = pool.allocated[outer_node["next"]]
        outer_node_address = filter_points(
            pool, outer_node["address"], outer_node["next"]
        )
        outer_node = pool.allocated[outer_node_address]

    return outer_node["address"], index


# find a bridge between vertices that connects hole with an outer ring and and link it
@nb.njit
def eliminate_hole(pool, index, hole, outer_node):
    outer_node_address = find_hole_bridge(pool, hole, outer_node)
    if outer_node_address != -1:
        outer_node = pool.allocated[outer_node_address]
        b_address, pool, index = split_polygon(pool, index, outer_node, hole)
        b = pool.allocated[b_address]
        b_next = pool.allocated[b["next"]]
        filter_points(pool, b["address"], b_next["address"])
    return pool, index


# David Eberly's algorithm for finding a bridge between hole and outer polygon
@nb.njit
def find_hole_bridge(pool, hole, outer_node) -> ListPointer:
    do = True
    p = outer_node
    hx = hole["x"]
    hy = hole["y"]
    qx = -np.inf
    m_address = -1

    # find a segment intersected by a ray from the hole's leftmost point to the left;
    # segment's endpoint with lesser x will be potential connection point
    while do or p["address"] != outer_node["address"]:
        do = False
        p_next = pool.allocated[p["next"]]
        if hy <= p["y"] and hy >= p_next["y"] and p_next["y"] - p["y"] != 0:
            x = p["x"] + (hy - p["y"]) * (p_next["x"] - p["x"]) / (p_next["y"] - p["y"])

            if x <= hx and x > qx:
                qx = x

                if x == hx:
                    if hy == p["y"]:
                        return p["address"]
                    if hy == p_next["y"]:
                        return p_next["address"]

                m = p if p["x"] < p_next["x"] else p_next
                m_address = m["address"]

        p = p_next

    if m_address == -1:
        return -1

    if hx == qx:
        return m["prev"]

    # look for points inside the triangle of hole point, segment intersection and endpoint;
    # if there are no points found, we have a valid connection;
    # otherwise choose the point of the minimum angle with the ray as connection point

    stop = m
    mx = m["x"]
    my = m["y"]
    tan_min = np.inf
    tan = None

    p = pool.allocated[p["next"]]

    while p["address"] != stop["address"]:
        hx_or_qx = hx if hy < my else qx
        qx_or_hx = qx if hy < my else hx

        if (
            hx >= p["x"]
            and p["x"] >= mx
            and point_in_triangle(hx_or_qx, hy, mx, my, qx_or_hx, hy, p["x"], p["y"])
        ):
            tan = abs(hy - p["y"]) / (hx - p["x"])  # tangential

            if (
                tan < tan_min or (tan == tan_min and p["x"] > m["x"])
            ) and locally_inside(pool, p, hole):
                m = p
                tan_min = tan

        p = pool.allocated[p["next"]]

    return m["address"]


# interlink polygon nodes in z-order
@nb.njit(inline="always")
def index_curve(pool, start, xmin, ymin, size) -> None:
    do = True
    p = start

    while do or p["address"] != start["address"]:
        do = False

        if p["z"] == -1:
            p["z"] = z_order(p["x"], p["y"], xmin, ymin, size)

        p["prev_z"] = p["prev"]
        p["next_z"] = p["next"]
        p = pool.allocated[p["next"]]

    p_prev_z = pool.allocated[p["prev_z"]]
    p_prev_z["next_z"] = -1
    p["prev_z"] = -1

    sort_linked(pool, p["address"])
    return


@nb.njit
def sort_linked(pool, _list):
    in_size = 1
    while True:
        p_address = _list
        _list = -1
        tail_address = -1
        n_merges = 0

        while p_address != -1:
            n_merges += 1
            q_address = p_address
            p_size = 0
            for _ in range(in_size):
                p_size += 1
                q = pool.allocated[q_address]
                q_address = q["next_z"]
                if q_address == -1:
                    break

            q_size = in_size
            while (p_size > 0) or ((q_size > 0) and (q_address != -1)):
                if p_size == 0:
                    e_address = q_address
                    q = pool.allocated[q_address]
                    q_address = q["next_z"]
                    q_size -= 1
                elif q_size == 0 or q_address == -1:
                    e_address = p_address
                    p = pool.allocated[p_address]
                    p_address = p["next_z"]
                    p_size -= 1
                else:
                    p = pool.allocated[p_address]
                    q = pool.allocated[q_address]
                    if p["z"] <= q["z"]:
                        e_address = p_address
                        p_address = p["next_z"]
                        p_size -= 1
                    else:
                        e_address = q_address
                        q_address = q["next_z"]
                        q_size -= 1

                if tail_address != -1:
                    tail = pool.allocated[tail_address]
                    tail["next_z"] = e_address
                else:
                    _list = e_address  # Update list head when tail is null

                e = pool.allocated[e_address]
                e["prev_z"] = tail_address
                tail_address = e_address

            p_address = q_address

        if tail_address != -1:
            tail = pool.allocated[tail_address]
            tail["next_z"] = -1

        if n_merges <= 1:
            return _list

        in_size *= 2


# z-order of a point given coords and size of the vertices bounding box
@nb.njit(inline="always")
def z_order(x, y, xmin, ymin, size):
    # coords are transformed into non-negative 15-bit integer range
    x = 32767 * int((x - xmin) / size)
    y = 32767 * int((y - ymin) / size)

    x = (x | (x << 8)) & 0x00FF00FF
    x = (x | (x << 4)) & 0x0F0F0F0F
    x = (x | (x << 2)) & 0x33333333
    x = (x | (x << 1)) & 0x55555555

    y = (y | (y << 8)) & 0x00FF00FF
    y = (y | (y << 4)) & 0x0F0F0F0F
    y = (y | (y << 2)) & 0x33333333
    y = (y | (y << 1)) & 0x55555555

    return x | (y << 1)


# find the leftmost node of a polygon ring
@nb.njit(inline="always")
def get_leftmost(pool, start):
    do = True
    p = start
    leftmost = start

    while do or p["address"] != start["address"]:
        do = False
        if p["x"] < leftmost["x"]:
            leftmost = p
        p = pool.allocated[p["next"]]

    return leftmost["address"], leftmost["x"]


# check if a point lies within a convex triangle
@nb.njit(inline="always")
def point_in_triangle(ax, ay, bx, by, cx, cy, px, py):
    return (
        (cx - px) * (ay - py) - (ax - px) * (cy - py) >= 0
        and (ax - px) * (by - py) - (bx - px) * (ay - py) >= 0
        and (bx - px) * (cy - py) - (cx - px) * (by - py) >= 0
    )


# check if a diagonal between two polygon nodes is valid (lies in polygon interior)
@nb.njit(inline="always")
def is_valid_diagonal(pool, a, b):
    anext = pool.allocated[a["next"]]
    aprev = pool.allocated[a["prev"]]
    bprev = pool.allocated[b["prev"]]
    bnext = pool.allocated[b["next"]]
    return (
        # Doesn't intersect other edges
        anext["i"] != b["i"]
        and aprev["i"] != b["i"]
        and not intersects_polygon(pool, a, b)
        and (
            # locally visible
            (
                locally_inside(pool, a, b)
                and locally_inside(pool, b, a)
                and middle_inside(pool, a, b)
                # does not create opposite-facing sectors
                and (area(aprev, a, bprev) != 0.0 or area(a, bprev, b) != 0.0)
            )
            # special zero-length case
            or (
                equals(a, b) and area(aprev, a, anext) > 0 and area(bprev, b, bnext) > 0
            )
        )
    )


# signed area of a triangle
@nb.njit(inline="always")
def area(p, q, r):
    return (q["y"] - p["y"]) * (r["x"] - q["x"]) - (q["x"] - p["x"]) * (r["y"] - q["y"])


# check if two points are equal
@nb.njit(inline="always")
def equals(p1, p2):
    return p1["x"] == p2["x"] and p1["y"] == p2["y"]


# check if two segments intersect
@nb.njit(inline="always")
def intersects(p1, q1, p2, q2) -> bool:
    o1 = sign(area(p1, q1, p2))
    o2 = sign(area(p1, q1, q2))
    o3 = sign(area(p2, q2, p1))
    o4 = sign(area(p2, q2, q1))
    
    if (o1 != o2) and (o3 != o4):
        return True

    if (o1 == 0 and on_segment(p1, p2, q1)):
        return True
    if (o2 == 0 and on_segment(p1, q2, q1)):
        return True
    if (o3 == 0 and on_segment(p2, p1, q2)):
        return True
    if (o4 == 0 and on_segment(p2, q1, q2)):
        return True
    
    return False


@nb.njit(inline="always")
def on_segment(p, q, r) -> bool:
    return (
        q["x"] <= max(p["x"], r["x"])
        and q["x"] >= min(p["x"], r["x"])
        and q["y"] <= max(p["y"], r["y"])
        and q["y"] >= min(p["y"], r["y"])
    )


@nb.njit(inline="always")
def sign(val) -> int:
    return (0.0 < val) - (val < 0.0)


# check if a polygon diagonal intersects any polygon segments
@nb.njit(inline="always")
def intersects_polygon(pool, a, b):
    do = True
    p = a

    while do or p["address"] != a["address"]:
        do = False
        p_next = pool.allocated[p["next"]]
        if (
            p["i"] != a["i"]
            and p_next["i"] != a["i"]
            and p["i"] != b["i"]
            and p_next["i"] != b["i"]
            and intersects(p, p_next, a, b)
        ):
            return True

        p = p_next

    return False


# check if a polygon diagonal is locally inside the polygon
@nb.njit(inline="always")
def locally_inside(pool, a, b):
    aprev = pool.allocated[a["prev"]]
    anext = pool.allocated[a["next"]]
    if area(aprev, a, anext) < 0:
        return area(a, b, anext) >= 0 and area(a, aprev, b) >= 0
    else:
        return area(a, b, aprev) < 0 or area(a, anext, b) < 0


# check if the middle point of a polygon diagonal is inside the polygon
@nb.njit(inline="always")
def middle_inside(pool, a, b):
    do = True
    p = a
    inside = False
    px = 0.5 * (a["x"] + b["x"])
    py = 0.5 * (a["y"] + b["y"])

    while do or p["address"] != a["address"]:
        do = False
        p_next = pool.allocated[p["next"]]
        if ((p["y"] > py) != (p_next["y"] > py)) and (
            px
            < (p_next["x"] - p["x"]) * (py - p["y"]) / (p_next["y"] - p["y"]) + p["x"]
        ):
            inside = not inside

        p = p_next

    return inside


# link two polygon vertices with a bridge; if the vertices belong to the same ring, it splits polygon into two;
# if one belongs to the outer ring and another to a hole, it merges it into a single ring
@nb.njit(inline="always")
def split_polygon(pool, index, a, b):
    # This is the only place where we don't know how many nodes will be created.
    # Make sure there is room by reserving for the two additional nodes.
    if index + 2 >= pool.size:
        array = pool.allocated
        new_size = pool.size + pool.blocksize
        # numba refuses to allocate NodeDType directly...
        allocated = np.full(new_size * ITEMSIZE, -1, dtype=np.uint8).view(NodeDType)
        allocated[: array.size] = array[:]
        pool = ObjectPool(allocated, new_size, pool.blocksize)

    a2 = pool.allocated[index]
    index = construct(pool, index, a["i"], a["x"], a["y"])
    b2 = pool.allocated[index]
    index = construct(pool, index, b["i"], b["x"], b["y"])

    an = pool.allocated[a["next"]]
    bp = pool.allocated[b["prev"]]

    a["next"] = b["address"]
    b["prev"] = a["address"]

    a2["next"] = an["address"]
    an["prev"] = a2["address"]

    b2["next"] = a2["address"]
    a2["prev"] = b2["address"]

    bp["next"] = b2["address"]
    b2["prev"] = bp["address"]

    return b2["address"], pool, index


# create a node and optionally link it with previous one (in a circular doubly linked _list)
@nb.njit(inline="always")
def insert_node(pool, index, i, x, y, last_address):
    p = pool.allocated[index]
    index = construct(pool, index, i, x, y)

    if last_address == -1:
        p["prev"] = p["address"]
        p["next"] = p["address"]
    else:
        last = pool.allocated[last_address]
        p["next"] = last["next"]
        p["prev"] = last["address"]
        _next = pool.allocated[last["next"]]
        _next["prev"] = p["address"]
        last["next"] = p["address"]

    return p["address"], index


@nb.njit(inline="always")
def remove_node(pool, p) -> None:
    _next = pool.allocated[p["next"]]
    _next["prev"] = p["prev"]
    _prev = pool.allocated[p["prev"]]
    _prev["next"] = p["next"]

    if p["prev_z"] != -1:
        p_prev_z = pool.allocated[p["prev_z"]]
        p_prev_z["next_z"] = p["next_z"]
    if p["next_z"] != -1:
        p_next_z = pool.allocated[p["next_z"]]
        p_next_z["prev_z"] = p["prev_z"]
    return


@nb.njit(inline="always")
def signed_area(vertices, start, end):
    area = 0.0
    x0, y0 = vertices[end - 1]
    for i in range(start, end):
        x1, y1 = vertices[i]
        area += (x0 - x1) * (y1 + y0)
        x0, y0 = x1, y1
    return area