Speedup triangulation of edges (napari#7268)

# References and relevant issues

# Description

Similar to napari#7256 this PR is a speedup of face triangulation by avoiding
a lot of reallocation of arrays.

After long investigation, I found that using best way to speed this up
is by avoiding the allocation of additional arrays by using numba.

For shapes that do not have bevels:

On main:


![obraz](https://github.com/user-attachments/assets/b6775d94-a88b-4a99-959f-f4b6f7f768aa)

In this PR:


![obraz](https://github.com/user-attachments/assets/9a3a302c-2eb1-4ecd-8c38-a07617512ede)

PR with cache:


![obraz](https://github.com/user-attachments/assets/6d7f00a8-4d55-4bb2-9f49-955bc746e01f)



As you may see, the time of `triangulate_edge` is reduced from 44s to
17s. Half of it is numba compilation.

With sharp edges that require bevel join, the speedup is even bigger:

Main:


![obraz](https://github.com/user-attachments/assets/e2a55f13-6ba7-4f1d-9ff9-56fa46fdf49c)

PR:

![obraz](https://github.com/user-attachments/assets/6279350e-1e50-441e-a538-4f01b5420d40)


PR with cache:


![obraz](https://github.com/user-attachments/assets/3c992928-b63c-42d7-be6d-d0ec7ec620f5)


This PR also fixes the bevel join 

On main:


![obraz](https://github.com/user-attachments/assets/9a24bc31-9777-4d60-81d5-ff6cc1b0d1a2)


With PR:


![obraz](https://github.com/user-attachments/assets/b167a4d5-9bce-442c-80c6-d5a0c90e5f8d)


In this PR I have added `examples/dev/triangle_edge.py` with a set of
helper functions to visualize triangulation. These are helpful for
debugging of triangulation problems.

# Algorithm explanation

Here I add some explanation to the code implementation of this PR. 

## Calculation of miter vector <a href="#user-content-miter"
id="miter">#</a>

This is a screenshot of miter join:


![obraz](https://github.com/user-attachments/assets/08d2dd78-bb02-4da8-baa5-7b0596e1c3c7)

And let's assume markings as on image:


![obraz](https://github.com/user-attachments/assets/d0cec9ff-3786-4edc-805a-d45e58f33aab)


The $\vec{BA}$ (blue) is half of the normal vector from vertex $n$ to
$n+1$. The $\vec{AC}$ is minus half of the normal vector from vertex
$n-1$ to $n$. The vector $\vec{AD}$ (green) is orthogonal to line $BA$
of length $0.5$

Based on these vectors, we could easily calculate $\vec{BC} = \vec{BA} +
\vec{AC}$. However, the length of $|GC|$ may be different from 1. We
would like to get the vector $\vec{BE}$

We will use here trigonometry and [Intercept
theorem](https://en.wikipedia.org/wiki/Intercept_theorem)

We know that $\frac{|GC|}{|AC|} = \sin{\alpha}$ and $|AC| = \frac{1}{2}$
so $|GC| = \frac{\sin{\alpha}}{2}$


From intercept theorem, we know that $\frac{|BC|}{|GC|} =
\frac{|BE|}{|FE|}$. The length of
$|FE| = \frac{1}{2}$. So $|BE| = |BC| \times \frac{1}{\sin{\alpha}}$.
That mens that $\vec{BE} = \vec{BC} \times \frac{1}{\sin{\alpha}}$

## Bevel limit <a href="#user-content-bevel-limit"
id="bevel-limit">#</a>

In general, we use miter join, however, if the angle is sharp, it may
result in strange artifacts (see napari#6706). So if the join will be too long
(The length of $\vec{BA}$ vector from above), we would like to change to
a bevel join.

So we need to calculate the length of $\vec{BE}$. We would like to avoid
calculation of square root in every turn of the loop.

From Pythagoras theorem we know that $|BE|^2 = |EF|^2 + |BF|^2 =
\frac{1}{2}^2 + |BF|^2$

We also know that $|GA| = \frac{\cos{\alpha}}{2}$ so $|BG| = \frac{1 -
\cos{\alpha}}{2}$. From miter vector calculation and intercept theorem
we know that $|BF| = \frac{|BG|}{\sin{\alpha}} = \frac{1 -
\cos{\alpha}}{2\sin{\alpha}}$

So we know the length of the vector 

$$|BE|^2 = \frac{1}{2}^2 + \left(\frac{1 -
\cos{\alpha}}{2\sin{\alpha}}\right)^2 = \frac{\sin^2{\alpha} + 1^2 -
2\cos{\alpha} + \cos^2{\alpha}}{4\sin^2{\alpha}}
= \frac{2 - 2 \cos{\alpha}}{4\sin^2{\alpha}}
$$

Using trigonometric identity

$$|BE|^2 = \frac{1}{2} \cdot \frac{1 - \cos{\alpha}}{1 - \cos^2{\alpha}}
=
\frac{1}{2} \cdot \frac{1 - \cos{\alpha}}{(1 - \cos{\alpha})(1 +
\cos{\alpha})} =
\frac{1}{2} \cdot \frac{1}{1+\cos{\alpha}}$$ 

So if we would like to check condition $|BE| > limit$ we could use the
following equation:

$$ \frac{1}{2} \cdot \frac{1}{1+\cos{\alpha}} > {limit}^2$$ 

So:

$$ \frac{1}{2 \cdot {limit}^2} - 1 > \cos{\alpha} $$ 

It means that if $\cos\alpha$ is below a given threshold, we should use
bevel join.

## Too long vector in bevel join <a href="#user-content-bevel-cut"
id="bevel-cut">#</a>

As I showed above, this PR is fixing the bevel join by fixing the
position of inner edge of the bevel join. You may see an example in this
image (see the longest yellow vector):


![obraz](https://github.com/user-attachments/assets/58e13821-a1b6-4dd0-a49b-ade97edf9e73)

The schematic of generating long edge:


![obraz](https://github.com/user-attachments/assets/329b2b61-3c9d-4050-9623-7053c1f46248)

As users may use shapes with different length edges, we cannot use a
hard-coded limit for the length of the offset vector that creates the
problems pointed on above image.

However, when calculating the normal vectors, we need to calculate the
edge width, so we may save them for later usage and limit the length of
such problematic offset vector to length of shortest adjusted edge.

In the Bevel limit section, we have an explicit formula for miter offset
length. However, it contains a square root. And calculation of a square
root is expensive on current CPU, so we would like to avoid it.

However, the described problem happens only when $\measuredangle BAC$ is
close to $180\degree$, but in such a situation, the $|\vec{BA} +
\vec{AC}|$ is close to 1 (as resulted vector is close to diameter).

So we could approximate the length of vector $\vec{BC}$ as 1. Then the
estimated length of vector $\vec{BE}$ is $\frac{1}{\sin{\alpha}}$

So if $\frac{1}{\sin{\alpha}}$ is bigger than the length of half of any
adjusted edges, we could calculate a new scaling factor with value

$$\pm 0.5 \cdot\min(vec1\textunderscore{}len,
vec2\textunderscore{}len))$$

# Additional notes

## Artifacts

There are still artifacts in rendering 


![obraz](https://github.com/user-attachments/assets/416df018-f293-425c-bd38-916ba8137934)

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Carol Willing <[email protected]>
Co-authored-by: Juan Nunez-Iglesias <[email protected]>

examples/dev/triangle_edge.py

@@ -0,0 +1,388 @@
+"""Visualize the triangulation algorithms used by the Shapes layer.
+
+This example uses napari layers to draw each of the components of a face and
+edge triangulation in a Shapes layer.
+
+Shapes layers don't "just" draw polygons, ellipses, and so on: each shape,
+as well as its borders, is broken down into a collection of triangles (this
+is called a *triangulation*), which are then sent to OpenGL for drawing:
+drawing triangles is one of the "visualization primitives" in OpenGL and most
+2D and 3D drawing frameworks.
+
+It turns out that converting arbitrary shapes into a collection of triangles
+can be quite tricky: very shallow angles cause errors in the algorithms, and
+can also make certain desirable properties (such as edges not overlapping with
+each other when a polygon makes a sharp turn) actually impossible to achieve
+with fast (single-pass) algorithms.
+
+This example draws the Shapes layer normally, but also overlays all the
+elements of the triangulation: the triangles themselves, and the normal vectors
+on each polygon vertex, from which the triangulation is computed.
+"""
+
+from dataclasses import dataclass, fields
+from functools import partial
+
+try:
+    import numba
+except ImportError:
+    # Create a dummy numba.njit to allow running this script without numba
+    import toolz as tz
+    class numba:
+        @tz.curry
+        def njit(func, cache=True, inline=None):
+            return func
+
+import numpy as np
+
+import napari
+from napari.layers import Shapes
+from napari.layers.shapes._shapes_utils import generate_2D_edge_meshes
+
+
+def generate_regular_polygon(n, radius=1):
+    """Generate vertices of a regular n-sided polygon centered at the origin.
+
+    Parameters
+    ----------
+    n : int
+        The number of sides (vertices).
+    radius : float, optional
+        The radius of the circumscribed circle.
+
+    Returns
+    -------
+    np.ndarray
+        An array of shape (n, 2) containing the vertex coordinates.
+    """
+    angles = np.linspace(0, 2 * np.pi, n, endpoint=False)
+    return np.column_stack((radius * np.cos(angles), radius * np.sin(angles)))
+
+
+def get_reference_edge_triangulation_points(shape: Shapes) -> np.ndarray:
+    """Get the non-accelerated points"""
+    shapes = shape._data_view.shapes
+    path_list = [(x.data, x._closed) for x in shapes]
+    mesh_list = [generate_2D_edge_meshes(path, closed)
+                 for path, closed in path_list]
+    mesh = tuple(np.concatenate(el, axis=0) for el in zip(*mesh_list))
+
+    return mesh[0] + mesh[1]
+
+
+def generate_order_vectors(path_, closed) -> np.ndarray:
+    """Generate the vectors tangent to the path.
+
+    Parameters
+    ----------
+    path_ : np.ndarray, shape (N, 2)
+        A list of 2D path coordinates.
+    closed : bool
+        Whether the coordinates represent a closed polygon or an open
+        line/path.
+
+    Returns
+    -------
+    vec : np.ndarray, shape (N, 2, 2)
+        A set of vectors, defined by a 2D position and a 2D projection.
+    """
+    raw_vecs = np.diff(path_, axis=0)
+    norm = np.linalg.norm(raw_vecs, axis=1, keepdims=True)
+    directions = raw_vecs / norm
+    vec = np.empty((path_.shape[0], 2, 2))
+    vec[:, 0] = path_
+    vec[:-1, 1] = directions
+    if closed:
+        # point from the last vertex towards the first vertex
+        vec[-1, 1] = (
+                (path_[0] - path_[-1]) / np.linalg.norm(path_[-1] - path_[0])
+                )
+    else:
+        # point back towards the penultimate vertex
+        vec[-1, 1] = -vec[-2, 1]
+    return vec
+
+
+def generate_miter_helper_vectors(
+        direction_vectors_array: np.ndarray
+        ) -> np.ndarray:
+    """Generate the miter helper vectors.
+
+    The miter helper vectors are a pair of vectors for each point in the path,
+    which together help define whether a bevel join will be needed. The first
+    vector is half of the normalized direction vector for that vertex, and the
+    second is *minus* half of the normalized direction vector for the previous
+    vertex. Their angle is the angle of the edges at that vertex.
+
+    Parameters
+    ----------
+    direction_vectors_array : array of shape (n, 2)
+        array of normalized (length 1) direction vectors for each point in the
+        path.
+
+    Returns
+    -------
+    array of shape (n, 2, 2)
+        array of miter helper vectors
+    """
+    half_direction = direction_vectors_array.copy()
+    half_direction[:, 1] *= 0.5
+    half_prev_direction = half_direction.copy()
+    half_prev_direction[:, 1] *= -1
+    half_prev_direction[:, 1] = np.roll(half_prev_direction[:, 1], 1, axis=0)
+    half_prev_direction[:, 0] += half_direction[:, 1]
+    return np.concatenate([half_direction, half_prev_direction], axis=0)
+
+
+@numba.njit
+def generate_orthogonal_vectors(direction_vectors: np.ndarray) -> np.ndarray:
+    """Generate the orthogonal vectors to the direction vectors.
+
+    The orthogonal vector starts at the middle of the direction vector and is
+    orthogonal to it in the positive orientation. Its length is half of the
+    direction vector, to align with the normalized edge width.
+
+    Parameters
+    ----------
+    direction_vectors : array, shape (n, 2, 2)
+        The direction vector start points (``direction_vectors[:, 0, :]``) and
+        directions (``direction_vectors[:, 1, :]``).
+
+    Returns
+    -------
+    orthogonal_vectors : array, shape(n, 2, 2)
+        The orthogonal vector start points and directions.
+    """
+    position = 0
+    vector = 1
+    y, x = 0, 1
+    half_direction = 0.5 * direction_vectors[:, 1, :]
+    orthogonal_vectors = direction_vectors.copy()
+    orthogonal_vectors[:, position] = (
+            direction_vectors[:, position] + half_direction
+            )
+
+    orthogonal_vectors[:, vector, y] = -half_direction[:, x]
+    orthogonal_vectors[:, vector, x] = half_direction[:, y]
+    return orthogonal_vectors
+
+
+@numba.njit
+def generate_miter_vectors(
+        mesh: tuple[np.ndarray, np.ndarray, np.ndarray]
+        ) -> np.ndarray:
+    """For each point on path, generate the vectors pointing to miter points.
+
+    Parameters
+    ----------
+    mesh : tuple[np.ndarray]
+        vertices, offsets, and triangles of the mesh corresponding to the edges
+        of a single shape.
+
+    Returns
+    -------
+    np.ndarray, shape (n, 2, 2)
+        Positions and projections of vectors corresponding to the miter points
+        offset from the path points.
+    """
+    vec_points = np.empty((mesh[0].shape[0], 2, 2))
+    vec_points[:, 0] = mesh[0]
+    vec_points[:, 1] = mesh[1]
+    return vec_points
+
+
+@numba.njit
+def generate_edge_triangle_borders(centers, offsets, triangles) -> np.ndarray:
+    """Generate 3 vectors that represent the borders of the triangle.
+
+    These are vectors to show the *ordering* of the triangles in the data.
+
+    Parameters
+    ----------
+    centers, offsets, triangles : np.ndarray of float
+        The mesh triangulation of the shape's edge path.
+
+    Returns
+    -------
+    borders : np.ndarray of float
+        Positions and projections corresponding to each triangle edge in
+        the triangulation of the path.
+    """
+    borders = np.empty((len(triangles)*3, 2, 2), dtype=centers.dtype)
+    position = 0
+    projection = 1
+    for i, triangle in enumerate(triangles):
+        a, b, c = triangle
+        a1 = centers[a] + offsets[a]
+        b1 = centers[b] + offsets[b]
+        c1 = centers[c] + offsets[c]
+        borders[i * 3, position] = a1
+        borders[i * 3, projection] = (b1 - a1)
+        borders[i * 3 + 1, position] = b1
+        borders[i * 3 + 1, projection] = (c1 - b1)
+        borders[i * 3 + 2, position] = c1
+        borders[i * 3 + 2, projection] = (a1 - c1)
+    return borders
+
+
+@numba.njit
+def generate_face_triangle_borders(vertices, triangles) -> np.ndarray:
+    """For each triangle in mesh generate 3 vectors that represent the borders of the triangle.
+    """
+    borders = np.empty((len(triangles)*3, 2, 2), dtype=vertices.dtype)
+    for i, triangle in enumerate(triangles):
+        a, b, c = triangle
+        a1 = vertices[a]
+        b1 = vertices[b]
+        c1 = vertices[c]
+        borders[i * 3, 0] = a1
+        borders[i * 3, 1] = (b1 - a1)
+        borders[i * 3 + 1, 0] = b1
+        borders[i * 3 + 1, 1] = (c1 - b1)
+        borders[i * 3 + 2, 0] = c1
+        borders[i * 3 + 2, 1] = (a1 - c1)
+    return borders
+
+
+@dataclass
+class Helpers:
+    """Simple class to hold all auxiliary vector data for a shapes layer."""
+    reference_join_points: np.ndarray
+    join_points: np.ndarray
+    direction_vectors: np.ndarray
+    miter_helper_vectors: np.ndarray
+    orthogonal_vectors: np.ndarray
+    miter_vectors: np.ndarray
+    triangles_vectors: np.ndarray
+    face_triangles_vectors: np.ndarray
+
+
+def asdict(dataclass_instance):
+    """Shallow copy version of `dataclasses.asdict`."""
+    return {f.name: getattr(dataclass_instance, f.name)
+            for f in fields(dataclass_instance)}
+
+
+def get_helper_data_from_shapes(shapes_layer: Shapes) -> Helpers:
+    """Function to generate all auxiliary data for a shapes layer."""
+    shapes = shapes_layer._data_view.shapes
+    mesh_list = [(s._edge_vertices, s._edge_offsets, s._edge_triangles)
+                 for s in shapes]
+    path_list = [(s.data, s._closed) for s in shapes]
+    mesh = tuple(np.concatenate(el, axis=0) for el in zip(*mesh_list))
+    face_mesh_list = [(s._face_vertices, s._face_triangles)
+                      for s in shapes
+                      if len(s._face_vertices) > 0]
+    order_vectors_list = [generate_order_vectors(path_, closed)
+                          for path_, closed in path_list]
+
+    helpers = Helpers(
+        reference_join_points=get_reference_edge_triangulation_points(
+                shapes_layer
+                ),
+        join_points=mesh[0] + mesh[1],
+        direction_vectors=np.concatenate(order_vectors_list, axis=0),
+        miter_helper_vectors=np.concatenate(
+                [generate_miter_helper_vectors(o) for o in order_vectors_list],
+                axis=0,
+                ),
+        orthogonal_vectors=np.concatenate(
+                [generate_orthogonal_vectors(o) for o in order_vectors_list],
+                axis=0,
+                ),
+        miter_vectors=np.concatenate(
+                [generate_miter_vectors(m) for m in mesh_list],
+                axis=0,
+                ),
+        triangles_vectors=np.concatenate(
+                [generate_edge_triangle_borders(*m) for m in mesh_list],
+                axis=0,
+                ),
+        face_triangles_vectors=np.concatenate(
+                [generate_face_triangle_borders(*m) for m in face_mesh_list],
+                axis=0,
+                ),
+    )
+
+    return helpers
+
+
+def update_helper_layers(viewer: napari.Viewer, source_layer: Shapes):
+    updated_helpers = get_helper_data_from_shapes(source_layer)
+    for name, data in asdict(updated_helpers).items():
+        viewer.layers[name].data = data
+
+
+def add_helper_layers(viewer: napari.Viewer, source_layer):
+    """Add helper layers to the viewer that track with the source shapes."""
+    helpers = get_helper_data_from_shapes(source_layer)
+    # sizes and colors are hardcoded based on vibes
+    sizes = [0.2, 0.1, 0.1, 0.06, 0.04, 0.05, 0.04, 0.04]
+    colors = ['yellow', 'white', 'red', 'blue',
+              'green', 'yellow', 'pink', 'magenta']
+    for (name, data), size, color in zip(
+            asdict(helpers).items(), sizes, colors
+            ):
+        if 'points' in name:
+            viewer.add_points(data, name=name, size=size, face_color=color)
+        else:  # all other fields are vectors
+            viewer.add_vectors(
+                    data,
+                    name=name,
+                    vector_style='arrow', edge_width=size, edge_color=color,
+                    )
+    source_layer.events.set_data.connect(
+        partial(update_helper_layers, viewer=viewer, source_layer=source_layer)
+    )
+
+
+path = np.array([[0,0], [0,1], [1,1], [1,0]]) * 10
+sparkle = np.array([[1, 1], [10, 0], [1, -1], [0, -10],
+                    [-1, -1], [-10, 0], [-1, 1], [0, 10]])
+fork = np.array([[2, 10], [0, -5], [-2, 10], [-2, -10], [2, -10]])
+
+polygons = [
+    # square
+    generate_regular_polygon(4, radius=1) * 10,
+    # decagon
+    generate_regular_polygon(10, radius=1) * 10 + np.array([[25, 0]]),
+    # triangle
+    generate_regular_polygon(3, radius=1) * 10 + np.array([[0, 25]]),
+    # two sharp prongs
+    fork + np.array([[25, 25]]),
+    # a four-sided star
+    sparkle + np.array([[50, 0]]),
+    # same star, but winding in the opposite direction
+    sparkle[::-1] + np.array([[50, 26]]),
+    # problem shape —
+    # lighting bolt with sharp angles and overlapping edge widths
+    np.array([[10.97627008, 14.30378733],
+              [12.05526752, 10.89766366],
+              [8.47309599, 12.91788226],
+              [8.75174423, 17.83546002],
+              [19.27325521, 7.66883038],
+              [15.83450076, 10.5778984]],
+             ) + np.array([[60, -15]]),
+    ]
+
+paths = [
+    # a simple backwards-c shape
+    path + np.array([[0, 50]]),
+    # an unclosed decagon
+    generate_regular_polygon(10, radius=1) * 10 + np.array([[25, 50]]),
+    # a three-point straight line (tests collinear points in path)
+    np.array([[0, -10], [0, 0], [0, 10]]) + np.array([[50, 50]]),
+    ]
+
+shapes = polygons + paths
+shape_types=['polygon'] * len(polygons) + ['path'] * len(paths)
+
+viewer = napari.Viewer()
+shapes_layer = viewer.add_shapes(shapes, shape_type=shape_types, name='shapes')
+
+add_helper_layers(viewer, source_layer=shapes_layer)
+viewer.layers.selection = {shapes_layer}
+viewer.reset_view()
+
+if __name__ == '__main__':
+    napari.run()

napari/_qt/layer_controls/qt_vectors_controls.py

@@ -125,7 +125,7 @@ def __init__(self, layer) -> None:
         self.widthSpinBox = QDoubleSpinBox()
         self.widthSpinBox.setKeyboardTracking(False)
         self.widthSpinBox.setSingleStep(0.1)
-        self.widthSpinBox.setMinimum(0.1)
+        self.widthSpinBox.setMinimum(0.01)
         self.widthSpinBox.setMaximum(np.inf)
         self.widthSpinBox.setValue(self.layer.edge_width)
         self.widthSpinBox.valueChanged.connect(self.change_width)