diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 0725374..bde39c1 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -26,7 +26,7 @@ jobs:
         run: sudo apt-get update && sudo apt-get install -y libcairo2-dev pkg-config
 
       - name: Install project
-        run: uv sync --extra dev --python ${{ matrix.python-version }}
+        run: uv sync --extra dev --extra image --extra intersecting_cylinders --python ${{ matrix.python-version }}
 
       - name: Run tests
         run: uv run --no-sync pytest -q -m "not slow"
diff --git a/eucare/alternating_flagstones.py b/eucare/alternating_flagstones.py
index 9c60058..5c0d12e 100644
--- a/eucare/alternating_flagstones.py
+++ b/eucare/alternating_flagstones.py
@@ -36,14 +36,15 @@
 
 from collections import defaultdict
 from dataclasses import dataclass, field
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, Any
 
 import numpy as np
+from numpy.typing import NDArray
 
 from .base import angle_to_axis, line_intersection, unit_vector
 from .conway import alternating_flagstone_graph
 from .cutting import cut_out_poly
-from .half import Face, Vertex
+from .half import Face, GeometricHEG, HalfEdge, Vertex
 from .overlap import CREASE_ASSIGNMENT, MOUNTAIN, VALLEY, color_creases
 from .rendering import inset_poly
 
@@ -202,7 +203,7 @@ def _assign_initial_creases(structure: AlternatingFlagstoneStructure) -> None:
 # Metric
 
 
-def _edge_length(h) -> float:
+def _edge_length(h: HalfEdge) -> float:
     return float(np.linalg.norm(h.orig["pos"] - h.dest["pos"]))
 
 
@@ -327,18 +328,18 @@ def optimize_alternating_flagstone(
     #   target_length: indices of the two flagstone-side endpoints whose
     #                  current distance is the *target* length.
     corner_index = {v: i for i, v in enumerate(corners)}
-    connections = []
-    target_length_indices = []
+    connections_list: list[list[int]] = []
+    target_length_indices_list: list[list[int]] = []
     for h in structure.original.halfedges_representing_edges():
         if h.on_border() or h.rev.on_border():
             continue
         va, _ = original_to_corner[(h.orig, h.face)]
         vb, _ = original_to_corner[(h.dest, h.rev.face)]
         vc, _ = original_to_corner[(h.dest, h.face)]
-        connections.append([corner_index[va], corner_index[vb]])
-        target_length_indices.append([corner_index[va], corner_index[vc]])
-    connections = torch.LongTensor(connections)
-    target_length_indices = torch.LongTensor(target_length_indices)
+        connections_list.append([corner_index[va], corner_index[vb]])
+        target_length_indices_list.append([corner_index[va], corner_index[vc]])
+    connections = torch.LongTensor(connections_list)
+    target_length_indices = torch.LongTensor(target_length_indices_list)
 
     star_group_indices = [[corner_index[c] for c in star_groups[star]] for star in star_vertices]
 
@@ -349,18 +350,18 @@ def optimize_alternating_flagstone(
     offsets = torch.nn.Parameter(torch.zeros(n_faces, 2).float())
     star_points = torch.nn.Parameter(initial_star_points.clone())
 
-    def rot_mat(a):
+    def rot_mat(a: Any) -> Any:
         s, c = torch.sin(a), torch.cos(a)
         return torch.stack([c, -s, s, c], dim=-1).view(-1, 2, 2)
 
-    def forward():
+    def forward() -> tuple[Any, Any]:
         rms = rot_mat(angles)
         out = []
         for coords, center, R, off in zip(face_coords0, rotation_centers, rms, offsets):
             out.append(((coords - center) @ R * scale) + center + off)
         return torch.cat(out), star_points
 
-    def loss_fn(face_coords, star_coords):
+    def loss_fn(face_coords: Any, star_coords: Any) -> Any:
         # connection loss
         a = face_coords[connections]
         l = (a[:, 1] - a[:, 0]).pow(2).sum(-1).sqrt()
@@ -388,7 +389,7 @@ def loss_fn(face_coords, star_coords):
     for v, p in zip(corners + star_vertices, position_view):
         v["pos"] = p
 
-    iterator = range(n_steps)
+    iterator: Any = range(n_steps)
     if progress:
         try:  # pragma: no cover
             from tqdm.auto import tqdm
@@ -433,18 +434,18 @@ def loss_fn(face_coords, star_coords):
 # Border processing
 
 
-def _translation_mat(t):
+def _translation_mat(t: NDArray) -> NDArray:
     m = np.eye(3)
     m[:2, 2] = t
     return m
 
 
-def _rotation_mat(a):
+def _rotation_mat(a: float | NDArray) -> NDArray:
     s, c = np.sin(a), np.cos(a)
     return np.array([[c, -s, 0], [s, c, 0], [0, 0, 1]])
 
 
-def _mirror_mat(line):
+def _mirror_mat(line: NDArray | list[NDArray]) -> NDArray:
     """3x3 affine that mirrors the plane across the line ``[p0, p1]``."""
     t1 = _translation_mat(-line[0])
     t2 = _translation_mat(line[0])
@@ -455,7 +456,7 @@ def _mirror_mat(line):
     return t2 @ r2 @ mx @ r1 @ t1
 
 
-def _apply_affine(m, v):
+def _apply_affine(m: NDArray, v: NDArray) -> NDArray:
     return (m @ np.concatenate([v, [1]]))[:2]
 
 
@@ -490,7 +491,9 @@ def extend_border(CP: "EuclideanPositionHEG") -> "EuclideanPositionHEG":
         try:
             v = h.rev.nex.dest
             pos = _apply_affine(_mirror_mat([h.orig["pos"], h.dest["pos"]]), v["pos"])
-            h_new, _ = out.subdivide_face(h.rev.face, h.orig, h.dest)
+            rev_face = h.rev.face
+            assert rev_face is not None
+            h_new, _ = out.subdivide_face(rev_face, h.orig, h.dest)
             h_new[CREASE_ASSIGNMENT] = h_new.rev[CREASE_ASSIGNMENT] = MOUNTAIN
             h_new_2, v_new = out.subdivide_edge(h, pos=pos)
             out.subdivide_face(None, h_new_2.nex.dest, v_new)
@@ -507,7 +510,7 @@ def extend_border(CP: "EuclideanPositionHEG") -> "EuclideanPositionHEG":
 # Cutting twist centres for a foldable CP
 
 
-def _edge_midpoint(h):
+def _edge_midpoint(h: HalfEdge) -> NDArray:
     return np.mean([v["pos"] for v in (h.orig, h.dest)], axis=0)
 
 
@@ -533,7 +536,7 @@ def cut_twist_centres(
     """
     CP_for_folding, (v_map, _, f_map) = structure.CP.copy(return_mappings=True)
 
-    polys = []
+    polys: list[NDArray] = []
     for star, group in structure.star_groups.items():
         star_c = v_map[star]
         if star_c not in CP_for_folding.vertices:
@@ -542,15 +545,15 @@ def cut_twist_centres(
 
         if star_c.on_border():
             h = next(h1 for h1 in star_c.incoming_iter() if h1.on_border()).rev
-            poly = [h.orig["pos"], h.dest["pos"]]
+            poly_pts: list[NDArray] = [h.orig["pos"], h.dest["pos"]]
             while True:
                 h = h.pre.rev
                 if h.on_border():
                     break
                 h = h.pre.rev
-                poly.append(h.dest["pos"])
-            poly.append(_edge_midpoint(h))
-            poly = np.stack(poly)
+                poly_pts.append(h.dest["pos"])
+            poly_pts.append(_edge_midpoint(h))
+            poly = np.stack(poly_pts)
         else:
             poly = np.stack([v["pos"] for v in star_c.vertex_iter() if v in group_c])
         polys.append(poly)
@@ -570,7 +573,7 @@ def cut_twist_centres(
     CP_for_folding.delete_subset(to_delete)
 
     for poly in polys:
-        cut_out_poly(CP_for_folding, inset_poly(poly, inset), delete_outside=True)
+        cut_out_poly(CP_for_folding, np.asarray(inset_poly(list(poly), inset)), delete_outside=True)
     CP_for_folding.recompute_lengths_and_angles()
     return CP_for_folding
 
@@ -579,9 +582,10 @@ def cut_twist_centres(
 # Curved-fold subdivision (for Origami Simulator)
 
 
-def _opposite_triangle_vertex(h):
+def _opposite_triangle_vertex(h: HalfEdge) -> Vertex | None:
     if h.on_border():
         return None
+    assert h.face is not None
     if h.face.order() != 3:
         return None
     return h.nex.dest
@@ -638,10 +642,13 @@ def subdivide_ridges_for_curved_fold(
             )
         new_vs = [G.subdivide_edge(h, pos=p)[1] for p in pts]
 
+        assert h.face is not None
         for vi in new_vs:
             G.subdivide_face(h.face, v1, vi, color_key=(1, 1, 0))
         if not h.rev.on_border():
+            rev_face = h.rev.face
+            assert rev_face is not None
             for vi in new_vs:
-                G.subdivide_face(h.rev.face, vi, v2, color_key=(1, 1, 0))
+                G.subdivide_face(rev_face, vi, v2, color_key=(1, 1, 0))
 
     return G
diff --git a/eucare/base.py b/eucare/base.py
index bcb91a0..7cd0b36 100644
--- a/eucare/base.py
+++ b/eucare/base.py
@@ -7,12 +7,15 @@
 
 from __future__ import annotations
 
-from typing import Callable
+from typing import TYPE_CHECKING, Callable
 
 import numpy as np
 from numba import jit
 from numpy.typing import ArrayLike, NDArray
 
+if TYPE_CHECKING:
+    from .geometries.base import Geometry
+
 pi = np.pi
 tau = 2 * np.pi
 
@@ -54,7 +57,7 @@ def edge_lengths(points: NDArray) -> NDArray:
     return np.linalg.norm(edge_vectors, axis=1)
 
 
-def edge_lengths_and_in_angles(points: NDArray, geometry) -> tuple[list[float], list[float]]:
+def edge_lengths_and_in_angles(points: NDArray, geometry: type[Geometry]) -> tuple[list[float], list[float]]:
     """Return ``(edge_lengths, interior_angles)`` for a polygon under *geometry*.
 
     *geometry* must expose ``distance(p, q)`` and ``angle(p, q, r)`` (see
@@ -128,7 +131,7 @@ def find_affine(line0: NDArray, line1: NDArray) -> NDArray:
     return np.concatenate([linear, offset[None]])
 
 
-def nearest_neighbor(data: NDArray, query: NDArray, return_index: bool = True):
+def nearest_neighbor(data: NDArray, query: NDArray, return_index: bool = True) -> NDArray | tuple[NDArray, int]:
     """Return the nearest point in *data* to *query* (brute-force).
 
     Args:
@@ -138,12 +141,14 @@ def nearest_neighbor(data: NDArray, query: NDArray, return_index: bool = True):
     """
     if len(data.shape) > len(query.shape):
         query = query[None]
-    index = np.argmin(np.linalg.norm(data - query, axis=-1))
-    return data[index], index if return_index else data[index]
+    index = int(np.argmin(np.linalg.norm(data - query, axis=-1)))
+    if return_index:
+        return data[index], index
+    return data[index]
 
 
 @jit(nopython=True)
-def signed_area(pts):
+def signed_area(pts: NDArray) -> float:
     """Signed area of a polygon (positive = CCW). Numba-accelerated.
 
     *pts* must have shape ``(n, 2)``; the polygon is closed implicitly.
@@ -154,7 +159,7 @@ def signed_area(pts):
 
 
 @jit(nopython=True)
-def orientation(pts, eps=0):
+def orientation(pts: NDArray, eps: float = 0) -> int:
     """Return ``+1`` (CCW), ``-1`` (CW), or ``0`` (degenerate) for a polygon.
 
     Numba-accelerated.  *eps* is the absolute-area threshold for degeneracy.
@@ -172,7 +177,7 @@ def euclidean_to_barycentric_map(tri: NDArray) -> Callable[[NDArray], NDArray]:
     """
     tri = np.array(tri, dtype=np.float32)
 
-    def inner(point):
+    def inner(point: NDArray) -> NDArray:
         mat = np.repeat(tri[None, :], 3, axis=0)
         mat[np.eye(3, dtype=bool)] = point
         coords = np.array([signed_area(pts) for pts in mat], dtype=np.float32)
@@ -184,7 +189,7 @@ def inner(point):
 def barycentric_to_euclidean_map(tri: NDArray) -> Callable[[NDArray], NDArray]:
     """Return a function converting barycentric coords back to 2D points w.r.t. *tri*."""
 
-    def inner(barycentric_coords):
+    def inner(barycentric_coords: NDArray) -> NDArray:
         return tri.T @ barycentric_coords
 
     return inner
diff --git a/eucare/classifiers.py b/eucare/classifiers.py
index f1d0fd2..60fc1f3 100755
--- a/eucare/classifiers.py
+++ b/eucare/classifiers.py
@@ -12,16 +12,22 @@
 
 from __future__ import annotations
 
+from typing import Any, Callable
+
 import numpy as np
+from numpy.typing import NDArray
 
 
 class Classifier:
     """Classify items by a hashable index, optionally tracking items and indices per class."""
 
+    saved_items: dict[Any, set[Any]] | None
+    saved_indices: dict[Any, Any] | None
+
     def __init__(self, save_items: bool = False, save_indices: bool = False) -> None:
         super(Classifier, self).__init__()
 
-        self.used_indices = set()
+        self.used_indices: set[Any] = set()
 
         # option to keep track of a dict mapping classes to items
         self.save_items = save_items
@@ -37,16 +43,18 @@ def __init__(self, save_items: bool = False, save_indices: bool = False) -> None
         else:
             self.saved_indices = None
 
-    def _get_index(self, item):
+    def _get_index(self, item: Any) -> Any:
         # the returned 'index' can be any hashable
         raise NotImplementedError
 
-    def classify(self, item):
+    def classify(self, item: Any) -> Any:
         """Return the equivalence class index for ``item`` and update saved items/indices."""
         index = self._get_index(item)
         if self.save_items:
+            assert self.saved_items is not None
             self.saved_items[index] = self.saved_items.get(index, set()).union({item})
         if self.save_indices:
+            assert self.saved_indices is not None
             self.saved_indices[item] = index
         return index
 
@@ -54,13 +62,13 @@ def classify(self, item):
 class CountingClassifier(Classifier):
     """Wrap a classifier to remap its indices to consecutive natural numbers."""
 
-    def __init__(self, other, *super_args, **super_kwargs):
+    def __init__(self, other: Classifier, *super_args: Any, **super_kwargs: Any) -> None:
         super(CountingClassifier, self).__init__(*super_args, **super_kwargs)
         self.non_counting_classifier = other
         self.current_count = 0
-        self.index_to_count = dict()
+        self.index_to_count: dict[Any, int] = dict()
 
-    def _get_index(self, item):
+    def _get_index(self, item: Any) -> int:
         index = self.non_counting_classifier.classify(item)
         if index not in self.index_to_count:
             self.index_to_count[index] = self.current_count
@@ -71,22 +79,22 @@ def _get_index(self, item):
 class RepresentationClassifier(Classifier):
     """Classify items by computing a representation and comparing it against known classes."""
 
-    def __init__(self, *super_args, **super_kwargs):
+    def __init__(self, *super_args: Any, **super_kwargs: Any) -> None:
         super(RepresentationClassifier, self).__init__(*super_args, **super_kwargs)
         self.current_count = 0
-        self.count_to_repr = dict()
+        self.count_to_repr: dict[int, Any] = dict()
         self.represented_first = False
 
-    def _compare_representations(self, query_rep, saved_rep):
+    def _compare_representations(self, query_rep: Any, saved_rep: Any) -> Any:
         return query_rep == saved_rep
 
-    def _represent_item(self, item):
+    def _represent_item(self, item: Any) -> Any:
         return item
 
-    def _represent_query_item(self, item):
+    def _represent_query_item(self, item: Any) -> Any:
         return self._represent_item(item)
 
-    def _get_index(self, item):
+    def _get_index(self, item: Any) -> int:
         query_rep = self._represent_query_item(item)
         if (
             not self.represented_first and self.current_count == 1
@@ -107,17 +115,22 @@ def _get_index(self, item):
 class NestedClassifier(Classifier):
     """Chain multiple classifiers from coarse to fine, producing a tuple index."""
 
-    def __init__(self, coarse_to_fine, *super_args, **super_kwargs):
+    def __init__(
+        self,
+        coarse_to_fine: list[Callable[[], Classifier]],
+        *super_args: Any,
+        **super_kwargs: Any,
+    ) -> None:
         # coarse_to_fine should be a list of classifier classes
         super(NestedClassifier, self).__init__(*super_args, **super_kwargs)
         self.coarse_to_fine = coarse_to_fine
-        self.nested_classfier_dict = dict()
+        self.nested_classfier_dict: dict[Any, Any] = dict()
         self.base_classifier = self.coarse_to_fine[0]()
 
-    def _get_index(self, item):
+    def _get_index(self, item: Any) -> tuple[Any, ...]:
         current_dict = self.nested_classfier_dict
         current_index = self.base_classifier.classify(item)
-        result = (current_index,)
+        result: tuple[Any, ...] = (current_index,)
         for cls in self.coarse_to_fine[1:]:
             if current_index not in current_dict:
                 current_dict[current_index] = dict(classifier=cls(), index_mapping=dict())
@@ -129,13 +142,13 @@ def _get_index(self, item):
         return result
 
 
-def lambda_classifier(func):
+def lambda_classifier(func: Callable[[Any], Any]) -> type[Classifier]:
     """Create a Classifier class that uses the given function as its index."""
 
     class LambdaClassifier(Classifier):
         """Classifier whose index is computed by the wrapped function."""
 
-        def _get_index(self, item):
+        def _get_index(self, item: Any) -> Any:
             return func(item)
 
     return LambdaClassifier
@@ -144,7 +157,7 @@ def _get_index(self, item):
 class LenClassifier(Classifier):
     """Classify items by their length."""
 
-    def _get_index(self, item):
+    def _get_index(self, item: Any) -> int:
         return len(item)
 
 
@@ -154,32 +167,38 @@ def _get_index(self, item):
 class SumClassifier(RepresentationClassifier):
     """Classify items by the sum of their elements (with tolerance)."""
 
-    def _compare_representations(self, query_rep, saved_rep):
+    def _compare_representations(self, query_rep: Any, saved_rep: Any) -> Any:
         return np.all(np.abs(query_rep - saved_rep) < tol)
 
-    def _represent_item(self, item):
+    def _represent_item(self, item: Any) -> Any:
         return np.sum(np.array(item))
 
 
 class UnorderedClassifier(RepresentationClassifier):
     """Classify items by their sorted elements, ignoring order."""
 
-    def _compare_representations(self, query_rep, saved_rep):
+    def _compare_representations(self, query_rep: Any, saved_rep: Any) -> Any:
         return np.all(query_rep == saved_rep)
 
-    def _represent_item(self, item):
+    def _represent_item(self, item: Any) -> NDArray[Any]:
         return np.sort(np.array(item))
 
 
 class CyclicClassifier(RepresentationClassifier):
     """Classify items up to cyclic permutation (and optionally reflection)."""
 
-    def __init__(self, tolerance=tol, allow_flip=False, *super_args, **super_kwargs):
+    def __init__(
+        self,
+        tolerance: float = tol,
+        allow_flip: bool = False,
+        *super_args: Any,
+        **super_kwargs: Any,
+    ) -> None:
         super(CyclicClassifier, self).__init__(*super_args, **super_kwargs)
         self.tolerance = tolerance
         self.allow_flip = allow_flip
 
-    def _compare_representations(self, query_rep, saved_rep):
+    def _compare_representations(self, query_rep: Any, saved_rep: Any) -> Any:
         if query_rep.shape != saved_rep.shape[1:]:
             return False
         if self.tolerance == 0:
@@ -190,7 +209,7 @@ def _compare_representations(self, query_rep, saved_rep):
                 <= self.tolerance
             )
 
-    def _represent_item(self, item):
+    def _represent_item(self, item: Any) -> NDArray[Any]:
         pts = self._represent_query_item(item)
         if not self.allow_flip:
             return np.stack([np.roll(pts, i, axis=0) for i in np.arange(len(pts))])
@@ -203,35 +222,41 @@ def _represent_item(self, item):
                 axis=0,
             )
 
-    def _represent_query_item(self, item):
+    def _represent_query_item(self, item: Any) -> NDArray[Any]:
         return np.array(item)
 
 
 class PreMapClassifier(Classifier):
     """Apply a function to each item before passing it to another classifier."""
 
-    def __init__(self, other, func, *super_args, **super_kwargs):
+    def __init__(
+        self,
+        other: Classifier,
+        func: Callable[[Any], Any],
+        *super_args: Any,
+        **super_kwargs: Any,
+    ) -> None:
         super(PreMapClassifier, self).__init__(*super_args, **super_kwargs)
         self.func = func
         self.other = other
 
-    def _get_index(self, item):
+    def _get_index(self, item: Any) -> Any:
         return self.other.classify(self.func(item))
 
 
-def _face_to_array(f) -> np.ndarray:
+def _face_to_array(f: Any) -> NDArray[Any]:
     """Represent a face as an (n, 2) array of (length, in_angle) pairs along its boundary."""
-    data = []
+    data: list[tuple[Any, Any]] = []
     for e in f.halfedge_iter():
         data.append((e["length"], e["in_angle"]))
         # data.append(np.array(e.orig['pos'], dtype=np.float32))
-    data = np.stack(data)
-    # data -= np.mean(data, axis=0, keepdims=True)
-    # print(data)
-    return data
+    arr = np.stack(data)
+    # arr -= np.mean(arr, axis=0, keepdims=True)
+    # print(arr)
+    return arr
 
 
-def congruency_classifier(allow_flip=False):
+def congruency_classifier(allow_flip: bool = False) -> CountingClassifier:
     """Return a classifier that groups faces by polygon congruence (edge lengths and angles)."""
     return CountingClassifier(
         PreMapClassifier(
@@ -244,9 +269,10 @@ def congruency_classifier(allow_flip=False):
 class AdjacencyClassifier(CyclicClassifier):
     """Classify faces by the cyclic sequence of a given attribute on their neighbors."""
 
-    def __init__(self, key, *super_args, **super_kwargs):
-        super(AdjacencyClassifier, self).__init__(tolerance=0, *super_args, **super_kwargs)
+    def __init__(self, key: str, *super_args: Any, **super_kwargs: Any) -> None:
+        super_kwargs.pop("tolerance", None)
+        super(AdjacencyClassifier, self).__init__(0, *super_args, **super_kwargs)
         self.key = key
 
-    def _represent_query_item(self, item):
+    def _represent_query_item(self, item: Any) -> NDArray[Any]:
         return np.array([(f[self.key] if f is not None else None, item[self.key]) for f in item.face_iter()])
diff --git a/eucare/colorization.py b/eucare/colorization.py
index bd391f4..f3aa0dc 100755
--- a/eucare/colorization.py
+++ b/eucare/colorization.py
@@ -6,15 +6,18 @@
 
 from __future__ import annotations
 
+from typing import Any
+
 from .classifiers import Classifier, congruency_classifier
+from .half import HalfEdgeGraph
 
 
-def colorize(graph, classifier: Classifier, key: str = "color_key") -> None:
+def colorize(graph: HalfEdgeGraph, classifier: Classifier, key: str = "color_key") -> None:
     """Assign ``face[key] = classifier.classify(face)`` for every face in *graph*."""
     for f in graph.faces:
         f[key] = classifier.classify(f)
 
 
-def congruency_colorize(graph, **kwargs) -> None:
+def congruency_colorize(graph: HalfEdgeGraph, **kwargs: Any) -> None:
     """Colour faces by polygon congruence (same edge lengths and interior angles)."""
     colorize(graph, congruency_classifier(), **kwargs)
diff --git a/eucare/conversions.py b/eucare/conversions.py
index 889beb3..f725ac9 100755
--- a/eucare/conversions.py
+++ b/eucare/conversions.py
@@ -12,6 +12,7 @@
 
 import logging
 from copy import copy
+from typing import Any, Literal, overload
 
 import networkx as nx
 import numpy as np
@@ -36,11 +37,28 @@ def _delete_dangling_edges_nx(nx_graph: nx.Graph) -> int:
     return n_deleted
 
 
+@overload
 def EHEG_from_nx(
     nxg: nx.Graph,
-    positions: dict | None = None,
+    positions: dict[Any, Any] | None = ...,
+    *,
+    return_v_lookup: Literal[True],
+) -> tuple[EuclideanPositionHEG, dict[Any, Vertex]]: ...
+
+
+@overload
+def EHEG_from_nx(
+    nxg: nx.Graph,
+    positions: dict[Any, Any] | None = ...,
+    return_v_lookup: Literal[False] = ...,
+) -> EuclideanPositionHEG: ...
+
+
+def EHEG_from_nx(
+    nxg: nx.Graph,
+    positions: dict[Any, Any] | None = None,
     return_v_lookup: bool = False,
-) -> "EuclideanPositionHEG | tuple[EuclideanPositionHEG, dict]":
+) -> EuclideanPositionHEG | tuple[EuclideanPositionHEG, dict[Any, Vertex]]:
     """Convert a planar undirected :class:`networkx.Graph` to an :class:`EuclideanPositionHEG`.
 
     Args:
@@ -71,7 +89,7 @@ def EHEG_from_nx(
         v["pos"] = positions[n]
         v_lookup[n] = v
     result.add_vertices(v_lookup.values())
-    h_lookup = dict()
+    h_lookup: dict[Vertex, dict[Vertex, HalfEdge]] = dict()
     # orig, dest
     for n in nxg.nodes():
         v = v_lookup[n]
@@ -116,8 +134,8 @@ def EHEG_from_nx(
     # detect 'outside' faces which should be None by their orientation
     # it is selected as the one with maximal negative area
     # (area 0 faces might have slightly negative areas due to numerical issues)
-    outside_face = None
-    current_min_area = 0
+    outside_face: Face | None = None
+    current_min_area: float = 0.0
     for f in frozenset(result.faces):
         area = f.area()
         if area < current_min_area:
diff --git a/eucare/conway/factories.py b/eucare/conway/factories.py
index 14d4415..5b2b62e 100644
--- a/eucare/conway/factories.py
+++ b/eucare/conway/factories.py
@@ -14,7 +14,7 @@ def dual_graph() -> GeometricConwayOperator:
     vf = (1, 0)
     v2 = (0, 1)
     ve = (0, 0)
-    G = nx.Graph()
+    G: nx.Graph = nx.Graph()
     nx.add_cycle(G, [v1, vf, v2, ve], delete=True)
     G.add_nodes_from([v1, v2], delete=True)
     G.add_nodes_from([ve], join=True)
@@ -32,7 +32,7 @@ def kis_graph() -> GeometricConwayOperator:
     vf = (1, 0)
     v2 = (0, 1)
     # construct nx graph with delete and join attributes
-    G = nx.Graph()
+    G: nx.Graph = nx.Graph()
     nx.add_cycle(G, [v1, vf, v2])
 
     # construct EHEG and ConwayOperator
@@ -47,7 +47,7 @@ def join_graph() -> GeometricConwayOperator:
     vf = (1, 0)
     v2 = (0, 1)
     # construct nx graph with delete and join attributes
-    G = nx.Graph()
+    G: nx.Graph = nx.Graph()
     nx.add_cycle(G, [v1, vf, v2])
     G.add_edge(v2, v1, delete=True)
     # construct EHEG and ConwayOperator
@@ -61,7 +61,7 @@ def meta_graph() -> GeometricConwayOperator:
     vf = (1, 0)
     v2 = (0, 1)
     ve = (0, 0)
-    G = nx.Graph()
+    G: nx.Graph = nx.Graph()
     nx.add_cycle(G, [v1, vf, v2, ve])
     G.add_nodes_from([v1, v2])
     G.add_nodes_from([ve])
@@ -78,7 +78,7 @@ def ortho_graph() -> GeometricConwayOperator:
     vf = (1, 0)
     v2 = (0, 1)
     ve = (0, 0)
-    G = nx.Graph()
+    G: nx.Graph = nx.Graph()
     nx.add_cycle(G, [v1, vf, v2, ve])
     G.add_edge(ve, vf)
     G.edges[v1, vf]["delete"] = True
@@ -99,11 +99,11 @@ def ambo_graph() -> GeometricConwayOperator:
     v1f = (1 / 2, -1 / 2)
     v2f = (1 / 2, 1 / 2)
     # construct nx graph with delete and join attributes
-    G = nx.Graph()
+    G: nx.Graph = nx.Graph()
     nx.add_cycle(G, [v1, v1f, vf, v2f, v2, v12], delete=True)
     G.add_nodes_from([v1, vf, v2], delete=True)
     G.add_nodes_from([v1f, v2f], join=True)
-    G.add_edges_from([[v12, v1f], [v12, v2f]])
+    G.add_edges_from([(v12, v1f), (v12, v2f)])
 
     # construct EHEG and ConwayOperator
     heg, v_lookup = EHEG_from_nx(G, return_v_lookup=True)
@@ -120,13 +120,13 @@ def goldberg2_graph() -> GeometricConwayOperator:
     v1f = (1 / 2, -1 / 2)
     v2f = (1 / 2, 1 / 2)
     # construct nx graph with delete and join attributes
-    G = nx.Graph()
+    G: nx.Graph = nx.Graph()
     nx.add_cycle(G, [v1, v1f, vf, v2f, v2, v12], delete=True)
     del G.edges[v1, v12]["delete"]
     del G.edges[v12, v2]["delete"]
     G.add_nodes_from([vf], delete=True)
     G.add_nodes_from([v1f, v2f], join=True)
-    G.add_edges_from([[v12, v1f], [v12, v2f]])
+    G.add_edges_from([(v12, v1f), (v12, v2f)])
 
     # construct EHEG and ConwayOperator
     heg, v_lookup = EHEG_from_nx(G, return_v_lookup=True)
@@ -142,12 +142,12 @@ def truncate_graph(t: float = 1 / 2) -> GeometricConwayOperator:
     v1ft = (t / 2, -1 + t / 2)
     v21t = (0, 1 - t)
     v2ft = (t / 2, 1 - t / 2)
-    G = nx.Graph()
+    G: nx.Graph = nx.Graph()
     nx.add_cycle(G, [v1, v1ft, vf, v2ft, v2, v21t, v12t], delete=True)
     del G.edges[v12t, v21t]["delete"]
     G.add_nodes_from([v1, vf, v2], delete=True)
     G.add_nodes_from([v1ft, v2ft], join=True)
-    G.add_edges_from([[v12t, v1ft], [v21t, v2ft]])
+    G.add_edges_from([(v12t, v1ft), (v21t, v2ft)])
 
     # construct EHEG and ConwayOperator
     heg, v_lookup = EHEG_from_nx(G, return_v_lookup=True)
@@ -160,9 +160,9 @@ def gyro_graph(g: tuple[float, float] = (1 / 4, -1 / 4)) -> GeometricConwayOpera
     vf = (1, 0)
     v2 = (0, 1)
     ve = (0, 0)
-    G = nx.Graph()
+    G: nx.Graph = nx.Graph()
     nx.add_cycle(G, [v1, vf, v2, ve], delete=True)
-    G.add_edges_from([[g, ve], [g, v1], [g, vf]])
+    G.add_edges_from([(g, ve), (g, v1), (g, vf)])
     G.add_node(ve, join=True)
 
     # construct EHEG and ConwayOperator
@@ -178,10 +178,10 @@ def starify_graph(t: float = 1 / 3) -> GeometricConwayOperator:
     v0 = (0, 0)
     v1f = (t, -1 + t)
     v2f = (t, 1 - t)
-    G = nx.Graph()
+    G: nx.Graph = nx.Graph()
     nx.add_cycle(G, [v1, v1f, vf, v2f, v2, v0], delete="True")
     G.add_nodes_from([v0], join=True)
-    G.add_edges_from([[v0, v2f], [v2f, v1f]])
+    G.add_edges_from([(v0, v2f), (v2f, v1f)])
 
     # construct EHEG and ConwayOperator
     heg, v_lookup = EHEG_from_nx(G, return_v_lookup=True)
@@ -196,14 +196,14 @@ def alternating_flagstone_graph(t: float = 1 / 3) -> GeometricConwayOperator:
     v0 = (0, 0)
     v1f = (t, -1 + t)
     v2f = (t, 1 - t)
-    G = nx.Graph()
+    G: nx.Graph = nx.Graph()
     nx.add_cycle(G, [v1, v1f, vf, v2f, v2, v0], delete=True)
     del G.edges[v1, v1f]["delete"]
     del G.edges[v2, v2f]["delete"]
     G.add_edge(v2f, v1, color_key=(1, 0, 0))
     G.add_nodes_from([v0], join=True)
     G.add_nodes_from([vf], delete=True)
-    G.add_edges_from([[v0, v2f], [v2f, v1f]])
+    G.add_edges_from([(v0, v2f), (v2f, v1f)])
 
     # construct EHEG and ConwayOperator
     heg, v_lookup = EHEG_from_nx(G, return_v_lookup=True)
@@ -227,15 +227,17 @@ def shrink_rotate_graph(t: float = 1 / 2) -> GeometricConwayOperator:
     v1ft = (t, -1 + t)
     v21t = (0, 1 - t)
     v2ft = (t, 1 - t)
-    G = nx.Graph()
+    G: nx.Graph = nx.Graph()
     nx.add_cycle(G, [v1, v1ft, vf, v2ft, v2, v21t, v12t], delete=True)
     G.add_nodes_from([v1, vf, v2], delete=True)
     G.add_nodes_from([v12t, v21t], join=True)
-    G.add_edges_from([[v12t, v1ft], [v21t, v2ft], [v1ft, v2ft]])
+    G.add_edges_from([(v12t, v1ft), (v21t, v2ft), (v1ft, v2ft)])
 
     # construct EHEG and ConwayOperator
     heg, v_lookup = EHEG_from_nx(G, return_v_lookup=True)
-    v_lookup[vf].get_outgoing_border().rev.face["shrink_rotate"] = True
+    inner_face = v_lookup[vf].get_outgoing_border().rev.face
+    assert inner_face is not None
+    inner_face["shrink_rotate"] = True
     return GeometricConwayOperator(heg, *(v_lookup[v] for v in [v1, vf, v2]))
 
 
@@ -247,10 +249,10 @@ def loft_graph(t: float = 1 / 2) -> GeometricConwayOperator:
     v2 = (0, 1)
     v1ft = (t, -1 + t)
     v2ft = (t, 1 - t)
-    G = nx.Graph()
+    G: nx.Graph = nx.Graph()
     nx.add_cycle(G, [v1, v1ft, v2ft, v2])
     G.add_nodes_from([vf], delete=True)
-    G.add_edges_from([[v1ft, vf], [vf, v2ft]], delete=True)
+    G.add_edges_from([(v1ft, vf), (vf, v2ft)], delete=True)
 
     # construct EHEG and ConwayOperator
     heg, v_lookup = EHEG_from_nx(G, return_v_lookup=True)
@@ -267,12 +269,12 @@ def lace_graph(t: float = 1 / 2, join: bool = False) -> GeometricConwayOperator:
     v1f = (t, -1 + t)
     v2f = (t, 1 - t)
 
-    G = nx.Graph()
+    G: nx.Graph = nx.Graph()
     nx.add_cycle(G, [v1, v1f, vf, v2f, v2], delete=True)
     if not join:
         del G.edges[v1, v2]["delete"]
 
-    G.add_edges_from([[vc, v1], [vc, v2], [vc, v1f], [vc, v2f]])
+    G.add_edges_from([(vc, v1), (vc, v2), (vc, v1f), (vc, v2f)])
     G.add_nodes_from([v1f, v2f], join=True)
     G.add_nodes_from([vf], delete=True)
 
@@ -292,10 +294,10 @@ def expand_graph(t: float = 1 / 2) -> GeometricConwayOperator:
     v1f = (t, -1 + t)
     v2f = (t, 1 - t)
 
-    G = nx.Graph()
+    G: nx.Graph = nx.Graph()
     nx.add_cycle(G, [v2, v21, v12, v1, v1f, vf, v2f], delete=True)
 
-    G.add_edges_from([[v1f, v12], [v2f, v21], [v1f, v2f]])
+    G.add_edges_from([(v1f, v12), (v2f, v21), (v1f, v2f)])
     G.add_nodes_from([v12, v21, v1f, v2f], join=True)
     G.add_nodes_from([vf], delete=True)
 
@@ -321,10 +323,10 @@ def flagstone_pvitelli_graph(t: float = 1 / 4) -> GeometricConwayOperator:
     vN1e = (0, -1 + 3 / 2 * t)
     vN2e = (0, 1 - 3 / 2 * t)
 
-    G = nx.Graph()
+    G: nx.Graph = nx.Graph()
     nx.add_cycle(G, [v2, vL2, vN2e, vN1e, vL1, v1, vL1f, vV1, vf, vV2, vL2f], delete=True)
     nx.add_cycle(G, [vL2, vN2, vN1, vL1, vV1, vV2])
-    G.add_edges_from([[vL1f, vL1], [vL2f, vL2], [vN1e, vN1], [vN2e, vN2], [vN1, vV1], [vN2, vV2]])
+    G.add_edges_from([(vL1f, vL1), (vL2f, vL2), (vN1e, vN1), (vN2e, vN2), (vN1, vV1), (vN2, vV2)])
     G.add_nodes_from([vL1f, vL2f, vN1e, vN2e], join=True)
 
     # label the points which will be joined that are closer to v1
@@ -339,7 +341,9 @@ def flagstone_pvitelli_graph(t: float = 1 / 4) -> GeometricConwayOperator:
     # construct EHEG and ConwayOperator
     heg, v_lookup = EHEG_from_nx(G, return_v_lookup=True)
 
-    v_lookup[vf].get_outgoing_border().rev.face["is_central_polygon"] = True
+    central_face = v_lookup[vf].get_outgoing_border().rev.face
+    assert central_face is not None
+    central_face["is_central_polygon"] = True
 
     return GeometricConwayOperator(heg, *(v_lookup[v] for v in [v1, vf, v2]))
 
diff --git a/eucare/conway/methods.py b/eucare/conway/methods.py
index 9995e81..c3f0a37 100644
--- a/eucare/conway/methods.py
+++ b/eucare/conway/methods.py
@@ -24,6 +24,7 @@
 from __future__ import annotations
 
 import inspect
+from typing import Any, Callable
 
 from ..half import GeometricHEG
 from . import factories
@@ -43,7 +44,7 @@
 _CALL_PARAM_NAMES = frozenset(p.name for p in _CALL_PARAMS)
 
 
-def _shorthand(factory):
+def _shorthand(factory: Callable[..., Any]) -> Callable[..., GeometricHEG]:
     """Build a ``GeometricHEG`` method that applies ``factory(...)`` to ``self``.
 
     The returned method exposes a synthesized signature ``(self, *factory_params,
@@ -52,7 +53,7 @@ def _shorthand(factory):
     ``help()`` and ``?`` even though the method body is only one line.
     """
 
-    def method(self, *args, **kwargs):
+    def method(self: GeometricHEG, *args: Any, **kwargs: Any) -> GeometricHEG:
         call_kwargs = {k: kwargs.pop(k) for k in _CALL_PARAM_NAMES if k in kwargs}
         return factory(*args, **kwargs)(self, **call_kwargs)
 
diff --git a/eucare/conway/operators.py b/eucare/conway/operators.py
index a56c376..f47a2fb 100644
--- a/eucare/conway/operators.py
+++ b/eucare/conway/operators.py
@@ -4,16 +4,23 @@
 
 from collections.abc import Callable
 from copy import copy
+from typing import Any, TypeVar
 
 import numpy as np
+from numpy.typing import NDArray
 
 from eucare.rendering import Rendering
 
 from ..base import euclidean_to_barycentric_map
+from ..geometries.base import Geometry
 from ..half import Face, GeometricHEG, HalfEdge, HalfEdgeGraph, Vertex
 from ..utils import invert_mapping
 
 
+G = TypeVar("G", bound="HalfEdgeGraph")
+G_geom = TypeVar("G_geom", bound="GeometricHEG")
+
+
 class TopologicalConwayOperator:
     """Apply a Conway operator to a half-edge graph by substituting a fundamental domain into each face triangle.
 
@@ -38,14 +45,15 @@ def __init__(self, graph: HalfEdgeGraph, v1: Vertex, vf: Vertex, v2: Vertex) ->
 
     def show(self) -> None:
         """Render the fundamental domain graph for visualization."""
+        assert isinstance(self.graph, GeometricHEG)
         self.graph.show(scale=300, line_width=0.03, render_faces=False)
 
-    def get_tri(self, h: HalfEdge) -> "np.ndarray | None":
+    def get_tri(self, h: HalfEdge) -> NDArray[np.floating[Any]] | None:
         """Return the triangle for half-edge ``h``, or None for purely topological operators."""
         return None
 
     def generate_graph_and_corners(
-        self, tri: "np.ndarray | None", h: HalfEdge | None = None
+        self, tri: NDArray[np.floating[Any]] | None, h: HalfEdge | None = None
     ) -> tuple[HalfEdgeGraph, tuple[Vertex, Vertex, Vertex]]:
         """Return a copy of the fundamental-domain graph and its three corner vertices.
 
@@ -67,12 +75,12 @@ def generate_graph_and_corners(
 
     def __call__(
         self,
-        graph: HalfEdgeGraph,
+        graph: G,
         faces: "set[Face] | Callable[[Face], bool] | None" = None,
         delete_on_border: bool = True,
         delete_inner_border: bool = False,
         copy_graph: bool = False,
-    ) -> HalfEdgeGraph:
+    ) -> G:
         """Apply the operator to ``graph``, optionally restricted to ``faces``.
 
         Args:
@@ -90,10 +98,12 @@ def __call__(
         Returns:
             The (possibly copied) graph after substitution.
         """
+        obj_map: dict[Any, Any] = dict()
         if copy_graph:
-            graph, (v_map, h_map, f_map) = graph.copy(return_mappings=True)
-            v_map, h_map, f_map = [invert_mapping(m) for m in (v_map, h_map, f_map)]
-            obj_map = dict()
+            graph, (v_map_raw, h_map_raw, f_map_raw) = graph.copy(return_mappings=True)
+            v_map = invert_mapping(v_map_raw)
+            h_map = invert_mapping(h_map_raw)
+            f_map = invert_mapping(f_map_raw)
             obj_map.update(v_map)
             obj_map.update(h_map)
             obj_map.update(f_map)
@@ -102,7 +112,6 @@ def __call__(
                     del obj["pre_conway"]
 
         # apply the operator to a set of halfedges in a graph
-        assert isinstance(graph, HalfEdgeGraph)
         if faces is None:
             faces = graph.faces
         elif callable(faces):
@@ -112,21 +121,24 @@ def __call__(
         affected_faces = {h.face for h in halfedges}
         assert None not in affected_faces, "Cannot apply Conway operator to boundary edge"  # Or can we?
         old_halfedges = frozenset(graph.halfedges)
-        v1_out_lookup = dict()
-        v2_out_lookup = dict()
-        vf_lookup = dict()
-        vf_set = set()
+        v1_out_lookup: dict[HalfEdge, HalfEdge] = dict()
+        v2_out_lookup: dict[HalfEdge, HalfEdge] = dict()
+        vf_lookup: dict[HalfEdge, Vertex] = dict()
+        vf_set: set[Vertex] = set()
 
         graphs_and_corners = [self.generate_graph_and_corners(self.get_tri(h), h) for h in halfedges]
 
         for gc, h in zip(graphs_and_corners, halfedges):
             orig_face = h.face
+            assert orig_face is not None
             con_graph, (v1, vf, v2) = gc
 
             # add reference to old face/vertex to new face/vertex
             for new_vertex, old_obj in [(v1, h.dest), (vf, h.face), (v2, h.orig)]:
                 if new_vertex.attributes.get("delete", False):
-                    new_vertex.get_outgoing_border().rev.face["pre_conway"] = old_obj
+                    border_face = new_vertex.get_outgoing_border().rev.face
+                    assert border_face is not None
+                    border_face["pre_conway"] = old_obj
                 else:
                     new_vertex["pre_conway"] = old_obj
 
@@ -183,6 +195,7 @@ def __call__(
                     current = next
             else:
                 if True:  # not delete_on_border or not h.rev.on_border(): #Fixme
+                    assert h.face is not None
                     for k in h.face.halfedge_iter():
                         if (not delete_inner_border) or h.rev.on_border():
                             k["delete"] = False
@@ -190,6 +203,7 @@ def __call__(
                         k.rev["border_delete"] = True
 
                 if h.rev.on_border():
+                    assert h.face is not None
                     graph.delete_face(h.face)
                 else:
                     HalfEdgeGraph.delete_edge(graph, h)
@@ -198,12 +212,14 @@ def __call__(
         # TODO: keep a record of which edges to delete in the end..
         # and a record of all affected vertices to update angles
 
-        to_delete = set()
-        to_process = copy(graph.halfedges)  # this is bad for performance: make everything work locally!
-        to_keep = set()
-        while to_process:
-            h = to_process.pop()
-            to_process.remove(h.rev)
+        to_delete: set[HalfEdge] = set()
+        to_process_set: set[HalfEdge] = copy(
+            graph.halfedges
+        )  # this is bad for performance: make everything work locally!
+        to_keep: set[HalfEdge] = set()
+        while to_process_set:
+            h = to_process_set.pop()
+            to_process_set.remove(h.rev)
             if h.attributes.get("delete", False):
                 # only delete edges if their reverse also wants to be deleted
                 if h.rev.attributes.get("delete", False):
@@ -219,8 +235,8 @@ def __call__(
             to_keep.union(to_delete) == graph.halfedges
         ), f"{graph.halfedges.difference(to_keep.union(to_delete))}, {to_keep.union(to_delete).difference(graph.halfedges)}"
 
-        faces_to_keep = set()
-        faces_to_maybe_remove = set()
+        faces_to_keep: set[Face] = set()
+        faces_to_maybe_remove: set[Face | None] = set()
         while to_keep:
             # find the new faces
             h = to_keep.pop()
@@ -253,6 +269,7 @@ def __call__(
             for e in list(graph.border_edges()):
                 if e in graph.halfedges:
                     if e.rev.attributes.get("border_delete", False):
+                        assert e.rev.face is not None
                         graph.delete_face(e.rev.face)
             # delete dangling faces
             while True:
@@ -283,7 +300,10 @@ def __call__(
 class GeometricConwayOperator(TopologicalConwayOperator):
     """Conway operator that assigns new vertex positions using barycentric coordinate interpolation."""
 
-    def __init__(self, *super_args: object, **super_kwargs: object) -> None:
+    graph: GeometricHEG
+    geometry: type[Geometry] | None
+
+    def __init__(self, graph: GeometricHEG, v1: Vertex, vf: Vertex, v2: Vertex) -> None:
         """Store the fundamental domain and convert its positions to barycentric coordinates.
 
         After construction, every vertex in ``self.graph`` carries barycentric
@@ -291,38 +311,42 @@ def __init__(self, *super_args: object, **super_kwargs: object) -> None:
         re-projected to Euclidean coordinates per target triangle inside
         :meth:`generate_graph_and_corners`.
         """
-        super(GeometricConwayOperator, self).__init__(*super_args, **super_kwargs)
+        super(GeometricConwayOperator, self).__init__(graph, v1, vf, v2)
         # convert euclidean to barycentric coordinates
         to_barycentric = euclidean_to_barycentric_map(np.array([self.v1["pos"], self.vf["pos"], self.v2["pos"]]))
         for v in self.graph.vertices:
             v["pos"] = to_barycentric(v["pos"])
         self.geometry = None
 
-    def get_tri(self, h: HalfEdge) -> np.ndarray:
+    def get_tri(self, h: HalfEdge) -> NDArray[np.floating[Any]]:
         """Return the triangle ``(h.dest, face midpoint, h.orig)`` for half-edge ``h``."""
+        assert h.face is not None
+        assert self.geometry is not None
         midpoint = h.face.get(
             "midpoint", self.geometry.center_of_mass(np.stack([v["pos"] for v in h.face.vertex_iter()]))
         )
         return np.array([h.dest["pos"], midpoint, h.orig["pos"]])
 
     def generate_graph_and_corners(
-        self, tri: np.ndarray, h: HalfEdge | None = None
+        self, tri: NDArray[np.floating[Any]] | None, h: HalfEdge | None = None
     ) -> tuple[HalfEdgeGraph, tuple[Vertex, Vertex, Vertex]]:
         """Return a copy of the domain with positions mapped from barycentric to Euclidean coordinates."""
         result, corners = super(GeometricConwayOperator, self).generate_graph_and_corners(tri, h)
 
+        assert tri is not None
+        assert self.geometry is not None
         to_euclidean = self.geometry.barycentric_to_euclidean_map(tri)
         # this could be vectorized
         for v in result.vertices:
             v["pos"] = to_euclidean(v["pos"])
         return result, corners
 
-    def __call__(
+    def __call__(  # type: ignore[override]
         self,
-        graph: GeometricHEG,
+        graph: G_geom,
         recompute_lengths_and_angles: bool = True,
-        **kwargs: object,
-    ) -> GeometricHEG:
+        **kwargs: Any,
+    ) -> G_geom:
         """Apply the geometric operator to ``graph``.
 
         Args:
@@ -334,7 +358,6 @@ def __call__(
         Returns:
             The transformed graph.
         """
-        assert isinstance(graph, GeometricHEG)
         self.geometry = graph.geometry
         result = super().__call__(graph, **kwargs)
         if recompute_lengths_and_angles:
@@ -343,14 +366,14 @@ def __call__(
 
     #: Canonical reference triangle for visualizing the fundamental domain.
     #: Matches the ``v1, vf, v2`` layout used by every factory in :mod:`eucare.conway.factories`.
-    _SHOW_REFERENCE_TRIANGLE: np.ndarray = np.array([[0.0, -1.0], [1.0, 0.0], [0.0, 1.0]])
+    _SHOW_REFERENCE_TRIANGLE: NDArray[np.floating[Any]] = np.array([[0.0, -1.0], [1.0, 0.0], [0.0, 1.0]])
 
     def get_fundamental_domain_graph_to_render(
         self,
         delete_color: tuple[float, float, float] = (0.85, 0.15, 0.15),
         join_color: tuple[float, float, float] = (0.15, 0.65, 0.20),
         keep_color: tuple[float, float, float] = (0.30, 0.30, 0.30),
-    ) -> HalfEdgeGraph:
+    ) -> GeometricHEG:
         """Build a copy of the fundamental-domain graph, colouring elements by their role.
 
         Vertices and edges flagged as ``delete`` are coloured *delete_color*
@@ -389,18 +412,18 @@ def get_fundamental_domain_graph_to_render(
                 h["color_key"] = keep_color
         return graph_copy
 
-    def render(self, **show_kwargs: object) -> Rendering:
+    def render(self, **show_kwargs: Any) -> Rendering:
         """Render the fundamental domain graph with styling from :meth:`get_fundamental_domain_graph_to_render`."""
 
         fundamental_domain_graph = self.get_fundamental_domain_graph_to_render()
-        kwargs = dict(
+        kwargs: dict[str, Any] = dict(
             line_width="50%",
             face_inset=0,
         )
         kwargs.update(show_kwargs)
         return fundamental_domain_graph.render(**kwargs)
 
-    def show(self, **kwargs: object) -> None:
+    def show(self, **kwargs: Any) -> None:
         """Render the fundamental domain graph with styling from :meth:`get_fundamental_domain_graph_to_render` and display it."""
         rendering = self.render(**kwargs)
         rendering.show()
diff --git a/eucare/cutting.py b/eucare/cutting.py
index 9b06a83..e9b669c 100755
--- a/eucare/cutting.py
+++ b/eucare/cutting.py
@@ -15,11 +15,12 @@
 from __future__ import annotations
 
 from collections import defaultdict
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, Any
 
 import numba
 import numpy as np
 from numba import jit, njit
+from numpy.typing import NDArray
 
 from eucare.half import rotate_by
 
@@ -34,7 +35,7 @@
 
 
 @jit(nopython=True)
-def pointinpolygon(x: float, y: float, poly: np.ndarray) -> bool:
+def pointinpolygon(x: float, y: float, poly: NDArray[Any]) -> bool:
     """Numba-jitted point-in-polygon test (ray casting).
 
     Args:
@@ -66,7 +67,7 @@ def pointinpolygon(x: float, y: float, poly: np.ndarray) -> bool:
 
 
 @njit(parallel=True)
-def parallelpointinpolygon(points: np.ndarray, polygon: np.ndarray) -> np.ndarray:
+def parallelpointinpolygon(points: NDArray[Any], polygon: NDArray[Any]) -> NDArray[Any]:
     """Vectorised :func:`pointinpolygon` over a batch of query points."""
     D = np.empty(len(points), dtype=numba.boolean)
     for i in numba.prange(0, len(D)):
@@ -75,7 +76,7 @@ def parallelpointinpolygon(points: np.ndarray, polygon: np.ndarray) -> np.ndarra
 
 
 @numba.jit(nopython=True)
-def polygon_line_segment_intersections(poly: np.ndarray, line_segment: np.ndarray, eps: float = 1e-12) -> list:
+def polygon_line_segment_intersections(poly: NDArray[Any], line_segment: NDArray[Any], eps: float = 1e-12) -> list[Any]:
     """Return the list of intersections of *line_segment* with each edge of closed polygon *poly*."""
     intersections = []
     p1 = poly[-1]
@@ -86,7 +87,9 @@ def polygon_line_segment_intersections(poly: np.ndarray, line_segment: np.ndarra
     return intersections
 
 
-def get_potential_intersections_2(segments1: np.ndarray, segments2: np.ndarray, epsilon: float = 1e-12) -> list:
+def get_potential_intersections_2(
+    segments1: NDArray[Any], segments2: NDArray[Any], epsilon: float = 1e-12
+) -> list[tuple[int, int]]:
     """Bounding-box prefilter: return ``(i, j)`` pairs whose AABBs overlap.
 
     Used as a pre-pass before the expensive exact intersection test.  Sweep-line
@@ -107,22 +110,22 @@ def get_potential_intersections_2(segments1: np.ndarray, segments2: np.ndarray,
     x_coords.sort(axis=1)
 
     # create tuples (position, index, is_start)
-    points = [(s[0][0] - epsilon, i, 1) for i, s in enumerate(segments)] + [
+    points_list: list[tuple[Any, int, int]] = [(s[0][0] - epsilon, i, 1) for i, s in enumerate(segments)] + [
         (s[1][0], i, 0) for i, s in enumerate(segments)
     ]
-    points = np.array(points, dtype=tuple)
+    points = np.array(points_list, dtype=tuple)
 
     # sort by first coordinate
     points = points[np.argsort(points[:, 0])]
-    active_labels_1 = set()
-    active_labels_2 = set()
-    possibly_intersecting = list()
+    active_labels_1: set[int] = set()
+    active_labels_2: set[int] = set()
+    possibly_intersecting: list[tuple[int, int]] = list()
     for i, is_start in points[:, 1:]:
 
         if i < ns1:  # segment in segments1
             if is_start:
                 for j in active_labels_2:
-                    if intervals_overlapping(segments[i, :, 1], segments[j, :, 1]):
+                    if intervals_overlapping(tuple(segments[i, :, 1]), tuple(segments[j, :, 1])):
                         possibly_intersecting.append((i, j - ns1))
                 active_labels_1.add(i)
             else:
@@ -130,7 +133,7 @@ def get_potential_intersections_2(segments1: np.ndarray, segments2: np.ndarray,
         else:  # segment in segments1
             if is_start:
                 for j in active_labels_1:
-                    if intervals_overlapping(segments[i, :, 1], segments[j, :, 1]):
+                    if intervals_overlapping(tuple(segments[i, :, 1]), tuple(segments[j, :, 1])):
                         possibly_intersecting.append((j, i - ns1))
                 active_labels_2.add(i)
             else:
@@ -138,7 +141,9 @@ def get_potential_intersections_2(segments1: np.ndarray, segments2: np.ndarray,
     return possibly_intersecting
 
 
-def get_ordered_crossings(segments1: np.ndarray, segments2: np.ndarray, eps: float = 1e-10) -> tuple:
+def get_ordered_crossings(
+    segments1: NDArray[Any], segments2: NDArray[Any], eps: float = 1e-10
+) -> tuple[NDArray[Any], list[Any], list[Any], list[set[int]], list[set[int]]]:
     """Compute exact intersections between two segment sets and group near-duplicates.
 
     Args:
@@ -154,26 +159,26 @@ def get_ordered_crossings(segments1: np.ndarray, segments2: np.ndarray, eps: flo
         ``crossings_to_segmentsX[k]`` is the set of segment indices in set X
         contributing to crossing ``k``.
     """
-    crossings = []
+    crossings_list: list[Any] = []
     # This will be a list mapping the index of a crossing to the pair (i, j) of indices of the corresponding edges.
-    crossings_to_edges = []
+    crossings_to_edges: list[tuple[int, int]] = []
 
     for i, j in get_potential_intersections_2(segments1, segments2, epsilon=eps):
         l1, l2 = segments1[i], segments2[j]
         intersections = line_segment_intersections(l1, l2, eps=eps)
         if not intersections:
             continue
-        crossings.extend(intersections)
+        crossings_list.extend(intersections)
         for _ in range(len(intersections)):
             crossings_to_edges.append((i, j))
-    crossings = np.array(crossings)
+    crossings = np.array(crossings_list)
     if len(crossings) == 0:
         return (
             np.zeros(0, dtype=np.int32),
             [],
             [],
-            [[] for _ in range(len(segments1))],
-            [[] for _ in range(len(segments2))],
+            [set() for _ in range(len(segments1))],
+            [set() for _ in range(len(segments2))],
         )
 
         # ------ group closeby crossings ------
@@ -185,10 +190,10 @@ def get_ordered_crossings(segments1: np.ndarray, segments2: np.ndarray, eps: flo
 
     # construct an array mapping crossings to all involved edges,
     # and on mapping edges to all crossings they are involved in.
-    filtered_crossings_to_segments1 = [set() for i in range(n_filtered_crossings)]
-    filtered_crossings_to_segments2 = [set() for i in range(n_filtered_crossings)]
-    segments1_to_crossings = [set() for i in range(len(segments1))]
-    segments2_to_crossings = [set() for i in range(len(segments2))]
+    filtered_crossings_to_segments1: list[set[int]] = [set() for i in range(n_filtered_crossings)]
+    filtered_crossings_to_segments2: list[set[int]] = [set() for i in range(n_filtered_crossings)]
+    segments1_to_crossings: list[Any] = [set() for i in range(len(segments1))]
+    segments2_to_crossings: list[Any] = [set() for i in range(len(segments2))]
     for i, edge_ids in enumerate(crossings_to_edges):
         i_filtered = clustering[i]
         filtered_crossings_to_segments1[i_filtered].add(edge_ids[0])
@@ -197,7 +202,7 @@ def get_ordered_crossings(segments1: np.ndarray, segments2: np.ndarray, eps: flo
         segments2_to_crossings[edge_ids[1]].add(clustering[i])
 
     # ------ get crossing orders ------
-    def order_crossings(segments, segments_to_crossings):
+    def order_crossings(segments: NDArray[Any], segments_to_crossings: list[Any]) -> None:
         segments_dirs = segments[:, 1] - segments[:, 0]
         segments_dirs /= np.linalg.norm(segments_dirs, axis=1, keepdims=True)
         for i, crossing_ids in enumerate(segments_to_crossings):
@@ -266,7 +271,7 @@ def cut_out_poly(
     # numba-jit'd ``line_segment_intersections`` requires a uniform float64
     # dtype across both segment arrays; cast defensively here so that callers
     # that produced positions via float32 paths (e.g. torch optimisers) work.
-    poly_segments = np.stack(list(rotate_by(poly, (0, 1)))).astype(np.float64)
+    poly_segments = np.stack(list(rotate_by(list(poly), (0, 1)))).astype(np.float64)
     es = list(G.halfedges_representing_edges())
     edge_segments = np.array([[e.orig["pos"], e.dest["pos"]] for e in es], dtype=np.float64)
     (
@@ -337,7 +342,7 @@ def cut_out_poly(
                             or np.all(
                                 parallelpointinpolygon(
                                     np.array(corners),
-                                    np.asarray(inset_poly(np.stack([vf["pos"] for vf in f.vertex_iter()]), -eps)),
+                                    np.asarray(inset_poly([vf["pos"] for vf in f.vertex_iter()], -eps)),
                                 )
                             )
                         )
@@ -397,11 +402,7 @@ def cut_out_poly(
                             or np.all(
                                 parallelpointinpolygon(
                                     np.array(corners),
-                                    np.asarray(
-                                        inset_poly(
-                                            np.stack([vf["pos"] for vf in reversed(list(f.vertex_iter()))]), -eps
-                                        )
-                                    ),
+                                    np.asarray(inset_poly([vf["pos"] for vf in reversed(list(f.vertex_iter()))], -eps)),
                                 )
                             )
                         )
diff --git a/eucare/example_graphs.py b/eucare/example_graphs.py
index 31a0000..b8efaad 100755
--- a/eucare/example_graphs.py
+++ b/eucare/example_graphs.py
@@ -3,8 +3,10 @@
 from __future__ import annotations
 
 import logging
+from typing import Any, Callable, Literal, Sequence, overload
 
 import numpy as np
+from numpy.typing import NDArray
 from sympy import N, elliptic_f
 
 from .example_tilesets import curved_zip, pgg_2x
@@ -17,7 +19,9 @@
 # TODO https://en.wikipedia.org/wiki/File:Planar_Fractalizing_Truncated_Hexagonal_Tiling_II.png
 
 
-def get_edge_with(graph: HalfEdgeGraph, func: "callable | None" = None, on_border: bool = False) -> HalfEdge:
+def get_edge_with(
+    graph: HalfEdgeGraph, func: Callable[[HalfEdge], bool] | None = None, on_border: bool = False
+) -> HalfEdge:
     """Return the first half-edge satisfying ``func``, optionally restricted to border edges.
 
     Args:
@@ -37,7 +41,9 @@ def get_edge_with(graph: HalfEdgeGraph, func: "callable | None" = None, on_borde
     raise LookupError("Cannot find edge with requested property")
 
 
-def get_vertex_with(graph: HalfEdgeGraph, func: "callable | None" = None, on_border: bool = False) -> Vertex:
+def get_vertex_with(
+    graph: HalfEdgeGraph, func: Callable[[Vertex], bool] | None = None, on_border: bool = False
+) -> Vertex:
     """Return the first vertex satisfying ``func``, optionally restricted to border vertices.
 
     Args:
@@ -61,8 +67,8 @@ def rosette(n: int = 8) -> EuclideanPositionHEG:
     """Construct a rosette pattern from ``n`` rhombus tiles around a central vertex."""
     assert isinstance(n, int)
     alpha = 2 * np.pi / n
-    G, edgedict = RhombusTile(alpha).make_graph(add_positions=True)
-    G = EuclideanPositionHEG(other=G)
+    base_graph, edgedict = RhombusTile(alpha).make_graph(add_positions=True)
+    G = EuclideanPositionHEG(other=base_graph)
     v = edgedict[0].dest
     while v.on_border():
         RhombusTile(alpha).attach_instruction(0)(G, v.get_outgoing_border())
@@ -98,16 +104,39 @@ def complete_closest_vertices(G: GeometricHEG, eps: float = 1e-6) -> None:
     assert isinstance(G, GeometricHEG)
     vertex_dists = {e.orig: G.geometry.distance_to_origin(e.orig["pos"]) for e in G.border_edge_iter()}
     d_min = np.min(list(vertex_dists.values()))
-    [complete_vertex(G, v) for v, d in vertex_dists.items() if d - d_min < eps and v.on_border()]
+    for v, d in vertex_dists.items():
+        if d - d_min < eps and v.on_border():
+            complete_vertex(G, v)
+
+
+@overload
+def from_tiles(
+    tiles: Sequence[ProtoTile],
+    rings: int = ...,
+    vertex_based: bool = ...,
+    base_tile: "int | ProtoTile | HalfEdgeGraph" = ...,
+    add_positions: Literal[True] = ...,
+) -> GeometricHEG: ...
+
+
+@overload
+def from_tiles(
+    tiles: Sequence[ProtoTile],
+    rings: int = ...,
+    vertex_based: bool = ...,
+    base_tile: "int | ProtoTile | HalfEdgeGraph" = ...,
+    *,
+    add_positions: Literal[False],
+) -> InAngleHEG: ...
 
 
 def from_tiles(
-    tiles: list[ProtoTile],
+    tiles: Sequence[ProtoTile],
     rings: int = 2,
     vertex_based: bool = True,
     base_tile: "int | ProtoTile | HalfEdgeGraph" = -1,
     add_positions: bool = True,
-) -> HalfEdgeGraph:
+) -> InAngleHEG:
     """Grow a tiling from a list of proto-tiles by expanding for the given number of rings.
 
     Args:
@@ -127,6 +156,7 @@ def from_tiles(
         base_tile = tiles[base_tile]
     if isinstance(base_tile, ProtoTile):
         base_tile = base_tile.make_graph(add_positions=add_positions)[0]
+    G: InAngleHEG
     if add_positions:
         assert (
             len({tile.geometry for tile in tiles}) == 1
@@ -147,13 +177,13 @@ def from_tiles(
     return G
 
 
-def pgg_2x_tiling(rings: int = 15) -> HalfEdgeGraph:
+def pgg_2x_tiling(rings: int = 15) -> InAngleHEG:
     """Construct a pgg wallpaper group tiling with the given number of rings."""
-    tiles = pgg_2x()
+    tiles = list(pgg_2x())
     return from_tiles(tiles, rings)
 
 
-def kised_soccer_ball() -> HalfEdgeGraph:
+def kised_soccer_ball() -> GeometricHEG:
     """Construct a kised soccer ball (icosahedron variant) on the sphere."""
     from eucare.conway import kis_graph
 
@@ -183,23 +213,23 @@ def hyperbolic_square_graph(
     """
     import eucare as ec
 
-    def wrapped_elliptic_f(z, m):
+    def wrapped_elliptic_f(z: complex, m: complex) -> complex:
         return complex(N(elliptic_f(z, m)))
 
-    def complex_to_array(z):
+    def complex_to_array(z: complex) -> NDArray[Any]:
         return np.array([z.real, z.imag])
 
-    def disk_to_square(z, w=1):
+    def disk_to_square(z: complex, w: complex = 1) -> complex:
         return np.sqrt(1j) * wrapped_elliptic_f(np.arcsin(w * z), -1)
         # return np.sqrt(2) * wrapped_elliptic_f(np.arcsin(np.sqrt(z+1)), np.sqrt(2)/2)
 
-    def disk_to_halfplane(z):
+    def disk_to_halfplane(z: complex) -> complex:
         return (z + 1j) / (1j * z + 1)
 
     if G is None:
-        G = ec.example_graphs.from_tiles(ec.example_graphs.curved_platonic(7, 3), 1)
+        G = ec.example_graphs.from_tiles(ec.example_tilesets.curved_platonic(7, 3), 1)
 
-    def map_to_square(G):
+    def map_to_square(G: GeometricHEG) -> GeometricHEG:
         G_square = G.copy()
         G_square.geometry = ec.geometries.EuclideanGeometry
         for v in G_square.vertices:
@@ -210,7 +240,7 @@ def map_to_square(G):
 
         return G_square
 
-    nv_before = None
+    nv_before: int | None = None
     nv_after = G.order
     while nv_before != nv_after:
         nv_before = nv_after
diff --git a/eucare/example_tilesets.py b/eucare/example_tilesets.py
index 26c292c..2b5523e 100755
--- a/eucare/example_tilesets.py
+++ b/eucare/example_tilesets.py
@@ -5,6 +5,7 @@
 import numpy as np
 
 from .geometries import EuclideanGeometry, PoincareDiskModel, SphereModel
+from .geometries.base import Geometry
 from .prototiles import PolygonalProtoTile, ProtoTile, RegularEuclideanTile, RegularProtoTile
 
 
@@ -21,7 +22,7 @@
 #         raise NotImplementedError
 
 
-def align_tiles(tile1: PolygonalProtoTile, label1, tile2: PolygonalProtoTile, label2) -> None:
+def align_tiles(tile1: PolygonalProtoTile, label1: object, tile2: PolygonalProtoTile, label2: object) -> None:
     """Set mutual gluing instructions between two tiles along the given edge labels."""
     tile1.edge_instructions[label1] = tile2.attach_instruction(label2)
     tile2.edge_instructions[label2] = tile1.attach_instruction(label1)
@@ -150,7 +151,7 @@ def u2_4_6_12__3_4_6_4() -> tuple[RegularEuclideanTile, ...]:
 # cairo_sq_B.edge_instructions['a'] = attatch_tile_instruction(cairo_tri, 1)
 
 
-def archimedean_vertex_to_geometry(face_orders: "list[int] | tuple[int, ...]") -> type:
+def archimedean_vertex_to_geometry(face_orders: "list[int] | tuple[int, ...]") -> type[Geometry]:
     """Determine the geometry (Euclidean, spherical, or hyperbolic) from vertex face orders.
 
     Args:
diff --git a/eucare/flat_foldable.py b/eucare/flat_foldable.py
index e1b5c55..c3b8f5f 100644
--- a/eucare/flat_foldable.py
+++ b/eucare/flat_foldable.py
@@ -14,6 +14,8 @@
 
 from __future__ import annotations
 
+from typing import Iterable
+
 import numpy as np
 
 from .half import HalfEdgeGraph, Vertex
@@ -35,12 +37,14 @@ def kawasaki_sum(v: Vertex) -> float:
     return np.sum(((angles + 2 * np.pi) % (2 * np.pi)) * (-1) ** np.arange(len(angles)))
 
 
-def max_kawasaki_sum(vertices) -> float:
+def max_kawasaki_sum(vertices: HalfEdgeGraph | Iterable[Vertex]) -> float:
     """Return the largest absolute Kawasaki sum over interior vertices.
 
     *vertices* may be a :class:`HalfEdgeGraph` (interior vertices are taken
     automatically) or any iterable of vertices.
     """
     if isinstance(vertices, HalfEdgeGraph):
-        vertices = [v for v in vertices.vertices if not v.on_border()]
-    return np.max([kawasaki_sum(v) for v in vertices])
+        vertex_list: Iterable[Vertex] = [v for v in vertices.vertices if not v.on_border()]
+    else:
+        vertex_list = vertices
+    return float(np.max([kawasaki_sum(v) for v in vertex_list]))
diff --git a/eucare/geometries/base.py b/eucare/geometries/base.py
index 2785ae8..0bb5b5e 100755
--- a/eucare/geometries/base.py
+++ b/eucare/geometries/base.py
@@ -14,18 +14,23 @@
 from __future__ import annotations
 
 from collections import Counter
+from typing import Any, Callable, Iterable
 
 import numpy as np
+from numpy.typing import NDArray
 from scipy.optimize import root_scalar
 
+Point = Any
+Transform = Callable[[Point], Point]
 
-def root_return(func):
+
+def root_return(func: Callable[..., Any]) -> Callable[..., float]:
     """Decorator that extracts the root from a root_scalar result or raises on failure."""
 
-    def inner(*args, **kwargs):
+    def inner(*args: Any, **kwargs: Any) -> float:
         result = func(*args, **kwargs)
         if result.converged:
-            return result.root
+            return float(result.root)
         else:
             raise ValueError("No root was found")
 
@@ -36,91 +41,91 @@ class Geometry:
     """Base class for 2D geometries (Euclidean, hyperbolic, spherical)."""
 
     @classmethod
-    def origin(cls):
+    def origin(cls) -> Point:
         """Return the origin of the geometry."""
         raise NotImplementedError
 
     @classmethod
-    def translation(cls, p1, p2):
+    def translation(cls, p1: Point, p2: Point) -> Transform:
         """Return a callable that translates points from p1 to p2."""
         raise NotImplementedError
 
     @classmethod
-    def _rotate_around_origin(cls, a1):
+    def _rotate_around_origin(cls, a1: float) -> Transform:
         """Return a callable that rotates points by angle a1 around the origin."""
         raise NotImplementedError
 
     @classmethod
-    def center_of_mass(cls, points, masses=None):
+    def center_of_mass(cls, points: NDArray[Any], masses: NDArray[Any] | None = None) -> Point:
         """Compute the center of mass of point masses at the given positions."""
         raise NotImplementedError
 
     @classmethod
-    def distance_to_origin(cls, p):
+    def distance_to_origin(cls, p: Point) -> float:
         """Compute the distance from point p to the origin."""
         raise NotImplementedError
 
     @classmethod
-    def angle_to_axis(cls, p):
+    def angle_to_axis(cls, p: Point) -> float:
         """Compute the angle from the standard ray to the ray from the origin through p."""
         raise NotImplementedError
 
     @classmethod
-    def point_along_axis(cls, x):
+    def point_along_axis(cls, x: float) -> Point:
         """Return the point on the standard axis at signed distance x from the origin."""
         raise NotImplementedError
 
     @classmethod
-    def to_euclidean(cls, pts):
+    def to_euclidean(cls, pts: NDArray[Any]) -> NDArray[Any]:
         """Convert points from this geometry's representation to Euclidean 2D coordinates."""
         raise NotImplementedError
 
     @classmethod
-    def invert(cls, p):
+    def invert(cls, p: Point) -> Point:
         """Return the point diametrically opposite to p through the origin."""
         return cls.translation(p, cls.origin())(cls.origin())
 
     @classmethod
-    def rotation(cls, p1, a1):
+    def rotation(cls, p1: Point, a1: float) -> Transform:
         """Return a callable that rotates points by angle a1 around point p1."""
         t1 = cls.translation(p1, cls.origin())
         r = cls._rotate_around_origin(a1)
         t2 = cls.translation(cls.origin(), p1)
 
-        def rotate(pts):
+        def rotate(pts: Point) -> Point:
             return t2(r(t1(pts)))
 
         return rotate
 
     @classmethod
-    def angle(cls, p1, p2, p3):
+    def angle(cls, p1: Point, p2: Point, p3: Point) -> float:
         """Return the angle from p1 to p3 with apex at p2."""
         p1, p3 = cls.translation(p2, cls.origin())(np.array([p1, p3]))
         a1, a3 = cls.angle_to_axis(p1), cls.angle_to_axis(p3)
-        return (a1 - a3) % (2 * np.pi)
+        return float((a1 - a3) % (2 * np.pi))
 
     @classmethod
-    def to_polar(cls, p):
+    def to_polar(cls, p: Point) -> tuple[float, float]:
         """Compute polar coordinates (distance, angle) of point p."""
         return cls.distance_to_origin(p), cls.angle_to_axis(p)
 
     @classmethod
-    def distance(cls, p1, p2):
+    def distance(cls, p1: Point, p2: Point) -> float:
         """Compute the distance between points p1 and p2."""
         return cls.distance_to_origin(cls.translation(p1, cls.origin())(p2))
 
     @classmethod
-    def from_polar(cls, r, a):
+    def from_polar(cls, r: float, a: float) -> Point:
         """Construct a point from polar coordinates (radius r, angle a)."""
         return cls.rotation(cls.origin(), a)(cls.point_along_axis(r))
 
     @classmethod
-    def unit_vector(cls, a):
+    def unit_vector(cls, a: float) -> Point:
         """Return the point at unit distance from the origin at angle a."""
         return cls.from_polar(r=1, a=a)
 
     @classmethod
-    def construct_next_poly_point(cls, a, b, angle, length):
+    def construct_next_poly_point(cls, a: Point, b: Point, angle: float, length: float) -> Point:
         """Construct point c such that angle(a, b, c) = angle and dist(b, c) = length."""
         a0 = cls.translation(b, cls.origin())(a)
         c0 = cls.from_polar(length, cls.angle_to_axis(a0) - angle)
@@ -128,33 +133,35 @@ def construct_next_poly_point(cls, a, b, angle, length):
         return c
 
     @classmethod
-    def regular_poly_in_angle(cls, n, r):
+    def regular_poly_in_angle(cls, n: int, r: float) -> float:
         """Compute the interior angle of a regular n-gon with circumradius r."""
         return 2 * cls.angle(cls.from_polar(r, -np.pi / n), cls.from_polar(r, np.pi / n), cls.origin())
 
     @classmethod
-    def regular_poly_side_length(cls, n, r):
+    def regular_poly_side_length(cls, n: int, r: float) -> float:
         """Compute the side length of a regular n-gon with circumradius r."""
         return cls.distance(cls.from_polar(r, -np.pi / n), cls.from_polar(r, np.pi / n))
 
     @classmethod
-    def platonic_side_length(cls, n, k):
+    def platonic_side_length(cls, n: int, k: int) -> float:
         """Find the side length of a regular n-gon with interior angle 2*pi/k."""
         raise NotImplementedError
 
     @classmethod
     @root_return
-    def platonic_side_length_to_radius(cls, n, l):
+    def platonic_side_length_to_radius(cls, n: int, l: float) -> Any:
         """Find the circumradius of a regular n-gon with side length l."""
         return root_scalar(lambda r: cls.regular_poly_side_length(n, r) - l, x0=0.1, x1=0.01)
 
     @classmethod
     @root_return
-    def archimedean_side_length(cls, faces_around_corner, **archimedean_side_length_root_kwargs):
+    def archimedean_side_length(
+        cls, faces_around_corner: Iterable[int], **archimedean_side_length_root_kwargs: Any
+    ) -> Any:
         """Find the common side length for an Archimedean vertex with the given face types."""
         multiplicities = Counter(faces_around_corner)
 
-        def length_to_angle_deficit(l):
+        def length_to_angle_deficit(l: float) -> float:
             return (
                 sum(
                     k * cls.regular_poly_in_angle(n, cls.platonic_side_length_to_radius(n, l))
@@ -171,7 +178,7 @@ def length_to_angle_deficit(l):
         )
 
     @classmethod
-    def archimedean_side_length_and_angles(cls, faces_around_corner):
+    def archimedean_side_length_and_angles(cls, faces_around_corner: Iterable[int]) -> tuple[float, dict[int, float]]:
         """Return the side length and a dict of interior angles for an Archimedean vertex."""
         length = cls.archimedean_side_length(faces_around_corner)
         return length, {
@@ -180,7 +187,7 @@ def archimedean_side_length_and_angles(cls, faces_around_corner):
         }
 
     @classmethod
-    def barycentric_to_euclidean_map(cls, tri):
+    def barycentric_to_euclidean_map(cls, tri: NDArray[Any]) -> Callable[[NDArray[Any]], Point]:
         """Return a callable mapping barycentric coordinates to points in the given triangle."""
         return lambda masses: cls.center_of_mass(tri, masses)
 
diff --git a/eucare/geometries/euclidean.py b/eucare/geometries/euclidean.py
index e0f03df..9759c9c 100755
--- a/eucare/geometries/euclidean.py
+++ b/eucare/geometries/euclidean.py
@@ -6,36 +6,39 @@
 
 from __future__ import annotations
 
+from typing import Any, Iterable
+
 import numpy as np
+from numpy.typing import NDArray
 
-from .base import Geometry
+from .base import Geometry, Transform
 
 
 class EuclideanGeometry(Geometry):
     """Flat 2D Euclidean geometry with standard vector operations."""
 
     @classmethod
-    def origin(cls):
+    def origin(cls) -> NDArray[Any]:
         return np.array([0, 0])
 
     @classmethod
-    def translation(cls, p1, p2):
-        def translate(p):
+    def translation(cls, p1: NDArray[Any], p2: NDArray[Any]) -> Transform:
+        def translate(p: NDArray[Any]) -> NDArray[Any]:
             return p + p2 - p1
 
         return translate
 
     @classmethod
-    def _rotate_around_origin(cls, a1):
+    def _rotate_around_origin(cls, a1: float) -> Transform:
         rot_mat = np.array([[np.cos(a1), np.sin(a1)], [-np.sin(a1), np.cos(a1)]])
 
-        def origin_rotate(p):
+        def origin_rotate(p: NDArray[Any]) -> NDArray[Any]:
             return p @ rot_mat
 
         return origin_rotate
 
     @classmethod
-    def center_of_mass(cls, points, masses=None):
+    def center_of_mass(cls, points: NDArray[Any], masses: NDArray[Any] | None = None) -> NDArray[Any]:
         assert len(points.shape) == 2 and points.shape[-1] == 2, f"{points.shape}"
         if masses is not None:
             masses = masses / np.sum(masses) * len(points)
@@ -43,23 +46,23 @@ def center_of_mass(cls, points, masses=None):
         return np.mean(points, axis=0)
 
     @classmethod
-    def distance_to_origin(cls, p):
-        return np.linalg.norm(p)
+    def distance_to_origin(cls, p: NDArray[Any]) -> float:
+        return float(np.linalg.norm(p))
 
     @classmethod
-    def angle_to_axis(cls, p):
+    def angle_to_axis(cls, p: NDArray[Any]) -> float:
         return np.arctan2(p[..., 1], p[..., 0])
 
     @classmethod
-    def point_along_axis(cls, x):
+    def point_along_axis(cls, x: float) -> NDArray[Any]:
         return np.array([x, 0])
 
     @classmethod
-    def to_euclidean(cls, pts):
+    def to_euclidean(cls, pts: NDArray[Any]) -> NDArray[Any]:
         return pts
 
     @classmethod
-    def archimedean_side_length(cls, faces_around_corner, eps=1e-6):
+    def archimedean_side_length(cls, faces_around_corner: Iterable[int], eps: float = 1e-6) -> float:
         euclidean_vertex_angle = sum(np.pi * (n - 2) / n for n in faces_around_corner)
         if abs(euclidean_vertex_angle - 2 * np.pi) < eps:
             return 1
diff --git a/eucare/geometries/hyperbolic.py b/eucare/geometries/hyperbolic.py
index e79bf6c..99622c7 100755
--- a/eucare/geometries/hyperbolic.py
+++ b/eucare/geometries/hyperbolic.py
@@ -8,12 +8,25 @@
 
 from __future__ import annotations
 
+from typing import Any, overload
+
 import numpy as np
+from numpy.typing import NDArray
 
 from .base import Geometry
 
+ComplexPoint = complex | NDArray[Any]
+
+
+@overload
+def apply_mobius(mat: NDArray[Any], points: complex) -> complex: ...
+
 
-def apply_mobius(mat, points):
+@overload
+def apply_mobius(mat: NDArray[Any], points: NDArray[Any]) -> NDArray[Any]: ...
+
+
+def apply_mobius(mat: NDArray[Any], points: ComplexPoint) -> ComplexPoint:
     """Apply a Mobius transformation given by a 2x2 matrix to complex-valued points."""
     return (mat[0, 0] * points + mat[0, 1]) / (mat[1, 0] * points + mat[1, 1])
 
@@ -21,53 +34,59 @@ def apply_mobius(mat, points):
 class MobiusTransform:
     """A Mobius transformation represented as a 2x2 complex matrix."""
 
-    def __init__(self, mat):
+    def __init__(self, mat: NDArray[Any] | list[list[complex]]) -> None:
         if not isinstance(mat, np.ndarray):
             mat = np.array(mat)
         assert mat.shape == (2, 2), f"{mat.shape}"
-        self.mat = mat
+        self.mat: NDArray[Any] = mat
+
+    @overload
+    def __call__(self, points: complex) -> complex: ...
+
+    @overload
+    def __call__(self, points: NDArray[Any]) -> NDArray[Any]: ...
 
-    def __call__(self, points):
+    def __call__(self, points: ComplexPoint) -> ComplexPoint:
         return apply_mobius(self.mat, points)
 
-    def __matmul__(self, other):
+    def __matmul__(self, other: MobiusTransform) -> MobiusTransform:
         assert isinstance(other, MobiusTransform), f"{type(other)}"
         return MobiusTransform(self.mat @ other.mat)
 
-    def __pow__(self, exponent):
+    def __pow__(self, exponent: int) -> MobiusTransform:
         return MobiusTransform(np.linalg.matrix_power(self.mat, exponent))
 
-    def __repr__(self):
+    def __repr__(self) -> str:
         return f"MobiusTransform({self.mat.tolist()})"
 
 
 # TODO: maybe have another class for the hyperboloid model
 
 
-def complex_to_real(z):
+def complex_to_real(z: NDArray[Any]) -> NDArray[Any]:
     """Convert complex numbers to real 2D coordinate arrays."""
     return np.stack([z.real, z.imag], axis=-1)
 
 
-def real_to_complex(x):
+def real_to_complex(x: NDArray[Any]) -> NDArray[Any]:
     """Convert real 2D coordinate arrays to complex numbers."""
     assert x.shape[-1] == 2
     return x[..., 0] + 1j * x[..., 1]
 
 
-def poincare_to_hyperboloid(z):
+def poincare_to_hyperboloid(z: NDArray[Any]) -> NDArray[Any]:
     """Map Poincare disk coordinates to the hyperboloid model."""
     pts = complex_to_real(z)
     squared_norm = (pts**2).sum(-1, keepdims=True)
     return np.concatenate([(1 + squared_norm), 2 * pts], axis=-1) / (1 - squared_norm)
 
 
-def hyperboloid_to_poincare(v):
+def hyperboloid_to_poincare(v: NDArray[Any]) -> NDArray[Any]:
     """Map hyperboloid model coordinates back to the Poincare disk."""
     return real_to_complex(v[..., 1:] / (1 + v[..., :1]))
 
 
-def hyperboloid_centroid(vs, ms=None, axis=None):
+def hyperboloid_centroid(vs: NDArray[Any], ms: NDArray[Any] | None = None, axis: int | None = None) -> NDArray[Any]:
     """Compute the centroid on the hyperboloid model, optionally weighted by masses."""
     if axis is None:
         assert len(vs.shape) == 2
@@ -78,7 +97,7 @@ def hyperboloid_centroid(vs, ms=None, axis=None):
     return mean
 
 
-def poincare_centroid(zs, ms=None, axis=None):
+def poincare_centroid(zs: NDArray[Any], ms: NDArray[Any] | None = None, axis: int | None = None) -> NDArray[Any]:
     """Compute the centroid of points in the Poincare disk via the hyperboloid model."""
     return hyperboloid_to_poincare(hyperboloid_centroid(poincare_to_hyperboloid(zs), ms, axis))
 
@@ -87,11 +106,11 @@ class PoincareDiskModel(Geometry):
     """Hyperbolic geometry using the Poincare disk model with complex coordinates."""
 
     @classmethod
-    def origin(cls):
+    def origin(cls) -> complex:
         return 0 + 0j
 
     @classmethod
-    def translation(cls, p1, p2):
+    def translation(cls, p1: complex, p2: complex) -> MobiusTransform:
         if p2 == 0:
             p1, p2 = 0, -p1
         if p1 == 0:
@@ -102,27 +121,27 @@ def translation(cls, p1, p2):
         return m3 @ m2 @ m1
 
     @classmethod
-    def rotation(cls, p1, a1):
+    def rotation(cls, p1: complex, a1: float) -> MobiusTransform:
         if p1 == 0:
             return MobiusTransform([[np.exp(1j * a1), 0], [0, 1]])
         return cls.translation(p1, 0) @ cls.rotation(0, a1) @ cls.translation(0, p1)
 
     @classmethod
-    def center_of_mass(cls, points, masses=None):
+    def center_of_mass(cls, points: NDArray[Any], masses: NDArray[Any] | None = None) -> NDArray[Any]:
         return poincare_centroid(points, masses)
 
     @classmethod
-    def distance_to_origin(cls, p):
-        return 2 * np.arctanh(np.linalg.norm(p))
+    def distance_to_origin(cls, p: ComplexPoint) -> float:
+        return float(2 * np.arctanh(np.linalg.norm(p)))
 
     @classmethod
-    def angle_to_axis(cls, p):
-        return np.arctan2(p.imag, p.real)
+    def angle_to_axis(cls, p: ComplexPoint) -> float:
+        return float(np.arctan2(p.imag, p.real))
 
     @classmethod
-    def point_along_axis(cls, x):
-        return np.sign(x) * np.tanh(np.abs(x) / 2)
+    def point_along_axis(cls, x: float) -> complex:
+        return complex(np.sign(x) * np.tanh(np.abs(x) / 2))
 
     @classmethod
-    def to_euclidean(cls, pts):
+    def to_euclidean(cls, pts: NDArray[Any]) -> NDArray[Any]:
         return complex_to_real(pts)
diff --git a/eucare/geometries/spherical.py b/eucare/geometries/spherical.py
index 92015fd..e95c1ca 100755
--- a/eucare/geometries/spherical.py
+++ b/eucare/geometries/spherical.py
@@ -7,16 +7,19 @@
 
 from __future__ import annotations
 
+from typing import Any, Iterable
+
 import numpy as np
+from numpy.typing import NDArray
 
-from .base import Geometry
+from .base import Geometry, Transform
 
 
-def _rot_x_mat(a1):
+def _rot_x_mat(a1: float) -> NDArray[Any]:
     return np.array([[1, 0, 0], [0, np.cos(a1), -np.sin(a1)], [0, np.sin(a1), np.cos(a1)]])
 
 
-def _rot_z_mat(a1):
+def _rot_z_mat(a1: float) -> NDArray[Any]:
     return np.array(
         [
             [
@@ -34,19 +37,19 @@ class SphereModel(Geometry):
     """Spherical geometry with points on the unit 2-sphere in R^3."""
 
     @classmethod
-    def origin(cls):
+    def origin(cls) -> NDArray[Any]:
         return np.array([1, 0, 0])
 
     @classmethod
-    def translation(cls, p1, p2):
-        def origin_translation_mat(p1):
+    def translation(cls, p1: NDArray[Any], p2: NDArray[Any]) -> Transform:
+        def origin_translation_mat(p1: NDArray[Any]) -> NDArray[Any]:
             a1 = np.arctan2(p1[2], p1[1])
             m1 = _rot_x_mat(-a1)
             m2 = _rot_z_mat(-np.arccos(p1[0]))
             m3 = _rot_x_mat(a1)
             return m3 @ m2 @ m1
 
-        def minus(p):
+        def minus(p: NDArray[Any]) -> NDArray[Any]:
             return np.array([p[0], *-p[1:]])
 
         m1 = origin_translation_mat(minus(p1))
@@ -55,22 +58,22 @@ def minus(p):
 
         mat = m3 @ m2 @ m1
 
-        def translate(p):
+        def translate(p: NDArray[Any]) -> NDArray[Any]:
             return p @ mat
 
         return translate
 
     @classmethod
-    def _rotate_around_origin(cls, a1):
+    def _rotate_around_origin(cls, a1: float) -> Transform:
         mat = _rot_x_mat(a1).T
 
-        def origin_rotate(p):
+        def origin_rotate(p: NDArray[Any]) -> NDArray[Any]:
             return p @ mat
 
         return origin_rotate
 
     @classmethod
-    def center_of_mass(cls, points, masses=None):
+    def center_of_mass(cls, points: NDArray[Any], masses: NDArray[Any] | None = None) -> NDArray[Any]:
         if masses is not None:
             masses = masses / np.sum(masses) * len(points)
             points = points * masses[..., None]
@@ -79,28 +82,28 @@ def center_of_mass(cls, points, masses=None):
         return result
 
     @classmethod
-    def distance_to_origin(cls, p):
-        return np.arccos(np.clip(p[0], -1, 1))
+    def distance_to_origin(cls, p: NDArray[Any]) -> float:
+        return float(np.arccos(np.clip(p[0], -1, 1)))
 
     @classmethod
-    def angle_to_axis(cls, p):
-        return np.arctan2(p[2], p[1])
+    def angle_to_axis(cls, p: NDArray[Any]) -> float:
+        return float(np.arctan2(p[2], p[1]))
 
     @classmethod
-    def point_along_axis(cls, x):
+    def point_along_axis(cls, x: float) -> NDArray[Any]:
         return np.array([np.cos(x), np.sin(x), 0])
 
     @classmethod
-    def archimedean_side_length(cls, faces_around_corner):
+    def archimedean_side_length(cls, faces_around_corner: Iterable[int], **kwargs: Any) -> float:
         return super().archimedean_side_length(faces_around_corner, bracket=[0.1, 2 * np.pi / max(faces_around_corner)])
 
     # --- Methods Specific to this geometry ---
 
     @classmethod
-    def stereographic_projection(cls, pts):
+    def stereographic_projection(cls, pts: NDArray[Any]) -> NDArray[Any]:
         """Apply stereographic projection with pole at (-1, 0, 0)."""
         return 2 * pts[..., 1:] / (pts[..., :1] + 1)
 
     @classmethod
-    def to_euclidean(cls, pts):
+    def to_euclidean(cls, pts: NDArray[Any]) -> NDArray[Any]:
         return cls.stereographic_projection(pts)
diff --git a/eucare/half.py b/eucare/half.py
index ce18e2a..5523e64 100755
--- a/eucare/half.py
+++ b/eucare/half.py
@@ -20,14 +20,16 @@
 from copy import copy, deepcopy
 from itertools import chain
 from math import pi
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, Any, Callable, Iterable, Iterator, Literal, TypeVar, overload
 
 if TYPE_CHECKING:
+    from .geometries.base import Geometry
     from .rendering import Rendering
 
 import matplotlib.pyplot as plt
 import networkx as nx
 import numpy as np
+from numpy.typing import NDArray
 
 from .base import edge_lengths_and_in_angles, signed_area
 from .geometries import EuclideanGeometry
@@ -40,40 +42,40 @@ class AttributeObject:
 
     def __init__(self) -> None:
         super(AttributeObject, self).__init__()
-        self.attributes: dict = dict()
+        self.attributes: dict[str, Any] = dict()
 
     def has_attributes(self) -> bool:
         """Return True if this object has any attributes set."""
         return bool(self.attributes)
 
-    def __getitem__(self, attr):
+    def __getitem__(self, attr: str) -> Any:
         return self.attributes[attr]
 
-    def __setitem__(self, attr, value) -> None:
+    def __setitem__(self, attr: str, value: Any) -> None:
         self.attributes[attr] = value
 
-    def __iter__(self):
+    def __iter__(self) -> Iterator[str]:
         return iter(self.attributes)
 
-    def __delitem__(self, key) -> None:
+    def __delitem__(self, key: str) -> None:
         del self.attributes[key]
 
-    def __contains__(self, item) -> bool:
+    def __contains__(self, item: object) -> bool:
         return item in self.attributes
 
-    def get(self, attr, *args, **kwargs):
+    def get(self, attr: str, *args: Any, **kwargs: Any) -> Any:
         """Return ``self.attributes.get(attr, *args, **kwargs)``."""
         return self.attributes.get(attr, *args, **kwargs)
 
-    def keys(self):
+    def keys(self) -> Iterable[str]:
         """Return a view of the attribute keys."""
         return self.attributes.keys()
 
-    def values(self):
+    def values(self) -> Iterable[Any]:
         """Return a view of the attribute values."""
         return self.attributes.values()
 
-    def items(self):
+    def items(self) -> Iterable[tuple[str, Any]]:
         """Return a view of the (key, value) attribute items."""
         return self.attributes.items()
 
@@ -85,7 +87,7 @@ class IdObject(AttributeObject):
     independent test runs).
     """
 
-    current_ids: dict = dict()
+    current_ids: dict[type, int] = dict()
 
     def __init__(self) -> None:
         super(IdObject, self).__init__()
@@ -107,9 +109,9 @@ def reset_ids(cls) -> None:
             IdObject.current_ids[cls] = 0
 
 
-def check_cyclic_iterator_consistency(iterator) -> None:
+def check_cyclic_iterator_consistency(iterator: Iterable[Any]) -> None:
     """Assert that every item yielded by ``iterator`` is unique (no cycles repeat)."""
-    items = set()
+    items: set[Any] = set()
     for item in iterator:
         assert item not in items, f"{iterator}, {item}, {items}"
         items.add(item)
@@ -128,10 +130,25 @@ class HalfEdge(IdObject):
         orig: The origin vertex.
         dest: The destination vertex.
         face: The face to the left, or ``None`` for border edges.
+
+    Note:
+        ``rev``/``nex``/``pre``/``orig``/``dest`` are typed as non-optional
+        because the data structure invariant requires them to be set once the
+        graph is fully constructed.  The constructor accepts ``None`` to
+        support incremental graph building; in that case the corresponding
+        attribute is simply left unset and reading it before it is assigned
+        raises :class:`AttributeError`.
     """
 
     printname = "HE"
 
+    rev: HalfEdge
+    nex: HalfEdge
+    pre: HalfEdge
+    orig: Vertex
+    dest: Vertex
+    face: Face | None
+
     def __init__(
         self,
         rev: HalfEdge | None = None,
@@ -142,22 +159,23 @@ def __init__(
         face: Face | None = None,
     ) -> None:
         super(HalfEdge, self).__init__()
-        # HalfEdge references
-        self.rev = rev
-        self.nex = nex
-        self.pre = pre
-
-        # Vertex references
-        self.orig = orig
-        self.dest = dest
-
-        # Face reference
+        if rev is not None:
+            self.rev = rev
+        if nex is not None:
+            self.nex = nex
+        if pre is not None:
+            self.pre = pre
+        if orig is not None:
+            self.orig = orig
+        if dest is not None:
+            self.dest = dest
+        # Face is genuinely optional: None on border edges.
         self.face = face
 
     # def __repr__(self):
     #    return f'{self}({self.orig},{self.dest})'
 
-    def __str__(self):
+    def __str__(self) -> str:
         return f'H{self["id"]}'
 
     def on_border(self) -> bool:
@@ -168,7 +186,7 @@ def check_consistency(self) -> None:
         """Assert that ``self.rev.rev is self``."""
         assert self.rev.rev is self, f"{self}, {self.rev}, {self.rev.rev}"
 
-    def midpoint(self) -> np.ndarray:
+    def midpoint(self) -> NDArray:
         """Return the midpoint of the underlying edge in Euclidean position space."""
         return np.mean([v["pos"] for v in (self.orig, self.dest)], axis=0)
 
@@ -179,15 +197,24 @@ class Vertex(IdObject):
     Links to one arbitrary outgoing half-edge (``any_outgoing``).  All incident
     half-edges and faces are reachable via ``outgoing_iter()`` and
     ``face_iter()``.
+
+    Note:
+        ``any_outgoing`` is typed as non-optional because every fully
+        constructed vertex points to some outgoing half-edge.  The constructor
+        accepts ``None`` to support incremental graph building; in that case
+        the attribute is left unset until assigned.
     """
 
     printname = "V"
 
+    any_outgoing: HalfEdge
+
     def __init__(self, any_outgoing: HalfEdge | None = None) -> None:
         super(Vertex, self).__init__()
-        self.any_outgoing = any_outgoing
+        if any_outgoing is not None:
+            self.any_outgoing = any_outgoing
 
-    def outgoing_iter(self):
+    def outgoing_iter(self) -> Iterator[HalfEdge]:
         """Yield outgoing half-edges in counter-clockwise order around this vertex."""
         initial = self.any_outgoing
         current = initial
@@ -197,7 +224,7 @@ def outgoing_iter(self):
             if current is initial:
                 break
 
-    def reverse_outgoing_iter(self):
+    def reverse_outgoing_iter(self) -> Iterator[HalfEdge]:
         """Yield outgoing half-edges in clockwise order around this vertex."""
         initial = self.any_outgoing
         current = initial
@@ -211,28 +238,28 @@ def order(self) -> int:
         """Return the number of edges incident to this vertex (its degree)."""
         return len(list(self.outgoing_iter()))
 
-    def incoming_iter(self):
+    def incoming_iter(self) -> Iterator[HalfEdge]:
         """Yield incoming half-edges in counter-clockwise order around this vertex."""
         for h in self.outgoing_iter():
             yield h.rev
 
-    def face_iter(self):
+    def face_iter(self) -> Iterator[Face | None]:
         """Yield adjacent faces (including ``None`` for border) in CCW order."""
         for h in self.outgoing_iter():
             yield h.face
 
-    def true_face_iter(self):
+    def true_face_iter(self) -> Iterator[Face]:
         """Yield only the non-``None`` adjacent faces."""
         for f in self.face_iter():
             if f is not None:
                 yield f
 
-    def vertex_iter(self):
+    def vertex_iter(self) -> Iterator[Vertex]:
         """Yield neighbouring vertices in counter-clockwise order around this vertex."""
         for h in self.outgoing_iter():
             yield h.dest
 
-    def common_faces_iter(self, other: "Vertex"):
+    def common_faces_iter(self, other: Vertex) -> Iterator[Face]:
         """Yield faces incident to both this vertex and ``other``."""
         fs = set(other.face_iter())
         for f in self.face_iter():
@@ -281,13 +308,22 @@ class Face(IdObject):
     Links to one arbitrary bounding half-edge (``any_side``).  Iterate
     over boundary half-edges with ``halfedge_iter()`` and vertices with
     ``vertex_iter()``.
+
+    Note:
+        ``any_side`` is typed as non-optional because every fully constructed
+        face references one of its boundary half-edges.  The constructor
+        accepts ``None`` to support incremental graph building; in that case
+        the attribute is left unset until assigned.
     """
 
+    any_side: HalfEdge
+
     def __init__(self, any_side: HalfEdge | None = None) -> None:
         super(Face, self).__init__()
-        self.any_side = any_side
+        if any_side is not None:
+            self.any_side = any_side
 
-    def halfedge_iter(self):
+    def halfedge_iter(self) -> Iterator[HalfEdge]:
         """Yield boundary half-edges in counter-clockwise order around this face."""
         initial = self.any_side
         current = initial
@@ -297,7 +333,7 @@ def halfedge_iter(self):
             if current is initial:
                 break
 
-    def reverse_halfedge_iter(self):
+    def reverse_halfedge_iter(self) -> Iterator[HalfEdge]:
         """Yield the reverse of each boundary half-edge (i.e. those facing outward)."""
         for h in self.halfedge_iter():
             yield h.rev
@@ -306,17 +342,17 @@ def order(self) -> int:
         """Return the number of edges (= sides) of this polygonal face."""
         return len(list(self.halfedge_iter()))
 
-    def vertex_iter(self):
+    def vertex_iter(self) -> Iterator[Vertex]:
         """Yield boundary vertices in counter-clockwise order around this face."""
         for h in self.halfedge_iter():
             yield h.orig
 
-    def face_iter(self):
+    def face_iter(self) -> Iterator[Face | None]:
         """Yield neighbouring faces (including ``None`` for border) in CCW order."""
         for h in self.halfedge_iter():
             yield h.rev.face
 
-    def true_face_iter(self):
+    def true_face_iter(self) -> Iterator[Face]:
         """Yield only the non-``None`` neighbouring faces."""
         for f in self.face_iter():
             if f is not None:
@@ -326,15 +362,15 @@ def on_border(self) -> bool:
         """Return True if any boundary half-edge of this face is adjacent to a border edge."""
         return any(h.rev.on_border() for h in self.halfedge_iter())
 
-    def outgoing_edge_at(self, v: "Vertex") -> HalfEdge:
+    def outgoing_edge_at(self, v: Vertex) -> HalfEdge:
         """Return the boundary half-edge of this face originating at vertex ``v``."""
         return next(h for h in self.halfedge_iter() if h.orig is v)
 
-    def incoming_edge_at(self, v: "Vertex") -> HalfEdge:
+    def incoming_edge_at(self, v: Vertex) -> HalfEdge:
         """Return the boundary half-edge of this face ending at vertex ``v``."""
         return next(h for h in self.halfedge_iter() if h.dest is v)
 
-    def common_halfedge_iter(self, other: "Face"):
+    def common_halfedge_iter(self, other: Face) -> Iterator[HalfEdge]:
         """Yield half-edges of this face whose reverse lies on ``other``."""
         assert isinstance(other, Face)
         hs = set(h.rev for h in other.halfedge_iter())
@@ -342,29 +378,32 @@ def common_halfedge_iter(self, other: "Face"):
             if h in hs:
                 yield h
 
-    def common_vertex_iter(self, other: "Face"):
+    def common_vertex_iter(self, other: Face) -> Iterator[Vertex]:
         """Yield vertices shared with ``other``."""
         vs = set(other.vertex_iter())
         for v in self.vertex_iter():
             if v in vs:
                 yield v
 
-    def common_face_iter(self, other: "Face"):
+    def common_face_iter(self, other: Face) -> Iterator[Face | None]:
         """Yield faces neighbouring both this face and ``other``."""
         fs = set(other.face_iter())
         for f in self.face_iter():
             if f in fs:
                 yield f
 
-    def midpoint(self) -> np.ndarray:
+    def midpoint(self) -> NDArray:
         """Return the (cached) face midpoint, falling back to the vertex centroid."""
-        return self.attributes.get("midpoint", np.mean([v["pos"] for v in self.vertex_iter()], axis=0))
+        cached = self.attributes.get("midpoint")
+        if cached is not None:
+            return cached
+        return np.mean([v["pos"] for v in self.vertex_iter()], axis=0)
 
-    def pseudo_incenter(self) -> np.ndarray:
+    def pseudo_incenter(self) -> NDArray:
         """Return the pseudo-incenter (true incenter for tangential polygons)."""
         return pseudo_incenter(self)
 
-    def recompute_lengths_and_angles(self, geometry) -> None:
+    def recompute_lengths_and_angles(self, geometry: type[Geometry]) -> None:
         """Recompute the ``length`` and ``in_angle`` attributes from current vertex positions."""
         points = np.stack([v["pos"] for v in self.vertex_iter()])
         lengths, angles = edge_lengths_and_in_angles(points, geometry)
@@ -401,7 +440,7 @@ def pseudo_incenter(f: Face | np.ndarray) -> np.ndarray:
     return incenter
 
 
-def pseudo_circumcenter(ps, return_radius=False) -> np.ndarray | tuple[np.ndarray, float]:
+def pseudo_circumcenter(ps: NDArray, return_radius: bool = False) -> NDArray | tuple[NDArray, float]:
     """Compute the "best fit" circumcenter of a polygon.
 
     Finds the interior point that minimizes the standard deviation of
@@ -444,6 +483,9 @@ def pseudo_circumcenter(ps, return_radius=False) -> np.ndarray | tuple[np.ndarra
     return center
 
 
+G = TypeVar("G", bound="HalfEdgeGraph")
+
+
 class HalfEdgeGraph:
     """Topology-only half-edge graph (DCEL).
 
@@ -453,7 +495,13 @@ class HalfEdgeGraph:
     ``.copy()`` first if you need the original.
     """
 
-    def __init__(self, other: "HalfEdgeGraph | None" = None) -> None:
+    halfedges: set[HalfEdge]
+    vertices: set[Vertex]
+    faces: set[Face]
+    _any_border: HalfEdge | None
+    simply_connected: bool
+
+    def __init__(self, other: HalfEdgeGraph | None = None) -> None:
         if other is not None:
             self.halfedges = copy(other.halfedges)
             self.vertices = copy(other.vertices)
@@ -488,7 +536,7 @@ def add_vertex(self, v: Vertex) -> None:
         """Register vertex ``v`` with this graph."""
         self.vertices.add(v)
 
-    def add_vertices(self, vs) -> None:
+    def add_vertices(self, vs: Iterable[Vertex]) -> None:
         """Register an iterable of vertices with this graph."""
         self.vertices.update(vs)
 
@@ -496,7 +544,7 @@ def add_face(self, f: Face) -> None:
         """Register face ``f`` with this graph."""
         self.faces.add(f)
 
-    def add_faces(self, fs) -> None:
+    def add_faces(self, fs: Iterable[Face]) -> None:
         """Register an iterable of faces with this graph."""
         self.faces.update(fs)
 
@@ -504,7 +552,7 @@ def add_halfedge(self, h: HalfEdge) -> None:
         """Register half-edge ``h`` with this graph."""
         self.halfedges.add(h)
 
-    def add_halfedges(self, hs) -> None:
+    def add_halfedges(self, hs: Iterable[HalfEdge]) -> None:
         """Register an iterable of half-edges with this graph."""
         self.halfedges.update(hs)
 
@@ -512,11 +560,11 @@ def delete_face(self, f: Face) -> None:
         """Remove a single face (and any newly-orphaned edges/vertices)."""
         self.delete_faces({f})
 
-    def fill_holes(self):
+    def fill_holes(self) -> None:
         """Fill all 'holes' inside the graph: Add faces to enclosed areas which are not yet a face"""
         raise NotImplementedError
 
-    def delete_faces(self, fs) -> None:
+    def delete_faces(self, fs: Iterable[Face]) -> None:
         """Delete all faces in ``fs`` and any half-edges/vertices they leave dangling."""
         self.faces.difference_update(set(fs))
         # 1. Determine edges that have to be deleted, remove them and their rev's from Graph
@@ -558,7 +606,7 @@ def delete_faces(self, fs) -> None:
                 if v.any_outgoing not in self.halfedges:  # no outgoing edge is still in graph
                     self.vertices.remove(v)
 
-    def delete_subset(self, *items) -> None:
+    def delete_subset(self, *items: Face | HalfEdge | Vertex | Iterable[Face | HalfEdge | Vertex]) -> None:
         """Delete a mixed collection of faces, half-edges, and vertices, repairing topology.
 
         Accepts individual ``Face``/``HalfEdge``/``Vertex`` objects or iterables
@@ -670,6 +718,7 @@ def delete_edge(self, h: HalfEdge) -> None:
         """
         if h.on_border():
             raise NotImplementedError
+        assert h.face is not None and h.rev.face is not None
         if h.face is h.rev.face:
             assert (
                 h.orig.order() == 1 or h.dest.order() == 1
@@ -711,6 +760,7 @@ def join_vertex(self, v: Vertex) -> None:
             h.pre = h.pre.pre
             h.pre.nex = h
             if not h.on_border():
+                assert h.face is not None
                 h.face.any_side = h
         self.vertices.remove(v)
 
@@ -722,6 +772,12 @@ def join_order_2_boundary_vertices(self) -> None:
                 to_join.append(v)
         for v in to_join:
             self.join_vertex(v)
+        # InAngleHEG and subclasses have ``recompute_lengths_and_angles``;
+        # plain HalfEdgeGraph does not, so this method only makes sense when
+        # called on a subclass that recomputes angles.
+        assert hasattr(
+            self, "recompute_lengths_and_angles"
+        ), "join_order_2_boundary_vertices requires an InAngleHEG subclass"
         self.recompute_lengths_and_angles()
 
     def join_edge(self, h: HalfEdge) -> Vertex:
@@ -733,8 +789,8 @@ def join_edge(self, h: HalfEdge) -> Vertex:
         Returns:
             The remaining vertex (originally ``h.orig``); ``h.dest`` is removed.
         """
-        assert h.on_border() or h.face.order() > 1
-        assert h.rev.on_border() or h.rev.face.order() > 1
+        assert h.on_border() or (h.face is not None and h.face.order() > 1)
+        assert h.rev.on_border() or (h.rev.face is not None and h.rev.face.order() > 1)
 
         v1, v2 = h.orig, h.dest
 
@@ -751,6 +807,7 @@ def join_edge(self, h: HalfEdge) -> Vertex:
             h2.pre.nex = h2.nex
             h2.nex.pre = h2.pre
             if not h2.on_border():
+                assert h2.face is not None
                 h2.face.any_side = h2.nex
             h2.dest.any_outgoing = h2.nex
 
@@ -803,7 +860,9 @@ def subdivide_edge(
 
         return h2, v
 
-    def subdivide_face(self, f: Face, v1: Vertex, v2: Vertex, **halfedge_attributes: object) -> tuple[HalfEdge, Face]:
+    def subdivide_face(
+        self, f: Face | None, v1: Vertex, v2: Vertex, **halfedge_attributes: object
+    ) -> tuple[HalfEdge, Face]:
         """Subdivide face ``f`` along a new edge from ``v1`` to ``v2``.
 
         Two new half-edges ``h12`` (``v1 -> v2``) and ``h21`` (``v2 -> v1``)
@@ -853,7 +912,7 @@ def halfedges_representing_edges(self) -> set[HalfEdge]:
             A set containing exactly one of ``{h, h.rev}`` for every
             undirected edge of the graph.
         """
-        result = set()
+        result: set[HalfEdge] = set()
         for h in self.halfedges:
             if h.rev not in result:
                 result.add(h)
@@ -861,7 +920,7 @@ def halfedges_representing_edges(self) -> set[HalfEdge]:
 
     def to_networkx_undirected(self) -> nx.Graph:
         """Return a :class:`networkx.Graph` of this graph's underlying undirected topology."""
-        result = nx.Graph()
+        result: nx.Graph = nx.Graph()
         result.add_edges_from([(h.orig, h.dest) for h in self.halfedges])
         return result
 
@@ -875,7 +934,7 @@ def get_any_border(self) -> HalfEdge:
                 return h
         raise LookupError("No border found.")
 
-    def border_edge_iter(self):
+    def border_edge_iter(self) -> Iterator[HalfEdge]:
         """Yield border half-edges. For simply connected graphs, in cyclic order."""
         if self.simply_connected:
             initial = self.get_any_border()
@@ -894,7 +953,7 @@ def border_edges(self) -> list[HalfEdge]:
         """Return all border half-edges as a list."""
         return list(self.border_edge_iter())
 
-    def border_vertex_iter(self):
+    def border_vertex_iter(self) -> Iterator[Vertex]:
         """Yield vertices lying on a border edge."""
         for h in self.border_edge_iter():
             yield h.orig
@@ -967,14 +1026,16 @@ def glue_e2e(self, e1: HalfEdge, e2: HalfEdge) -> None:
             # handle face stuff now
             if e1.face is e2.face is None:
                 pass
-            elif e1.nex is e2 and e2.nex is e1:
-                self.faces.remove(e1.face)
-            elif e1.nex is e2:
-                e1.face.any_side = e2.nex
-            elif e2.nex is e1:
-                e1.face.any_side = e1.nex
             else:
-                raise ValueError(f"Cannot glue: Edges {[e1, e2]} of face {e1.face} are not adjacent.")
+                assert e1.face is not None
+                if e1.nex is e2 and e2.nex is e1:
+                    self.faces.remove(e1.face)
+                elif e1.nex is e2:
+                    e1.face.any_side = e2.nex
+                elif e2.nex is e1:
+                    e1.face.any_side = e1.nex
+                else:
+                    raise ValueError(f"Cannot glue: Edges {[e1, e2]} of face {e1.face} are not adjacent.")
 
         # glue vertices
         for v1_out, v2_out in ((e1, e2.nex), (e1.nex, e2)):
@@ -994,12 +1055,21 @@ def glue_graph_e2e(self, graph: "HalfEdgeGraph", e1: HalfEdge, e2: HalfEdge) ->
         self.add_graph(graph)
         self.glue_e2e(e1, e2)
 
-    def close_vertex(self, v: Vertex) -> None:
-        """Close a boundary corner at ``v`` by gluing its two incident border edges."""
+    def close_vertex(self, v: Vertex, reverse: bool = False) -> Vertex:
+        """Close a boundary corner at ``v`` by gluing its two incident border edges.
+
+        Args:
+            v: Boundary vertex whose two adjacent border edges are to be sewn.
+            reverse: Swaps which side's vertex is kept (overridden by subclasses).
+
+        Returns:
+            The remaining vertex after gluing.
+        """
         # get edges to be glued
         e1 = v.get_outgoing_border()
         e2 = e1.pre
         self.glue_e2e(e1, e2)
+        return e1.rev.orig
 
     def twocolorable(self) -> bool:
         """Return True if every interior vertex has even order (necessary for face 2-coloring)."""
@@ -1019,7 +1089,7 @@ def twocolor_faces(self, key: str = "color_key", initial_face: Face | None = Non
         if initial_face is None:
             initial_face = next(iter(self.faces))
         yet_to_color = copy(self.faces)
-        frontier = {(initial_face, False)}
+        frontier: set[tuple[Face, bool]] = {(initial_face, False)}
         while frontier:
             face, label = frontier.pop()
             if face not in yet_to_color:
@@ -1027,10 +1097,16 @@ def twocolor_faces(self, key: str = "color_key", initial_face: Face | None = Non
             yet_to_color.remove(face)
             face[key] = label
             for f in face.face_iter():
-                frontier.add((f, not label))
+                if f is not None:
+                    frontier.add((f, not label))
         assert not yet_to_color, "Graph is not connected!"
 
-    def execute_edge_instruction(self, h: HalfEdge, instruction=None, key: str | None = None) -> None:
+    def execute_edge_instruction(
+        self,
+        h: HalfEdge,
+        instruction: Callable[[HalfEdgeGraph, HalfEdge], None] | None = None,
+        key: str | None = None,
+    ) -> None:
         """Run a tile-gluing or growth ``instruction`` on the border half-edge ``h``.
 
         If ``instruction`` is omitted, it is read from ``h[key]`` (default key
@@ -1043,12 +1119,20 @@ def execute_edge_instruction(self, h: HalfEdge, instruction=None, key: str | Non
             assert key is None, "Please specify not more than one of [key, instruction]."
         instruction(self, h)
 
-    def execute_all_edge_instructions(self, instruction=None, key: str | None = None) -> None:
+    def execute_all_edge_instructions(
+        self,
+        instruction: Callable[[HalfEdgeGraph, HalfEdge], None] | None = None,
+        key: str | None = None,
+    ) -> None:
         """Execute :meth:`execute_edge_instruction` on every border half-edge."""
         for h in self.border_edges():
             self.execute_edge_instruction(h, instruction, key)
 
-    def show_spring_layout(self, figsize: tuple[float, float] = (15, 15), emph_func: "callable | None" = None) -> None:
+    def show_spring_layout(
+        self,
+        figsize: tuple[float, float] = (15, 15),
+        emph_func: Callable[[HalfEdge], bool] | None = None,
+    ) -> None:
         """Display the underlying undirected graph using NetworkX's spring layout.
 
         Args:
@@ -1057,12 +1141,12 @@ def show_spring_layout(self, figsize: tuple[float, float] = (15, 15), emph_func:
                 Defaults to highlighting edges flagged with the ``'delete'`` attribute.
         """
 
-        G = nx.Graph()
+        G: nx.Graph = nx.Graph()
         G.add_edges_from([(h.orig, h.dest) for h in self.halfedges])
 
         if emph_func is None:
 
-            def emph_func(h):
+            def emph_func(h: HalfEdge) -> bool:
                 return h.attributes.get("delete", False)
 
         G.add_edges_from([(h.orig, h.dest) for h in self.halfedges if emph_func(h)], color="r")
@@ -1133,7 +1217,7 @@ def check_consistency(self) -> None:
                 "faces: %s, %s", referenced_faces.difference(self.faces), self.faces.difference(referenced_faces)
             )
 
-            reference_dict = {
+            reference_dict: dict[object, set[tuple[object, str]]] = {
                 obj: set()
                 for obj in referenced_halfedges.union(referenced_vertices).union(referenced_faces).union([None])
             }
@@ -1154,9 +1238,24 @@ def check_consistency(self) -> None:
                 logger.error("%s referenced by %s.", obj, reference_dict[obj])
             raise RuntimeError("Graph consistency check failed. See log for details.")
 
+    @overload
+    def copy(
+        self: G,
+        deepcopy_attributes: bool = ...,
+        *,
+        return_mappings: Literal[True],
+    ) -> tuple[G, tuple[dict[Vertex, Vertex], dict[HalfEdge, HalfEdge], dict[Face, Face]]]: ...
+
+    @overload
+    def copy(
+        self: G,
+        deepcopy_attributes: bool = ...,
+        return_mappings: Literal[False] = ...,
+    ) -> G: ...
+
     def copy(
         self, deepcopy_attributes: bool = False, return_mappings: bool = False
-    ) -> "HalfEdgeGraph | tuple[HalfEdgeGraph, tuple[dict, dict, dict]]":
+    ) -> HalfEdgeGraph | tuple[HalfEdgeGraph, tuple[dict[Vertex, Vertex], dict[HalfEdge, HalfEdge], dict[Face, Face]]]:
         """Return an independent copy of this graph.
 
         Args:
@@ -1169,24 +1268,24 @@ def copy(
             The copied graph -- a tuple ``(graph, (v_map, e_map, f_map))`` when ``return_mappings`` is True.
         """
 
-        def copy_with_attributes(obj):
+        def copy_with_attributes(obj: Any) -> Any:
             cls = type(obj)
-            new = cls.__new__(cls)
+            new = cls.__new__(cls)  # type: ignore[call-overload]
             new.attributes = copy(obj.attributes)
             return new
 
-        def copy_with_attributes_deep(obj):
+        def copy_with_attributes_deep(obj: Any) -> Any:
             cls = type(obj)
-            new = cls.__new__(cls)
+            new = cls.__new__(cls)  # type: ignore[call-overload]
             new.attributes = deepcopy(obj.attributes)
             return new
 
         copy_func = copy_with_attributes_deep if deepcopy_attributes else copy_with_attributes
 
         # init mappings from old to new vertex/halfedge/face objects
-        v_map, e_map, f_map = [
-            {obj: copy_func(obj) for obj in container} for container in (self.vertices, self.halfedges, self.faces)
-        ]
+        v_map: dict[Vertex, Vertex] = {v: copy_func(v) for v in self.vertices}
+        e_map: dict[HalfEdge, HalfEdge] = {e: copy_func(e) for e in self.halfedges}
+        f_map: dict[Face, Face] = {f: copy_func(f) for f in self.faces}
 
         # copy other potential attributes of graph (e.g. tau for InAngleHEG)
         cls = type(self)
@@ -1203,14 +1302,16 @@ def copy_with_attributes_deep(obj):
         for f, f_new in f_map.items():
             f_new.any_side = e_map[f.any_side]
 
-        f_map[None] = None  # to handle border
+        # face map extended with None for border edges
+        face_map_with_border: dict[Face | None, Face | None] = dict(f_map.items())
+        face_map_with_border[None] = None
         for e, e_new in e_map.items():
             e_new.orig = v_map[e.orig]
             e_new.dest = v_map[e.dest]
             e_new.nex = e_map[e.nex]
             e_new.pre = e_map[e.pre]
             e_new.rev = e_map[e.rev]
-            e_new.face = f_map[e.face]
+            e_new.face = face_map_with_border[e.face]
 
         if not return_mappings:
             return result
@@ -1349,21 +1450,31 @@ class GeometricHEG(InAngleHEG):
     used for position recomputation and length/angle calculations.
     """
 
-    def __init__(self, geometry: object = EuclideanGeometry, **super_kwargs: object) -> None:
+    geometry: type[Geometry]
+
+    def __init__(
+        self,
+        geometry: type[Geometry] = EuclideanGeometry,
+        angle_sum: float | None = None,
+        eps: float | None = None,
+        other: HalfEdgeGraph | None = None,
+    ) -> None:
         """Create a graph backed by the given geometry.
 
         Args:
             geometry: A geometry backend exposing ``distance``, ``angle``,
                 ``construct_next_poly_point``, and ``to_euclidean``. Defaults
                 to :class:`EuclideanGeometry`.
-            **super_kwargs: Forwarded to :class:`InAngleHEG`.
+            angle_sum: Forwarded to :class:`InAngleHEG`.
+            eps: Forwarded to :class:`InAngleHEG`.
+            other: Forwarded to :class:`InAngleHEG`.
         """
-        super(GeometricHEG, self).__init__(**super_kwargs)
+        super().__init__(angle_sum=angle_sum, eps=eps, other=other)
         self.geometry = geometry
 
-    def positions_coincide(self, p1: np.ndarray, p2: np.ndarray) -> bool:
+    def positions_coincide(self, p1: NDArray, p2: NDArray) -> bool:
         """Return True if positions ``p1`` and ``p2`` are within ``eps``."""
-        return np.linalg.norm(p1 - p2) < self.eps
+        return bool(np.linalg.norm(p1 - p2) < self.eps)
 
     def lengths_equal(self, l1: float, l2: float) -> bool:
         """Return True if lengths ``l1`` and ``l2`` differ by at most ``eps``."""
@@ -1383,7 +1494,7 @@ def join_edge(self, h: HalfEdge) -> Vertex:
         v["pos"] = new_pos
         return v
 
-    def recompute_positions(self, edge_to_start: HalfEdge | None = None, faces: "set[Face] | None" = None) -> None:
+    def recompute_positions(self, edge_to_start: HalfEdge | None = None, faces: Iterable[Face] | None = None) -> None:
         """Recompute vertex positions from edge lengths and interior angles.
 
         Performs a heap-ordered BFS over ``faces``, propagating positions
@@ -1397,31 +1508,34 @@ def recompute_positions(self, edge_to_start: HalfEdge | None = None, faces: "set
                 computation. Defaults to the longest interior edge.
             faces: Iterable of faces to process. Defaults to all faces.
         """
-        if faces is None:
-            faces = self.faces
-        if len(faces) == 0:
+        face_set: set[Face] = self.faces if faces is None else set(faces)
+        if len(face_set) == 0:
             return  # nothing to do
         if edge_to_start is None:
             # select longest edge
             hs = [h for h in self.halfedges if not h.on_border()]
             lengths = [h["length"] for h in hs]
-            edge_to_start = hs[np.argmax(lengths)]
+            edge_to_start = hs[int(np.argmax(lengths))]
         if edge_to_start.on_border():
             raise ValueError(f"edge_to_start must not be on border, got {edge_to_start}")
+        assert edge_to_start.face is not None
 
         # delete old positions
-        vertices = set.union(set(), *(f.vertex_iter() for f in faces)).difference(
-            {edge_to_start.orig, edge_to_start.dest}
-        )
+        face_vertices: set[Vertex] = set()
+        for f in face_set:
+            face_vertices.update(f.vertex_iter())
+        vertices = face_vertices.difference({edge_to_start.orig, edge_to_start.dest})
         for v in vertices:
             if "pos" in v.attributes:
                 del v["pos"]
 
         # compute positions for new vertices, face by face
-        yet_to_process = [(0, 0, (edge_to_start.face, edge_to_start.nex))]  # (error, iteration, (face, edge))
+        yet_to_process: list[tuple[int, int, tuple[Face, HalfEdge]]] = [
+            (0, 0, (edge_to_start.face, edge_to_start.nex))
+        ]  # (error, iteration, (face, edge))
         heapq.heapify(yet_to_process)
-        processed_faces = set()
-        err_dict = {edge_to_start.orig: 0, edge_to_start.dest: 0}
+        processed_faces: set[Face] = set()
+        err_dict: dict[Vertex, int] = {edge_to_start.orig: 0, edge_to_start.dest: 0}
         i = 0
         while yet_to_process:
             err, _, (f, initial) = heapq.heappop(yet_to_process)
@@ -1441,7 +1555,8 @@ def recompute_positions(self, edge_to_start: HalfEdge | None = None, faces: "set
                     )
                     err_dict[e.dest] = err_dict[e.orig] + 1
                 opposite_face = e.rev.face
-                if opposite_face in faces and opposite_face not in processed_faces:
+                if opposite_face in face_set and opposite_face not in processed_faces:
+                    assert opposite_face is not None
                     i += 1
                     heapq.heappush(
                         yet_to_process, (max(err_dict[e.orig], err_dict[e.dest]), i, (opposite_face, e.rev.nex))
@@ -1450,7 +1565,7 @@ def recompute_positions(self, edge_to_start: HalfEdge | None = None, faces: "set
                 if e is initial:
                     break
 
-    def construct_next_point(self, a: np.ndarray, b: np.ndarray, angle: float, length: float) -> np.ndarray:
+    def construct_next_point(self, a: NDArray, b: NDArray, angle: float, length: float) -> NDArray:
         """Construct the point ``c`` such that ``angle(a, b, c) == angle`` and ``|bc| == length``."""
         return self.geometry.construct_next_poly_point(a, b, angle, length)
         # next_angle = angle_to_axis(b - a) + np.pi - angle
@@ -1464,12 +1579,30 @@ def recompute_lengths_and_angles(self) -> None:
         for h in self.border_edges():
             h["length"] = h.rev["length"]
 
+    @overload
+    def get_position_view(
+        self,
+        vertices: list[Vertex] | None = ...,
+        *,
+        return_vertices: Literal[True] = ...,
+        position_key: str = ...,
+    ) -> tuple[NDArray, list[Vertex]]: ...
+
+    @overload
+    def get_position_view(
+        self,
+        vertices: list[Vertex] | None = ...,
+        *,
+        return_vertices: Literal[False],
+        position_key: str = ...,
+    ) -> NDArray: ...
+
     def get_position_view(
         self,
         vertices: list[Vertex] | None = None,
         return_vertices: bool = True,
         position_key: str = "pos",
-    ) -> "np.ndarray | tuple[np.ndarray, list[Vertex]]":
+    ) -> NDArray | tuple[NDArray, list[Vertex]]:
         """Return a single ``(N, d)`` array view of all vertex (and curve) positions.
 
         Mutating the returned array updates the underlying vertex/curve
@@ -1571,8 +1704,8 @@ def render(
         render_edges: bool = True,
         render_vertices: bool = True,
         for_cutting: bool = False,
-        **renderer_kwargs: object,
-    ) -> "Rendering":
+        **renderer_kwargs: Any,
+    ) -> Rendering:
         """Render the graph with Cairo and return an in-memory :class:`Rendering`.
 
         Args:
@@ -1605,7 +1738,7 @@ def render(
             for_cutting=for_cutting,
         )
 
-    def show(self, **style: object) -> None:
+    def show(self, **style: Any) -> None:
         """Render and display the graph (inline in Jupyter, a window in scripts).
 
         Args:
@@ -1613,7 +1746,7 @@ def show(self, **style: object) -> None:
         """
         self.render(**style).show()
 
-    def save(self, path: str, **style: object) -> None:
+    def save(self, path: str, **style: Any) -> None:
         """Render the graph and write it to *path*.
 
         ``path`` with no extension writes both ``path.svg`` and ``path.png``.
@@ -1629,14 +1762,14 @@ def central_face(self) -> Face:
         if self.geometry is not EuclideanGeometry:
             raise NotImplementedError
         fs = list(self.faces)
-        return fs[np.argmin([np.linalg.norm(f.midpoint()) for f in fs])]
+        return fs[int(np.argmin([np.linalg.norm(f.midpoint()) for f in fs]))]
 
     def central_vertex(self) -> Vertex:
         """Return the vertex closest to the origin (Euclidean only)."""
         if self.geometry is not EuclideanGeometry:
             raise NotImplementedError
         vs = list(self.vertices)
-        return vs[np.argmin([np.linalg.norm(v["pos"]) for v in vs])]
+        return vs[int(np.argmin([np.linalg.norm(v["pos"]) for v in vs]))]
 
 
 class EuclideanPositionHEG(GeometricHEG):
@@ -1645,7 +1778,7 @@ class EuclideanPositionHEG(GeometricHEG):
     Vertices carry ``pos`` attributes (numpy arrays of shape (2,)).  Provides epsilon-based vertex merging and position-aware graph operations.
     """
 
-    def __init__(self, **super_kwargs: object) -> None:
+    def __init__(self, **super_kwargs: Any) -> None:
         """Create a Euclidean-geometry half-edge graph."""
         super().__init__(geometry=EuclideanGeometry, **super_kwargs)
 
@@ -1653,7 +1786,7 @@ def __init__(self, **super_kwargs: object) -> None:
 # ------------------------------------------------ cyclic graph example ------------------------------------------------
 
 
-def rotate_by(list_like: "list | tuple", offset: "int | tuple[int, ...]") -> "list | zip":
+def rotate_by(list_like: list[Any] | tuple[Any, ...], offset: int | tuple[int, ...]) -> list[Any] | zip[Any]:
     """Cyclically rotate ``list_like`` by ``offset`` positions.
 
     Args:
@@ -1672,7 +1805,7 @@ def rotate_by(list_like: "list | tuple", offset: "int | tuple[int, ...]") -> "li
         return zip(*(rotate_by(list_like, off) for off in offset))
 
 
-def any_element(s):
+def any_element(s: Iterable[Any]) -> Any:
     """Return an arbitrary element from the iterable ``s``."""
     return next(iter(s))
 
@@ -1755,19 +1888,18 @@ def make_polygon_graph(positions: np.ndarray) -> "EuclideanPositionHEG":
         A :class:`EuclideanPositionHEG` whose only inner face is the polygon.
     """
     vs = [Vertex() for _ in range(len(positions))]
-    G = CyclicHalfedgeGraph(vs)
+    cyclic = CyclicHalfedgeGraph(vs)
     for v, p in zip(vs, positions):
         v["pos"] = p
-    G = EuclideanPositionHEG(other=G)
-    return G
+    return EuclideanPositionHEG(other=cyclic)
 
 
 class RegularNGon(CyclicHalfedgeGraph, InAngleHEG):
     """A regular *n*-gon as a half-edge graph with positions and angles."""
 
-    def __init__(self, n: int, *super_args: object, **super_kwargs: object) -> None:
+    def __init__(self, n: int, **super_kwargs: Any) -> None:
         """Construct a regular ``n``-gon with interior angle ``(n-2)/n * pi``."""
-        super(RegularNGon, self).__init__(vs=[Vertex() for _ in range(n)], *super_args, **super_kwargs)
+        super(RegularNGon, self).__init__(vs=[Vertex() for _ in range(n)], **super_kwargs)
         f = any_element(self.faces)
         for e in f.halfedge_iter():
             e["in_angle"] = (n - 2) / n * pi
diff --git a/eucare/image_to_graph.py b/eucare/image_to_graph.py
index 4bd9090..dce866c 100755
--- a/eucare/image_to_graph.py
+++ b/eucare/image_to_graph.py
@@ -9,11 +9,13 @@
 from __future__ import annotations
 
 import colorsys
+from typing import Any
 
 import mahotas as mh
 import networkx as nx
 import numpy as np
 from matplotlib import pyplot as plt
+from numpy.typing import NDArray
 from scipy import ndimage as ndi
 from skimage import color, io
 from skimage.transform import downscale_local_mean
@@ -23,11 +25,11 @@
 from .overlap import group_closeby
 
 
-def hsv_to_rgb(h, s, v):
+def hsv_to_rgb(h: float, s: float, v: float) -> NDArray[np.float32]:
     return np.array(colorsys.hsv_to_rgb(h, s, v), dtype=np.float32)
 
 
-def get_distinct_colors(n, min_sat=0.5, min_val=0.5):
+def get_distinct_colors(n: int, min_sat: float = 0.5, min_val: float = 0.5) -> NDArray[np.float32]:
     huePartition = 1.0 / (n + 1)
     hues = np.arange(0, n) * huePartition
     saturations = np.random.rand(n) * (1 - min_sat) + min_sat
@@ -35,7 +37,11 @@ def get_distinct_colors(n, min_sat=0.5, min_val=0.5):
     return np.stack([hsv_to_rgb(h, s, v) for h, s, v in zip(hues, saturations, values)], axis=0)
 
 
-def colorize_segmentation(seg, ignore_label=None, ignore_color=(0, 0, 0)):
+def colorize_segmentation(
+    seg: NDArray[Any],
+    ignore_label: int | None = None,
+    ignore_color: tuple[float, float, float] = (0, 0, 0),
+) -> NDArray[Any]:
     assert isinstance(seg, np.ndarray)
     assert seg.dtype.kind in ("u", "i")
     if ignore_label is not None:
@@ -49,7 +55,13 @@ def colorize_segmentation(seg, ignore_label=None, ignore_color=(0, 0, 0)):
     return result
 
 
-def plot_image(image, figheight=5, title=None, colorbar=False, **kwargs):
+def plot_image(
+    image: NDArray[Any],
+    figheight: float = 5,
+    title: str | None = None,
+    colorbar: bool = False,
+    **kwargs: Any,
+) -> None:
     plt.figure(figsize=(figheight * image.shape[1] / image.shape[0], figheight))
     im = plt.imshow(image, **kwargs)
     if colorbar:
@@ -117,7 +129,7 @@ def endPoints(skel: np.ndarray) -> np.ndarray:
     return ep
 
 
-def pruning(skeleton, size=None):
+def pruning(skeleton: NDArray[Any], size: int | None = None) -> NDArray[Any]:
     """remove iteratively end points "size"
     times from the skeleton
     """
@@ -203,14 +215,14 @@ def image_to_graph(
     edge_labels, n_edges = ndi.label(skeleton ^ branching_points)
 
     plot_image(
-        colorize_segmentation(ndi.grey_dilation(edge_labels, size=3).astype(np.int32), ignore_label=0),
+        colorize_segmentation(ndi.grey_dilation(edge_labels, size=(3, 3)).astype(np.int32), ignore_label=0),
         figheight=5,
         title="edges",
     )
     plt.show()
 
     # construct the graph: for every branching point, find the edges connected to it
-    edge_dict = {}
+    edge_dict: dict[int, set[int]] = {}
     for bp in range(1, n_branch_points + 1):
         adjacent_edges = np.unique(edge_labels[ndi.binary_dilation(branch_point_labels == bp)])
         adjacent_edges = adjacent_edges[adjacent_edges > 0]
@@ -251,7 +263,7 @@ def image_to_graph(
         (new_bp_labels[i - 1], new_bp_labels[j - 1]) for i, j in edges if new_bp_labels[i - 1] != new_bp_labels[j - 1]
     )
 
-    graph = nx.Graph()
+    graph: nx.Graph = nx.Graph()
     graph.add_edges_from(new_edges)
 
     # convert to eucare graph
diff --git a/eucare/instructions.py b/eucare/instructions.py
index effd42e..0bf19bd 100755
--- a/eucare/instructions.py
+++ b/eucare/instructions.py
@@ -2,14 +2,14 @@
 
 from __future__ import annotations
 
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, Callable
 
 if TYPE_CHECKING:
     from .half import HalfEdge, HalfEdgeGraph
     from .prototiles import ProtoTile
 
 
-def attatch_tile_instruction(proto_tile: "ProtoTile", label: object = None) -> "callable":
+def attatch_tile_instruction(proto_tile: ProtoTile, label: object = None) -> Callable[[HalfEdgeGraph, HalfEdge], None]:
     """Return a callable that builds a fresh tile graph from ``proto_tile`` and glues it to a given edge.
 
     The newly created tile graph may itself carry instructions on its border
diff --git a/eucare/intersecting_cylinders/mesh3d.py b/eucare/intersecting_cylinders/mesh3d.py
index ebd7ccf..13cf1c4 100644
--- a/eucare/intersecting_cylinders/mesh3d.py
+++ b/eucare/intersecting_cylinders/mesh3d.py
@@ -72,7 +72,7 @@
 from .profiles import Profile
 
 if TYPE_CHECKING:
-    from ..half import EuclideanPositionHEG
+    from ..half import EuclideanPositionHEG, Face
 
 
 def _spike_depth_from_profile(
@@ -196,12 +196,12 @@ def _build_ortho_with_tangent_points(
     return G_ortho
 
 
-def _vertex_circle_radii(G_ortho: "EuclideanPositionHEG") -> dict:
+def _vertex_circle_radii(G_ortho: "EuclideanPositionHEG") -> dict[half.Vertex, float]:
     """Return ``{original_vertex: r_v}`` where ``r_v = |v - t|`` is the radius
     of the blue circle centered at the original vertex (its distance to any of
     its incident edge-tangent points after the position fix).
     """
-    radii: dict = {}
+    radii: dict[half.Vertex, float] = {}
     for v_o in G_ortho.vertices:
         if "pre_conway" not in v_o.attributes:
             continue
@@ -220,7 +220,9 @@ def _vertex_circle_radii(G_ortho: "EuclideanPositionHEG") -> dict:
     return radii
 
 
-def _classify_ortho_quad(face) -> tuple | None:
+def _classify_ortho_quad(
+    face: "Face",
+) -> tuple[half.Vertex, half.Vertex, half.Vertex, half.Vertex] | None:
     """Return ``(v_corner, c_corner, t1_corner, t2_corner)`` for a 4-corner ortho
     quad, or ``None`` if it does not have the expected (V, E, F, E) structure.
 
diff --git a/eucare/intersecting_cylinders/pipeline.py b/eucare/intersecting_cylinders/pipeline.py
index a5ce4a6..975f358 100644
--- a/eucare/intersecting_cylinders/pipeline.py
+++ b/eucare/intersecting_cylinders/pipeline.py
@@ -6,7 +6,7 @@
 
 from __future__ import annotations
 
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, Any, TypeVar
 
 import numpy as np
 
@@ -16,18 +16,21 @@
 if TYPE_CHECKING:
     from ..half import EuclideanPositionHEG
 
+_K = TypeVar("_K")
+_V = TypeVar("_V")
+
 # RGB tuples used throughout the pipeline as ``color_key`` values.
 _BLACK = (0.0, 0.0, 0.0)
 _BLUE = (0.0, 0.0, 1.0)
 _RED = (1.0, 0.0, 0.0)
 
 
-def _reverse_mapping(d: dict) -> dict:
+def _reverse_mapping(d: dict[_K, _V]) -> dict[_V, _K]:
     """Return ``{v: k for k, v in d.items()}``."""
     return {value: key for key, value in d.items()}
 
 
-def _from_pre(obj) -> bool:
+def _from_pre(obj: Any) -> bool:
     """True iff ``obj`` carries a ``pre_conway`` attribute (set by Conway operators)."""
     return "pre_conway" in obj.attributes
 
@@ -128,6 +131,7 @@ def make_intersecting_cylinders(
             p = v["pos"]
             p2 = v2["pos"]
 
+            assert h_orig.face is not None
             c = h_orig.face.midpoint()
             hc = base.project_to_line(np.stack([p, p2]), c)
 
diff --git a/eucare/intersecting_cylinders/profiles.py b/eucare/intersecting_cylinders/profiles.py
index 0ebbc50..92e423d 100644
--- a/eucare/intersecting_cylinders/profiles.py
+++ b/eucare/intersecting_cylinders/profiles.py
@@ -21,7 +21,7 @@
 from __future__ import annotations
 
 from dataclasses import dataclass
-from typing import Callable
+from typing import Any, Callable
 
 import numpy as np
 from numpy.typing import NDArray
@@ -107,7 +107,7 @@ def from_function(
             shrink_factor=shrink_factor,
         )
 
-    def plot(self, ax=None) -> None:
+    def plot(self, ax: Any = None) -> None:
         """Plot the simplified ``(t, y)`` polyline."""
         import matplotlib.pyplot as plt
 
diff --git a/eucare/intersecting_cylinders/show_circle_packings.py b/eucare/intersecting_cylinders/show_circle_packings.py
index 6af3520..e120c02 100644
--- a/eucare/intersecting_cylinders/show_circle_packings.py
+++ b/eucare/intersecting_cylinders/show_circle_packings.py
@@ -56,13 +56,13 @@ def build_dual_circle_packings(G: "EuclideanPositionHEG") -> "EuclideanPositionH
 
     for v in G_ortho.vertices.union(G_ortho.faces):
         pre = v.get("pre_conway")
-        if isinstance(pre, half.Face):
+        if isinstance(pre, half.Face) and isinstance(v, half.Vertex):
             v["color_key"] = _FACE_COLOR
             v["vertex_radius"] = float(np.linalg.norm(v["pos"] - v.any_outgoing.dest["pos"]))
             for h in v.outgoing_iter():
                 h["line_width"] = 0.0
                 h.rev["line_width"] = 0.0
-        elif isinstance(pre, half.Vertex):
+        elif isinstance(pre, half.Vertex) and isinstance(v, half.Vertex):
             v["color_key"] = _VERTEX_COLOR
             v["vertex_radius"] = float(np.linalg.norm(v["pos"] - v.any_outgoing.dest["pos"]))
             for h in v.outgoing_iter():
@@ -78,7 +78,7 @@ def build_dual_circle_packings(G: "EuclideanPositionHEG") -> "EuclideanPositionH
     return G_ortho
 
 
-def show_dual_circle_packings(G: "EuclideanPositionHEG", **show_kwargs: Any):
+def show_dual_circle_packings(G: "EuclideanPositionHEG", **show_kwargs: Any) -> Any:
     """Render the two dual circle packings of a tiling.
 
     This is a thin convenience wrapper around :func:`build_dual_circle_packings`
diff --git a/eucare/intersecting_cylinders/triangle_twist.py b/eucare/intersecting_cylinders/triangle_twist.py
index 398b92e..c96eb23 100644
--- a/eucare/intersecting_cylinders/triangle_twist.py
+++ b/eucare/intersecting_cylinders/triangle_twist.py
@@ -7,7 +7,9 @@
 
 from __future__ import annotations
 
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, Any
+
+import numpy as np
 
 from .. import base, half
 
@@ -20,7 +22,7 @@
 _YELLOW = (1.0, 1.0, 0.0)
 
 
-def _set_color(color, *objs) -> None:
+def _set_color(color: tuple[float, float, float], *objs: Any) -> None:
     """Assign ``color_key = color`` to every element in ``objs`` (recursive over iterables)."""
     for obj in objs:
         if isinstance(obj, (half.Face, half.HalfEdge, half.Vertex)):
@@ -43,15 +45,16 @@ def convert_to_triangle_twist(G: "EuclideanPositionHEG", v: "Vertex") -> None:
     """
     G.delete_subset([h for h in v.outgoing_iter() if h["color_key"] == _RED])
 
-    tasks = []
+    tasks: list[tuple[half.HalfEdge, half.Face, half.Vertex, np.ndarray]] = []
     for h in list(v.outgoing_iter()):
         v2 = h.rev.nex.nex.dest
         _set_color(_GREEN, v2)
         direction = v["pos"] - h.rev.nex.dest["pos"]
         f = h.rev.face
+        assert f is not None
         pos = base.line_intersection(
-            [v2["pos"], v2["pos"] + direction],
-            [h.orig["pos"], h.dest["pos"]],
+            np.stack([v2["pos"], v2["pos"] + direction]),
+            np.stack([h.orig["pos"], h.dest["pos"]]),
         )
         tasks.append((h, f, v2, pos))
 
@@ -61,7 +64,9 @@ def convert_to_triangle_twist(G: "EuclideanPositionHEG", v: "Vertex") -> None:
         _set_color(_RED, h2)
 
     for h in v.outgoing_iter():
-        h2, _ = G.subdivide_face(h.face, h.dest, h.pre.orig)
+        f_h = h.face
+        assert f_h is not None
+        h2, _ = G.subdivide_face(f_h, h.dest, h.pre.orig)
         _set_color(_BLUE, h2)
 
     G.delete_subset([v])
diff --git a/eucare/io.py b/eucare/io.py
index a73f690..507af0d 100755
--- a/eucare/io.py
+++ b/eucare/io.py
@@ -3,19 +3,20 @@
 from __future__ import annotations
 
 import os
+from typing import Any, Callable
 
 import numpy as np
 import yaml
 
 import eucare as ec
 
-from .half import Face, HalfEdge, HalfEdgeGraph, Vertex
+from .half import AttributeObject, Face, HalfEdge, HalfEdgeGraph, Vertex
 
 
 def graph_to_dict(
     G: HalfEdgeGraph,
     attributes_to_save: tuple[str, ...] = ("pos", "length", "in_angle", "color_key"),
-) -> dict:
+) -> dict[str, Any]:
     """Serialise a half-edge graph to a JSON/YAML-friendly nested dict.
 
     Vertices, half-edges, and faces are each given an opaque string label
@@ -36,25 +37,27 @@ def graph_to_dict(
     halfedge_labels = {h: f"h{i}" for i, h in enumerate(G.halfedges)}
     face_labels = {f: f"f{i}" for i, f in enumerate(G.faces)}
 
-    labels = {None: None}
+    labels: dict[Any, str | None] = {None: None}
     labels.update(vertex_labels)
     labels.update(halfedge_labels)
     labels.update(face_labels)
 
-    def represent_attributes(obj):
-        result = {}
+    def represent_attributes(obj: AttributeObject) -> dict[str, Any]:
+        result: dict[str, Any] = {}
         for attr in attributes_to_save:
             if attr in obj.attributes:
-                value = obj[attr]
+                value: Any = obj[attr]
                 if isinstance(value, np.ndarray):
                     value = value.tolist()
-                if np.isscalar(value):
+                if isinstance(value, (int, float, np.floating, np.integer)):
                     value = float(value)
                 result[attr] = value
         return result
 
-    def add_attributes(func):
-        def wrapped(obj):
+    def add_attributes(
+        func: Callable[[Any], dict[str, Any]],
+    ) -> Callable[[Any], dict[str, Any]]:
+        def wrapped(obj: Any) -> dict[str, Any]:
             result = func(obj)
             attrs = represent_attributes(obj)
             if attrs:
@@ -64,11 +67,11 @@ def wrapped(obj):
         return wrapped
 
     @add_attributes
-    def represent_vertex(v):
+    def represent_vertex(v: Vertex) -> dict[str, Any]:
         return dict(any_outgoing=labels[v.any_outgoing])
 
     @add_attributes
-    def represent_halfedge(h):
+    def represent_halfedge(h: HalfEdge) -> dict[str, Any]:
         return dict(
             orig=labels[h.orig],
             dest=labels[h.dest],
@@ -79,7 +82,7 @@ def represent_halfedge(h):
         )
 
     @add_attributes
-    def represent_face(f):
+    def represent_face(f: Face) -> dict[str, Any]:
         return dict(any_side=labels[f.any_side])
 
     vertex_dict = {label: represent_vertex(v) for v, label in vertex_labels.items()}
@@ -90,15 +93,15 @@ def represent_face(f):
     return graph_dict
 
 
-def dict_to_graph(graph_dict: dict) -> ec.half.EuclideanPositionHEG:
+def dict_to_graph(graph_dict: dict[str, Any]) -> ec.half.EuclideanPositionHEG:
     """Inverse of :func:`graph_to_dict`: reconstruct a graph from its serialised dict.
 
     The returned graph is always an :class:`EuclideanPositionHEG` regardless of
     the source graph's class (TODO: persist the class).
     """
 
-    def unwrap_attributes(obj_dict):
-        result = {}
+    def unwrap_attributes(obj_dict: dict[str, Any]) -> dict[str, Any]:
+        result: dict[str, Any] = {}
         for key, value in obj_dict.pop("attributes", {}).items():
             if isinstance(value, list):
                 try:
@@ -108,12 +111,12 @@ def unwrap_attributes(obj_dict):
             result[key] = value
         return result
 
-    lookup = {None: None}
+    lookup: dict[str | None, Any] = {None: None}
     # create the halfedges
     for label in graph_dict["halfedges"]:
         lookup[label] = HalfEdge()
 
-    vs = set()
+    vs: set[Vertex] = set()
     for label, v_dict in graph_dict["vertices"].items():
         attrs = unwrap_attributes(v_dict)
         v_dict["any_outgoing"] = lookup[v_dict["any_outgoing"]]
@@ -122,7 +125,7 @@ def unwrap_attributes(obj_dict):
         lookup[label] = v
         vs.add(v)
 
-    fs = set()
+    fs: set[Face] = set()
     for label, f_dict in graph_dict["faces"].items():
         attrs = unwrap_attributes(f_dict)
         f_dict["any_side"] = lookup[f_dict["any_side"]]
@@ -131,7 +134,7 @@ def unwrap_attributes(obj_dict):
         lookup[label] = f
         fs.add(f)
 
-    hs = set()
+    hs: set[HalfEdge] = set()
     for label, h_dict in graph_dict["halfedges"].items():
         attrs = unwrap_attributes(h_dict)
         h = lookup[label]
diff --git a/eucare/layout.py b/eucare/layout.py
index cb00635..eca559f 100755
--- a/eucare/layout.py
+++ b/eucare/layout.py
@@ -46,6 +46,7 @@ def _angle_to_height(angle: float) -> float:
 def rotate_graph(G: "GeometricHEG", angle: float) -> None:
     """Rotate every vertex of *G* by *angle* in place."""
     ps = G.get_position_view(return_vertices=False)
+    assert isinstance(ps, np.ndarray)
     ps[:] = ps @ np.array([[np.cos(angle), np.sin(angle)], [-np.sin(angle), np.cos(angle)]])
 
 
diff --git a/eucare/marching_cubes.py b/eucare/marching_cubes.py
index cd0e0ae..e743e4b 100644
--- a/eucare/marching_cubes.py
+++ b/eucare/marching_cubes.py
@@ -9,9 +9,11 @@
 
 import logging
 from numbers import Number
+from typing import Any, Callable
 
 import numpy as np
 import sdf
+from numpy.typing import NDArray
 
 logger = logging.getLogger(__name__)
 
@@ -32,7 +34,12 @@
 class CartesianGrid:
     """Axis aligned cartesian grid in n-dimensional euclidean space"""
 
-    def __init__(self, ndim: int, origin: np.ndarray = None, cell_size: np.ndarray = None) -> None:
+    def __init__(
+        self,
+        ndim: int,
+        origin: NDArray[Any] | None = None,
+        cell_size: NDArray[Any] | None = None,
+    ) -> None:
         self.ndim = ndim
         self.origin = np.zeros(self.ndim) if origin is None else origin
         self.cell_size = np.ones(self.ndim, dtype=np.float32) if cell_size is None else cell_size
@@ -131,8 +138,10 @@ def subdivide_indices(self, indices: np.ndarray, subdivisions: int | tuple | np.
 
 
 def oct_tree_marching_cubes(
-    f: callable, step: np.ndarray, bounds: np.ndarray | None = None
-) -> tuple[np.ndarray, CartesianGrid]:
+    f: Callable[[NDArray[Any]], NDArray[Any]],
+    step: NDArray[Any],
+    bounds: NDArray[Any] | None = None,
+) -> tuple[NDArray[Any], CartesianGrid]:
     """Locate boundary cells of an SDF using an oct-tree refinement strategy.
 
     Builds successively coarser grids until the bounding box fits in one
@@ -170,8 +179,8 @@ def oct_tree_marching_cubes(
     level = 0
     block_indices = [np.arange(n) for n in box_size // step]
 
-    def block_indices_n_cells(block_indices):
-        return np.product([len(indices) for indices in block_indices])
+    def block_indices_n_cells(block_indices: list[NDArray[Any]]) -> int:
+        return int(np.prod([len(indices) for indices in block_indices]))
 
     n_cells_naive = block_indices_n_cells(block_indices)
     logger.info("n_cells_naive=%d", n_cells_naive)
diff --git a/eucare/overlap.py b/eucare/overlap.py
index 6fdbc5e..f05ee38 100755
--- a/eucare/overlap.py
+++ b/eucare/overlap.py
@@ -8,16 +8,18 @@
 import sys
 import tempfile
 from collections import defaultdict
+from collections.abc import Iterable, Iterator
 from contextlib import contextmanager
 from copy import copy
 from functools import cmp_to_key
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, Any, Callable
 
 import networkx as nx
 import numpy as np
 import pulp
 from fastcluster import linkage_vector
 from numba import jit
+from numpy.typing import NDArray
 from scipy.cluster.hierarchy import fcluster
 from tqdm.auto import tqdm
 
@@ -27,7 +29,7 @@
 from .rendering import BORDER_COLOR, MOUNTAIN_COLOR, VALLEY_COLOR, CairoRenderer, SvgwriteRenderer, multi_show
 
 if TYPE_CHECKING:
-    from .half import EuclideanPositionHEG, Face
+    from .half import EuclideanPositionHEG, Face, HalfEdge
 
 logger = logging.getLogger(__name__)
 
@@ -37,7 +39,7 @@ def intervals_overlapping(interval1: tuple[float, float], interval2: tuple[float
     return not (interval1[0] > interval2[1]) and not (interval1[1] < interval2[0])
 
 
-def get_potential_intersections(segments: np.ndarray, epsilon: float = 1e-12) -> list[tuple[int, int]]:
+def get_potential_intersections(segments: NDArray[Any], epsilon: float = 1e-12) -> list[tuple[int, int]]:
     """Return pairs of segment indices that may intersect, using a sweep-line approach on bounding boxes."""
     # list of start and end points, with index of corresponding segment and flag whether it is a start or an end point.
     segments = np.array(segments).copy()
@@ -50,19 +52,19 @@ def get_potential_intersections(segments: np.ndarray, epsilon: float = 1e-12) ->
     x_coords.sort(axis=1)
 
     # create tuples (position, index, is_start)
-    points = [(s[0][0] - epsilon, i, 1) for i, s in enumerate(segments)] + [
+    point_list: list[tuple[Any, int, int]] = [(s[0][0] - epsilon, i, 1) for i, s in enumerate(segments)] + [
         (s[1][0], i, 0) for i, s in enumerate(segments)
     ]
-    points = np.array(points, dtype=tuple)
+    points: NDArray[Any] = np.array(point_list, dtype=tuple)
 
     # sort by first coordinate
     points = points[np.argsort(points[:, 0])]
-    active_labels = set()
-    possibly_intersecting = list()
+    active_labels: set[int] = set()
+    possibly_intersecting: list[tuple[int, int]] = list()
     for i, is_start in points[:, 1:]:
         if is_start:
             for j in active_labels:
-                if intervals_overlapping(segments[i, :, 1], segments[j, :, 1]):
+                if intervals_overlapping(tuple(segments[i, :, 1]), tuple(segments[j, :, 1])):
                     possibly_intersecting.append((i, j))
             active_labels.add(i)
         else:
@@ -71,19 +73,19 @@ def get_potential_intersections(segments: np.ndarray, epsilon: float = 1e-12) ->
 
 
 @jit(nopython=True)
-def _on_segment(p, q, r):
+def _on_segment(p: NDArray[Any], q: NDArray[Any], r: NDArray[Any]) -> bool:
     """Check if q lies on segment pr, assuming the three points are collinear."""
     return (min(p[0], r[0]) <= q[0] <= max(p[0], r[0])) and (min(p[1], r[1]) <= q[1] <= max(p[1], r[1]))
 
 
 @jit(nopython=True)
-def _det(a, b):
+def _det(a: Any, b: Any) -> float:
     """Compute the 2x2 determinant of vectors a and b."""
     return a[0] * b[1] - a[1] * b[0]
 
 
 @jit(nopython=True)
-def line_segment_intersections(s1, s2, eps=1e-12):
+def line_segment_intersections(s1: NDArray[Any], s2: NDArray[Any], eps: float = 1e-12) -> list[NDArray[Any]]:
     """Return a list of intersection points between two line segments, handling collinear cases."""
 
     if s1.dtype is not np.float64:
@@ -140,22 +142,23 @@ def line_segment_intersections(s1, s2, eps=1e-12):
     return result
 
 
-def fast_group_closeby(pts: np.ndarray, eps: float) -> np.ndarray:
+def fast_group_closeby(pts: NDArray[Any], eps: float) -> NDArray[np.int32]:
     """Cluster nearby points using single-linkage clustering on cityblock distance."""
     Z = linkage_vector(pts, method="single", metric="cityblock")
     unsorted = fcluster(Z, eps, criterion="distance")
-    result = []
+    result: list[int] = []
     i = 0
-    id_map = {}
+    id_map: dict[int, int] = {}
     for idx in unsorted:
-        if idx not in id_map:
-            id_map[idx] = i
+        idx_int = int(idx)
+        if idx_int not in id_map:
+            id_map[idx_int] = i
             i += 1
-        result.append(id_map[idx])
+        result.append(id_map[idx_int])
     return np.array(result, dtype=np.int32)
 
 
-def faster_group_closeby_nx(arr: np.ndarray, eps: float) -> np.ndarray:
+def faster_group_closeby_nx(arr: NDArray[Any], eps: float) -> NDArray[np.int32]:
     """Cluster nearby points via approximate single-linkage using grid-based connected components.
 
     Points less than eps apart receive the same label. Return an array mapping
@@ -164,7 +167,7 @@ def faster_group_closeby_nx(arr: np.ndarray, eps: float) -> np.ndarray:
     if len(arr) == 0:
         return np.zeros(0, dtype=np.int32)
     scaled = arr / eps / 2
-    G = nx.Graph()
+    G: nx.Graph = nx.Graph()
     for offset in itertools.product([0, 0.5], repeat=arr.shape[-1]):
         _, index, inverse_ids, counts = np.unique(
             np.floor(scaled + offset), axis=0, return_counts=True, return_inverse=True, return_index=True
@@ -221,9 +224,9 @@ def overlap_graph(G: "EuclideanPositionHEG", eps: float = 1e-10) -> "EuclideanPo
     # TODO implement sweeping line algorithm https://www.geeksforgeeks.org/given-a-set-of-line-segments-find-if-any-two-segments-intersect/
 
     # This will be a list of the positions of crossings between line segments. This includes the original vertices.
-    crossings = []
+    crossings_list: list[NDArray[Any]] = []
     # This will be a list mapping the index of a crossing to the pair (i, j) of indices of the corresponding edges.
-    crossings_to_edges = []
+    crossings_to_edges: list[tuple[int, int]] = []
 
     # iterate over all pairs
     potential_intersections = get_potential_intersections(line_segments, epsilon=eps)
@@ -233,10 +236,10 @@ def overlap_graph(G: "EuclideanPositionHEG", eps: float = 1e-10) -> "EuclideanPo
         intersections = line_segment_intersections(l1, l2, eps=eps)
         if not intersections:
             continue
-        crossings.extend(intersections)
+        crossings_list.extend(intersections)
         for _ in range(len(intersections)):
             crossings_to_edges.append((i, j))
-    crossings = np.array(crossings)
+    crossings: NDArray[Any] = np.array(crossings_list)
     timer.round("crossings")
 
     # ------ group closeby crossings ------
@@ -254,38 +257,38 @@ def overlap_graph(G: "EuclideanPositionHEG", eps: float = 1e-10) -> "EuclideanPo
 
     # construct an array mapping crossings to all involved edges,
     # and on mapping edges to all crossings they are involved in.
-    filtered_crossings_to_edges = [set() for i in range(n_filtered_crossings)]  # TODO: get rid of
-    edges_to_crossings = [set() for i in range(len(edges))]
+    filtered_crossings_to_edges: list[set[int]] = [set() for i in range(n_filtered_crossings)]  # TODO: get rid of
+    edges_to_crossings: list[set[int]] = [set() for i in range(len(edges))]
     for i, edge_ids in enumerate(crossings_to_edges):
         filtered_crossings_to_edges[clustering[i]].update(edge_ids)
-        for e in edge_ids:
-            edges_to_crossings[e].add(clustering[i])
+        for edge_idx in edge_ids:
+            edges_to_crossings[edge_idx].add(clustering[i])
     timer.round("filtered crossings")
 
     # ------ get crossing orders, construct nx graph ------
     logger.info("getting crossing orders..")
     # nodes = np.arange(n_filtered_crossings)
-    nx_edges = []
-    edge_to_ordered_ids = []
+    nx_edges_list: list[NDArray[Any]] = []
+    edge_to_ordered_ids: list[NDArray[Any]] = []
     # for each edge (=line segment), order the crossings on it from its orig to dest.
     for i in tqdm(range(len(line_segments))):
-        e = edges[i]
+        he = edges[i]
         crossing_ids = np.array(list(edges_to_crossings[i]))
         crossing_positions = filtered_crossings[crossing_ids]
         progression_along_edge = (
-            (crossing_positions - e.orig["pos"][None]) * (e.dest["pos"] - e.orig["pos"])[None]
+            (crossing_positions - he.orig["pos"][None]) * (he.dest["pos"] - he.orig["pos"])[None]
         ).sum(-1)
         order = np.argsort(progression_along_edge)
         ordered_ids = crossing_ids[order]
         edge_to_ordered_ids.append(ordered_ids)
         # add line segments between crossings along the edge to the nx graph
-        nx_edges.append(np.stack([ordered_ids[:-1], ordered_ids[1:]], axis=-1))
-    nx_edges = np.concatenate(nx_edges)
+        nx_edges_list.append(np.stack([ordered_ids[:-1], ordered_ids[1:]], axis=-1))
+    nx_edges: NDArray[Any] = np.concatenate(nx_edges_list)
     timer.round("crossing orders")
 
     logger.info("constructing nx graph..")
-    nx_graph = nx.Graph()
-    nx_graph.add_edges_from(nx_edges)
+    nx_graph: nx.Graph = nx.Graph()
+    nx_graph.add_edges_from([(int(a), int(b)) for a, b in nx_edges])
     nx_positions = {i: pos for i, pos in enumerate(filtered_crossings)}
     timer.round("nx graph constructed.")
     logger.info("order of nx graph: %d", nx_graph.order())
@@ -311,7 +314,9 @@ def overlap_graph(G: "EuclideanPositionHEG", eps: float = 1e-10) -> "EuclideanPo
         f["original_faces"] = set()
 
     # get face orientations
-    face_orientations = {}
+    from .half import Face as _Face
+
+    face_orientations: dict[_Face, bool] = {}
     for f in G.faces:
         face_orientations[f] = f.area() > 0
 
@@ -319,9 +324,9 @@ def overlap_graph(G: "EuclideanPositionHEG", eps: float = 1e-10) -> "EuclideanPo
         ordered_ids_0 = edge_to_ordered_ids[i]
         if len(ordered_ids_0) == 0:
             raise RuntimeError(f"orig and dest of every edge should be crossings, but edge {edges[i]} has no crossings")
-        for e in (edges[i], edges[i].rev):
+        for he2 in (edges[i], edges[i].rev):
             ordered_ids = ordered_ids_0
-            if e is edges[i].rev:
+            if he2 is edges[i].rev:
                 # reverse order of crossings (=nodes in nx graph) for reversed edge
                 ordered_ids = ordered_ids[::-1]
             # this for loop sets adds the orig of the original edge to the 'original_vertices' attribute
@@ -331,11 +336,13 @@ def overlap_graph(G: "EuclideanPositionHEG", eps: float = 1e-10) -> "EuclideanPo
             #  along the edge is no vertex in the overlap graph. Show an example!
             for idx in ordered_ids:
                 if idx in v_lookup and "original_vertices" in v_lookup[idx]:
-                    v_lookup[idx]["original_vertices"].add(e.orig)
+                    v_lookup[idx]["original_vertices"].add(he2.orig)
                     break
             # to assign original edges and face_groups, revert the order of the edge if the face is oriented negatively.
-            if not e.on_border() and not face_orientations[e.face]:
-                ordered_ids = ordered_ids[::-1]
+            if not he2.on_border():
+                assert he2.face is not None
+                if not face_orientations[he2.face]:
+                    ordered_ids = ordered_ids[::-1]
             # iterate along segments between crossings on the original edge
             for k, l in zip(ordered_ids[:-1], ordered_ids[1:]):
                 if k not in v_lookup or l not in v_lookup:
@@ -343,16 +350,17 @@ def overlap_graph(G: "EuclideanPositionHEG", eps: float = 1e-10) -> "EuclideanPo
                 # find the halfedge in the overlap graph between crossings k and l
 
                 new_edge = next(h for h in v_lookup[k].outgoing_iter() if h.dest is v_lookup[l])
-                new_edge["original_edges"].add(e)
-                if not e.on_border():
-                    if frozenset((e, e.rev)) not in new_edge["original_face_groups"]:
-                        new_edge["original_face_groups"][frozenset((e, e.rev))] = set()
-                    new_edge["original_face_groups"][frozenset((e, e.rev))].add(e.face)
+                new_edge["original_edges"].add(he2)
+                if not he2.on_border():
+                    assert he2.face is not None
+                    if frozenset((he2, he2.rev)) not in new_edge["original_face_groups"]:
+                        new_edge["original_face_groups"][frozenset((he2, he2.rev))] = set()
+                    new_edge["original_face_groups"][frozenset((he2, he2.rev))].add(he2.face)
 
     # convert the face groups from a dict to a list of sets, as the information which edge the individual original faces
     # were coming from is no longer required.
-    for e in overlap_G.halfedges:
-        e["original_face_groups"] = [frozenset(group) for group in e["original_face_groups"].values()]
+    for he3 in overlap_G.halfedges:
+        he3["original_face_groups"] = [frozenset(group) for group in he3["original_face_groups"].values()]
 
     # ------ assign original faces ------
     logger.info("assigning original faces..")
@@ -360,24 +368,28 @@ def overlap_graph(G: "EuclideanPositionHEG", eps: float = 1e-10) -> "EuclideanPo
     # the folded model).
     initial_edge = next(overlap_G.border_edge_iter()).rev
     assert not initial_edge.on_border()
-    frontier = [(initial_edge, {e_orig.face for e_orig in initial_edge["original_edges"] if not e_orig.on_border()})]
+    frontier: list[tuple[HalfEdge, set[_Face]]] = [
+        (initial_edge, {e_orig.face for e_orig in initial_edge["original_edges"] if not e_orig.on_border()})
+    ]
     yet_to_assign = set(overlap_G.faces)
     while yet_to_assign:
         current_halfedge, original_faces = frontier.pop()
         assert not current_halfedge.on_border()
         current_face = current_halfedge.face
+        assert current_face is not None
         if current_face not in yet_to_assign:
             continue
         current_face["original_faces"] = original_faces
         yet_to_assign.remove(current_face)
-        for e in current_face.halfedge_iter():
-            if not e.rev.on_border() and e.rev.face in yet_to_assign:
+        for he4 in current_face.halfedge_iter():
+            if not he4.rev.on_border() and he4.rev.face in yet_to_assign:
+                assert he4.rev.face is not None
                 frontier.append(
                     (
-                        e.rev,
+                        he4.rev,
                         (
-                            original_faces - {e_orig.face for e_orig in e["original_edges"] if not e_orig.on_border()}
-                        ).union({e_orig.face for e_orig in e.rev["original_edges"] if not e_orig.on_border()}),
+                            original_faces - {e_orig.face for e_orig in he4["original_edges"] if not e_orig.on_border()}
+                        ).union({e_orig.face for e_orig in he4.rev["original_edges"] if not e_orig.on_border()}),
                     )
                 )
     timer.round("complete")
@@ -394,9 +406,9 @@ def overlap_graph(G: "EuclideanPositionHEG", eps: float = 1e-10) -> "EuclideanPo
 SORTED_ORIGINAL_FACES = "sorted_original_faces"
 
 
-def find_triplet_overlap_areas(G: "EuclideanPositionHEG") -> dict[frozenset, float]:
+def find_triplet_overlap_areas(G: "EuclideanPositionHEG") -> dict[frozenset[Any], float]:
     """Compute the total overlap area for each triplet of original faces."""
-    result = defaultdict(float)
+    result: dict[frozenset[Any], float] = defaultdict(float)
     bar = tqdm(G.faces, desc="finding triplet overlap areas")
     for f in bar:
         area = f.area()
@@ -406,12 +418,14 @@ def find_triplet_overlap_areas(G: "EuclideanPositionHEG") -> dict[frozenset, flo
     return result
 
 
-def find_fold_over_facet_lengths(G: "EuclideanPositionHEG") -> dict[tuple[frozenset, "Face"], float]:
+def find_fold_over_facet_lengths(G: "EuclideanPositionHEG") -> dict[tuple[frozenset[Any], "Face"], float]:
     """Compute total edge length where a fold passes over a facet."""
-    result = defaultdict(float)
+    result: dict[tuple[frozenset[Any], Face], float] = defaultdict(float)
     for e in tqdm(G.halfedges, desc="finding fold over facet lengths"):
         if e.on_border() or e.rev.on_border():
             continue
+        assert e.face is not None
+        assert e.rev.face is not None
         over_both = e.face[ORIGINAL_FACES].intersection(e.rev.face[ORIGINAL_FACES])
         edge_length = e[LENGTH]
         for group in e[FACES_OF_FOLDS_ON_EDGE]:
@@ -422,9 +436,9 @@ def find_fold_over_facet_lengths(G: "EuclideanPositionHEG") -> dict[tuple[frozen
     return result
 
 
-def find_conincident_fold_lengths(G: "EuclideanPositionHEG") -> dict[frozenset, float]:
+def find_conincident_fold_lengths(G: "EuclideanPositionHEG") -> dict[frozenset[Any], float]:
     """Compute total edge length where two folds coincide."""
-    result = defaultdict(float)
+    result: dict[frozenset[Any], float] = defaultdict(float)
     for e in tqdm(G.halfedges, desc="finding coincident fold lengths"):
         if e.on_border():
             continue
@@ -436,12 +450,12 @@ def find_conincident_fold_lengths(G: "EuclideanPositionHEG") -> dict[frozenset,
     return dict(result)
 
 
-def cache_all(cache=None):
+def cache_all(cache: dict[Any, Any] | None = None) -> Callable[[Callable[..., Any]], Callable[..., Any]]:
     """Decorator factory that memoizes all calls in a shared cache dict."""
     cache = dict() if cache is None else cache
 
-    def wrapper(func):
-        def wrapped(*args):
+    def wrapper(func: Callable[..., Any]) -> Callable[..., Any]:
+        def wrapped(*args: Any) -> Any:
             if args in cache:
                 return cache[args]
             result = func(*args)
@@ -456,14 +470,14 @@ def wrapped(*args):
 SOLVER_ORDER = [pulp.CPLEX, pulp.GLPK, pulp.PULP_CBC_CMD]
 
 
-def infer_additional_over_under_pairs(over_under_pairs: list, facet_triplets: set) -> list:
+def infer_additional_over_under_pairs(over_under_pairs: list[Any], facet_triplets: set[Any]) -> list[Any]:
     """Transitively close over/under relations: if A over B and B over C within a triplet, infer A over C."""
     # over_dict[f] is set of all facets that f lies over.
     over_dict = defaultdict(set)
     for over, under in over_under_pairs:
         over_dict[over].add(under)
 
-    def n_pairs():
+    def n_pairs() -> int:
         return sum(len(over_set) for over_set in over_dict.values())
 
     current = n_pairs()
@@ -487,13 +501,13 @@ def n_pairs():
 
 def find_folded_face_order(
     G: "EuclideanPositionHEG",
-    over_under_pairs=(),
-    solver=None,
+    over_under_pairs: Iterable[Any] = (),
+    solver: Any = None,
     double_fold_weight: float = 0,
     allow_slack: bool = True,
     problem_file: str | None = None,
     quiet: bool = False,
-) -> dict:
+) -> dict[Any, Any]:
     """Determine the stacking order of overlapping faces by solving an ILP.
 
     Assign 'sorted_original_faces' to each face of the overlap graph G and return
@@ -521,6 +535,7 @@ def find_folded_face_order(
         len(coincident_fold_lengths),
     )
 
+    over_under_pairs = list(over_under_pairs)
     n_over_under_before = len(over_under_pairs)
     over_under_pairs = infer_additional_over_under_pairs(over_under_pairs, set(triplet_overlap_areas.keys()))
 
@@ -532,15 +547,15 @@ def find_folded_face_order(
 
     n_vars = 0
 
-    def get_varname():
+    def get_varname() -> str:
         nonlocal n_vars
         n_vars += 1
         return "x" + str(n_vars)
 
-    over_dict = dict()
+    over_dict: dict[Any, Any] = dict()
 
     @cache_all(over_dict)
-    def over(face1, face2):
+    def over(face1: Any, face2: Any) -> Any:
         opposite = over_dict.get((face2, face1), None)
         if opposite is not None:
             return 1 - opposite
@@ -556,7 +571,7 @@ def over(face1, face2):
     prob += 0, ""  # empty objective function
     objective = 0
 
-    def add_constraint(constraint):
+    def add_constraint(constraint: Any) -> None:
         if isinstance(constraint, pulp.LpConstraint):
             nonlocal prob
             prob += constraint, ""
@@ -605,7 +620,7 @@ def add_constraint(constraint):
     else:
         logger.info("Skipping ILP since everything is already determined..")
 
-    def comparison_func(a, b):
+    def comparison_func(a: Any, b: Any) -> int:
         result = over_dict.get((a, b), 1 - over_dict.get((b, a), 0.5))
         if isinstance(result, (pulp.LpVariable, pulp.LpAffineExpression)):
             result = result.value()
@@ -618,17 +633,19 @@ def comparison_func(a, b):
 
     logger.info("determining face order..")
     for f in G.faces:
-        original_faces = list(f[ORIGINAL_FACES]).copy()
-        original_faces.sort(key=cmp_to_key(comparison_func))
-        f[SORTED_ORIGINAL_FACES] = original_faces
+        original_faces_list = list(f[ORIGINAL_FACES]).copy()
+        original_faces_list.sort(key=cmp_to_key(comparison_func))
+        f[SORTED_ORIGINAL_FACES] = original_faces_list
 
     logger.info("assigning creases..")
-    crease_assignment = dict()
-    original_faces = {f_orig for f in G.faces for f_orig in f[ORIGINAL_FACES]}
-    original_edges = {e for f in original_faces for e in f.halfedge_iter()}
+    crease_assignment: dict[Any, Any] = dict()
+    original_faces_set = {f_orig for f in G.faces for f_orig in f[ORIGINAL_FACES]}
+    original_edges = {e for f in original_faces_set for e in f.halfedge_iter()}
     for e in original_edges:
         if e.on_border() or e.rev.on_border():
             continue
+        assert e.face is not None
+        assert e.rev.face is not None
         crease_assignment[e] = comparison_func(e.face, e.rev.face) * (1 if e.face["color_key"] else -1)
 
     return crease_assignment
@@ -661,13 +678,16 @@ def get_over_under_pairs_from_creases(
     G: "EuclideanPositionHEG", two_coloring_key: str = "color_key"
 ) -> list[list["Face"]]:
     """Derive over/under face pairs from mountain/valley crease assignments on a two-colored graph."""
-    over_under_pairs = []
+    over_under_pairs: list[list[Face]] = []
     for e in G.halfedges:
         crease_type = e.attributes.get(CREASE_ASSIGNMENT, None)
         if crease_type in (MOUNTAIN, VALLEY) and not (e.on_border() or e.rev.on_border()):
+            assert e.face is not None
             e_above = e if e.face[two_coloring_key] else e.rev
             if crease_type is MOUNTAIN:
                 e_above = e_above.rev
+            assert e_above.face is not None
+            assert e_above.rev.face is not None
             over_under_pairs.append([e_above.face, e_above.rev.face])
     logger.info("number of pairs: %d", len(over_under_pairs))
     return over_under_pairs
@@ -680,8 +700,8 @@ def get_over_under_pairs_from_creases(
 def face_order_to_clean_graph(
     G: "EuclideanPositionHEG",
     side: str = TOP,
-    top_color: tuple = (0.5, 0.5, 0.9),
-    bottom_color: tuple = (0.8, 0.8, 0.8),
+    top_color: tuple[float, ...] = (0.5, 0.5, 0.9),
+    bottom_color: tuple[float, ...] = (0.8, 0.8, 0.8),
 ) -> "EuclideanPositionHEG":
     """Extract a clean renderable graph showing only the top or bottom layer of a folded model."""
     view = G.copy()
@@ -701,12 +721,14 @@ def face_order_to_clean_graph(
             key = not key
         f["color_key"] = color_key_mapping[key]
 
-    to_delete = [
-        e
-        for e in view.halfedges
-        if not (e.on_border() or e.rev.on_border())
-        and e.face["sorted_original_faces"][layer] is e.rev.face["sorted_original_faces"][layer]
-    ]
+    to_delete = []
+    for e in view.halfedges:
+        if e.on_border() or e.rev.on_border():
+            continue
+        assert e.face is not None
+        assert e.rev.face is not None
+        if e.face["sorted_original_faces"][layer] is e.rev.face["sorted_original_faces"][layer]:
+            to_delete.append(e)
     view.delete_subset(to_delete)
     view.recompute_lengths_and_angles()
 
@@ -781,7 +803,7 @@ def fold_complete(
 
 
 @contextmanager
-def _quiet_progress(quiet: bool):
+def _quiet_progress(quiet: bool) -> Iterator[None]:
     """Context manager: while ``quiet`` is True, silence this module's tqdm bars and logger."""
     if not quiet:
         yield
@@ -790,11 +812,11 @@ def _quiet_progress(quiet: bool):
     module = sys.modules[__name__]
     original_tqdm = module.tqdm
 
-    def silent_tqdm(iterable=None, *args, **kwargs):
+    def silent_tqdm(iterable: Any = None, *args: Any, **kwargs: Any) -> Any:
         kwargs["disable"] = True
         return original_tqdm(iterable, *args, **kwargs)
 
-    module.tqdm = silent_tqdm
+    module.tqdm = silent_tqdm  # type: ignore[attr-defined]
     previous_level = logger.level
     previous_disabled = logger.disabled
     logger.setLevel(logging.CRITICAL + 1)
@@ -802,7 +824,7 @@ def silent_tqdm(iterable=None, *args, **kwargs):
     try:
         yield
     finally:
-        module.tqdm = original_tqdm
+        module.tqdm = original_tqdm  # type: ignore[attr-defined]
         logger.setLevel(previous_level)
         logger.disabled = previous_disabled
 
@@ -836,12 +858,12 @@ def __init__(
             if value is not None:
                 self[name] = value
 
-    def __getattr__(self, name: str):
+    def __getattr__(self, name: str) -> Any:
         if name in self._FIELDS:
             return self.get(name)
         raise AttributeError(name)
 
-    def __setattr__(self, name: str, value) -> None:
+    def __setattr__(self, name: str, value: Any) -> None:
         if name in self._FIELDS:
             if value is None:
                 self.pop(name, None)
@@ -882,9 +904,9 @@ def show(
         if render_settings.get("line_width", "auto") == "auto" and self.CP is not None:
             render_settings["line_width"] = CairoRenderer.auto_line_width(self.CP)
 
-        graphs: list = []
-        titles: list[str] = []
-        per_subplot: list[dict] = []
+        graphs: list[Any] = []
+        titles: list[str | None] = []
+        per_subplot: list[dict[str, Any]] = []
 
         if self.CP is not None:
             graphs.append(self.CP)
@@ -934,20 +956,20 @@ def show(
             **render_settings,
         )
 
-    def save(self, path: str, **save_results_kwargs):
+    def save(self, path: str, **save_results_kwargs: Any) -> None:
         """Convenience wrapper around :func:`save_results`."""
         save_results(self, path=path, **save_results_kwargs)
 
 
 def save_results(
-    results,
+    results: Any,
     path: str = "results",
-    render_settings: dict | None = None,
+    render_settings: dict[str, Any] | None = None,
     min_foldable_length: float | None = None,
     bbox: tuple[float, float] | None = None,
-    extra_info: str = None,
+    extra_info: str | None = None,
     opacity: float = 0.15,
-):
+) -> None:
     """Save rendered crease pattern, folded views, and a plotter-ready SVG to a directory.
 
     The plotter SVG is scaled to fit within bbox (in cm) or so that the shortest
@@ -1040,7 +1062,9 @@ def save_results(
     logger.info(text)
 
 
-def remove_duplicates(G: "EuclideanPositionHEG", eps: float = 1e-6, exclude_edges=()) -> "EuclideanPositionHEG":
+def remove_duplicates(
+    G: "EuclideanPositionHEG", eps: float = 1e-6, exclude_edges: Iterable["HalfEdge"] = ()
+) -> "EuclideanPositionHEG":
     """Merge duplicate vertices and edges within eps distance, returning a clean graph.
 
     The faces of the resulting graph are built from scratch based on the new edges, so this only works for planar graphs.
@@ -1053,15 +1077,13 @@ def remove_duplicates(G: "EuclideanPositionHEG", eps: float = 1e-6, exclude_edge
     node_mapping = {i: j for i, j in enumerate(index[inverse])}
     v_index = {v: node_mapping[i] for i, v in enumerate(vs)}
 
-    nxG = nx.Graph()
+    nxG: nx.Graph = nx.Graph()
     nxG.add_nodes_from(v_index.values())
     nx_positions = {i: pos[i] for i in node_mapping.values()}
+    exclude_set = set(exclude_edges)
+    exclude_set_with_rev = exclude_set.union({e.rev for e in exclude_set})
     nxG.add_edges_from(
-        [
-            (v_index[e.orig], v_index[e.dest])
-            for e in G.halfedges_representing_edges()
-            if e not in set(exclude_edges).union({e.rev for e in exclude_edges})
-        ]
+        [(v_index[e.orig], v_index[e.dest]) for e in G.halfedges_representing_edges() if e not in exclude_set_with_rev]
     )
 
     G2 = EHEG_from_nx(nxG, nx_positions)
diff --git a/eucare/plotting.py b/eucare/plotting.py
index b3a3b22..66e9fff 100755
--- a/eucare/plotting.py
+++ b/eucare/plotting.py
@@ -2,6 +2,8 @@
 
 from __future__ import annotations
 
+from typing import Any
+
 import matplotlib.collections as mc
 import matplotlib.pyplot as plt
 import numpy as np
@@ -9,17 +11,17 @@
 from numpy.typing import NDArray
 
 
-def plot_lines(lines: NDArray, ax: Axes | None = None, **kwargs) -> None:
+def plot_lines(lines: NDArray, ax: Axes | None = None, **kwargs: Any) -> None:
     """Add a :class:`LineCollection` of *lines* (shape ``(n, 2, 2)``) to *ax*.
 
     Extra keyword arguments are forwarded to :class:`matplotlib.collections.LineCollection`.
     """
-    lc = mc.LineCollection(lines, **kwargs)
+    lc = mc.LineCollection(list(lines), **kwargs)
     ax = plt.gca() if ax is None else ax
     ax.add_collection(lc)
 
 
-def plot_polygon(points: NDArray, **kwargs) -> None:
+def plot_polygon(points: NDArray, **kwargs: Any) -> None:
     """Plot the closed polygon defined by *points* (shape ``(n, 2)``)."""
     plot_lines(
         np.stack(
diff --git a/eucare/prototiles.py b/eucare/prototiles.py
index da11619..22f84cb 100755
--- a/eucare/prototiles.py
+++ b/eucare/prototiles.py
@@ -2,18 +2,25 @@
 
 from __future__ import annotations
 
+from typing import Any, Callable
+
 import numpy as np
+from numpy.typing import NDArray
 
 from .base import angle_to_axis, unit_vector
 from .geometries import EuclideanGeometry, PoincareDiskModel, SphereModel
-from .half import CyclicHalfedgeGraph, Face, HalfEdge, InAngleHEG, Vertex
+from .geometries.base import Geometry
+from .half import CyclicHalfedgeGraph, Face, HalfEdge, HalfEdgeGraph, InAngleHEG, Vertex
 from .instructions import attatch_tile_instruction
 
 
 class ProtoTile:
     """Abstract base class for tile prototypes that can produce half-edge graphs."""
 
-    def make_graph(self) -> tuple[InAngleHEG, dict]:
+    geometry: type[Geometry]
+    points: NDArray[Any]
+
+    def make_graph(self, add_positions: bool = False) -> tuple[CyclicHalfedgeGraph, dict[Any, HalfEdge]]:
         # return a HEG and a list of edges
         raise NotImplementedError
 
@@ -23,12 +30,12 @@ class PolygonalProtoTile(ProtoTile):
 
     def __init__(
         self,
-        in_angles: list[float],
-        edge_lengths: list[float],
-        edge_labels: list | None = None,
-        vertex_labels: list | None = None,
+        in_angles: list[float] | NDArray[Any],
+        edge_lengths: list[float] | NDArray[Any],
+        edge_labels: list[Any] | None = None,
+        vertex_labels: list[Any] | None = None,
         face_label: object = None,
-        edge_instructions: dict | None = None,
+        edge_instructions: dict[Any, Callable[[HalfEdgeGraph, HalfEdge], None]] | None = None,
     ) -> None:
         assert len(in_angles) == len(edge_lengths)
         self.order = len(in_angles)
@@ -41,9 +48,9 @@ def __init__(
         # edge_instructions can either be a list for all edges, or a dict, mapping a label to an instruction
         self.edge_instructions = edge_instructions if edge_instructions is not None else dict()
 
-    def make_graph(self, add_positions: bool = False) -> tuple[CyclicHalfedgeGraph, dict]:
+    def make_graph(self, add_positions: bool = False) -> tuple[CyclicHalfedgeGraph, dict[Any, HalfEdge]]:
         """Build a half-edge graph for this tile and return it with an edge label dict."""
-        outer_edge_dict = dict()
+        outer_edge_dict: dict[Any, HalfEdge] = dict()
 
         outer_edges = [HalfEdge() for _ in range(self.order)]
         for e, label in zip(outer_edges, self.edge_labels):
@@ -78,12 +85,12 @@ def make_graph(self, add_positions: bool = False) -> tuple[CyclicHalfedgeGraph,
         graph = CyclicHalfedgeGraph(f=f, vs=vertices, inner_hs=inner_edges, outer_hs=outer_edges)
         return graph, outer_edge_dict
 
-    def attach_instruction(self, label=None):
+    def attach_instruction(self, label: object = None) -> Callable[[HalfEdgeGraph, HalfEdge], None]:
         """Return a glue instruction that attaches this tile at the given edge label."""
         assert label is None or label in self.edge_labels, f"{label}, {self.edge_labels}"
         return attatch_tile_instruction(self, label)
 
-    def __str__(self):
+    def __str__(self) -> str:
         return (
             f"PolygonalProtoTile("
             f"lenghts={self.edge_lengths}, "
@@ -97,7 +104,7 @@ def __str__(self):
 class RegularProtoTile(PolygonalProtoTile):
     """A regular polygon tile with uniform angles and edge lengths, supporting any geometry."""
 
-    def __init__(self, n: int, in_angle: float, edge_length: float, **super_kwargs) -> None:
+    def __init__(self, n: int, in_angle: float, edge_length: float, **super_kwargs: Any) -> None:
         super().__init__([in_angle] * n, [edge_length] * n, **super_kwargs)
         euclidean_angle_deficit = np.pi * (n - 2) - in_angle * n
         eps = 1e-6
@@ -108,14 +115,14 @@ def __init__(self, n: int, in_angle: float, edge_length: float, **super_kwargs)
         else:
             self.geometry = PoincareDiskModel
 
-        geo = self.geometry
+        geo: type[Geometry] = self.geometry
         # calculate points
-        points = [geo.origin(), geo.from_polar(edge_length, 0)]
+        points_list = [geo.origin(), geo.from_polar(edge_length, 0)]
         for _ in range(n - 2):
-            points.append(geo.construct_next_poly_point(points[-2], points[-1], in_angle, edge_length))
-        com = geo.center_of_mass(np.array(points))
+            points_list.append(geo.construct_next_poly_point(points_list[-2], points_list[-1], in_angle, edge_length))
+        com = geo.center_of_mass(np.array(points_list))
         translate = geo.translation(com, geo.origin())
-        points = np.array([translate(point) for point in points])
+        points = np.array([translate(point) for point in points_list])
         angle = -geo.angle_to_axis(points[0])
         rotate = geo.rotation(geo.origin(), angle)
         self.points = np.array([rotate(point) for point in points])
@@ -124,7 +131,7 @@ def __init__(self, n: int, in_angle: float, edge_length: float, **super_kwargs)
 class EuclideanProtoTile(PolygonalProtoTile):
     """A Euclidean tile defined by explicit 2D vertex positions."""
 
-    def __init__(self, points=None, **super_kwargs):
+    def __init__(self, points: NDArray[Any] | list[Any] | None = None, **super_kwargs: Any) -> None:
         self.geometry = EuclideanGeometry
         # points should have shape (n, 2)
         points = np.array(points)
@@ -141,7 +148,7 @@ def __init__(self, points=None, **super_kwargs):
 class RegularEuclideanTile(EuclideanProtoTile):
     """A regular n-gon in Euclidean geometry with unit edge length."""
 
-    def __init__(self, n, **super_kwargs):
+    def __init__(self, n: int, **super_kwargs: Any) -> None:
         points = unit_vector(np.linspace(0, 2 * np.pi, n, endpoint=False) + np.pi / n) / np.sin(np.pi / n) / 2
         super(RegularEuclideanTile, self).__init__(points=points, **super_kwargs)
 
@@ -149,7 +156,7 @@ def __init__(self, n, **super_kwargs):
 class RhombusTile(EuclideanProtoTile):
     """A rhombus tile with a given acute angle alpha (default pi/3)."""
 
-    def __init__(self, alpha=None, **super_kwargs):
+    def __init__(self, alpha: float | None = None, **super_kwargs: Any) -> None:
         alpha = np.pi / 3 if alpha is None else alpha
         pts = np.concatenate(
             [np.zeros((1, 2)), np.cumsum(unit_vector([-alpha / 2, alpha / 2, np.pi - alpha / 2]), axis=0)]
@@ -157,7 +164,7 @@ def __init__(self, alpha=None, **super_kwargs):
         super(RhombusTile, self).__init__(points=pts, **super_kwargs)
 
 
-def complete_vertex_with_rhombus(graph, vertex):
+def complete_vertex_with_rhombus(graph: InAngleHEG, vertex: Vertex) -> None:
     """Fill the remaining angle at a border vertex by gluing in a rhombus tile."""
     assert isinstance(graph, InAngleHEG)
     edge = vertex.get_outgoing_border()
diff --git a/eucare/rendering.py b/eucare/rendering.py
index 0883c90..d4fcfcc 100755
--- a/eucare/rendering.py
+++ b/eucare/rendering.py
@@ -10,16 +10,17 @@
 import logging
 import os
 from io import BytesIO
-from typing import TYPE_CHECKING, Iterable
+from typing import TYPE_CHECKING, Any, Iterable
 
 import cairo
 import numpy as np
+from numpy.typing import NDArray
 
 from .base import angle_to_axis, unit_vector
 from .half import rotate_by
 
 if TYPE_CHECKING:
-    from .half import HalfEdgeGraph
+    from .half import Face, HalfEdge, HalfEdgeGraph, Vertex
 
 logger = logging.getLogger(__name__)
 try:
@@ -174,7 +175,7 @@ def inset_poly(pts: list, dist: float) -> list:
 _seed_offset = np.random.randint(2**16)
 
 
-def random_color(seed: object = None) -> np.ndarray:
+def random_color(seed: object = None) -> NDArray:
     """Return a random RGB triple in ``[0, 1]^3``, optionally seeded by a hashable."""
     if seed is not None:
         np.random.seed((hash(seed) + _seed_offset) % 2**32)
@@ -191,18 +192,21 @@ def is_color(obj: object) -> bool:
             except ValueError:
                 return False
         return False
-    return isinstance(obj, Iterable) and len(obj) in (3, 4) and all([isinstance(c, (int, float)) for c in obj])
+    if not isinstance(obj, Iterable):
+        return False
+    seq = list(obj)
+    return len(seq) in (3, 4) and all(isinstance(c, (int, float)) for c in seq)
 
 
 def multi_show(
-    graphs: Iterable,
-    titles: list[str] | None = None,
+    graphs: Iterable[Any],
+    titles: list[str | None] | None = None,
     ncols: int | None = None,
     figsize: tuple[float, float] | None = None,
     suptitle: str | None = None,
     cell_size: float = 4.0,
-    per_subplot_kwargs: list[dict] | None = None,
-    **show_kwargs: object,
+    per_subplot_kwargs: list[dict[str, Any]] | None = None,
+    **show_kwargs: Any,
 ) -> None:
     """Render multiple half-edge graphs side-by-side in a matplotlib grid.
 
@@ -282,6 +286,9 @@ class CairoRenderer:
     back to a translucent blue (faces) or grey (edges).
     """
 
+    surface: cairo.SVGSurface
+    dc: cairo.Context[Any]
+
     def __init__(
         self,
         width: int | None = None,
@@ -308,20 +315,22 @@ def __init__(
             curve_position_key: Half-edge attribute holding curved-fold polylines
                 (overrides the straight ``orig -> dest`` line).
         """
-        if width is None and height is None:
-            width, height = 512, 512
-        self.width = width if width is not None else height
-        self.height = height if height is not None else width
-        self.scale = scale
-        self.line_width = line_width
+        if width is None:
+            width = height if height is not None else 512
+        if height is None:
+            height = width
+        self.width = width
+        self.height = height
+        self.scale: float | str = scale
+        self.line_width: float | str = line_width
         self.vertex_radius = vertex_radius
         self.face_inset = face_inset
         self.position_key = position_key
         self.curve_position_key = curve_position_key
-        self.surface = None
-        self.dc = None
+        # ``surface`` and ``dc`` are set by ``render_graph``; reading them
+        # before that call raises ``AttributeError``.
 
-    def render_face(self, face, color_key: str = "color_key"):
+    def render_face(self, face: Face, color_key: str = "color_key") -> cairo.SVGSurface:
         """Draw a single face's filled polygon, honouring ``face[color_key]``."""
         dc = self.dc
         points = []
@@ -331,7 +340,7 @@ def render_face(self, face, color_key: str = "color_key"):
             else:
                 points.extend(h[self.curve_position_key][:-1])
         inset = self.face_inset
-        if inset != 0:
+        if inset is not None and inset != 0:
             points = inset_poly(points, inset)
         dc.move_to(*points[-1])
         for point in points:
@@ -352,7 +361,7 @@ def render_face(self, face, color_key: str = "color_key"):
         dc.stroke()
         return self.surface
 
-    def set_source_color(self, color) -> None:
+    def set_source_color(self, color: Any) -> None:
         """Set the current cairo source colour, accepting RGB/RGBA tuples or any hashable seed."""
         if not is_color(color):
             color = random_color(color)
@@ -363,7 +372,13 @@ def set_source_color(self, color) -> None:
         else:
             self.dc.set_source_rgba(*color)
 
-    def render_edge(self, edge, color_key: str = "color_key", last_pos=None, tol: float = 1e-6):
+    def render_edge(
+        self,
+        edge: HalfEdge,
+        color_key: str = "color_key",
+        last_pos: NDArray | None = None,
+        tol: float = 1e-6,
+    ) -> cairo.SVGSurface | tuple[cairo.SVGSurface, NDArray]:
         """Draw a single half-edge as a line or a tapered ribbon.
 
         If both endpoints carry a ``line_width`` attribute, the edge is drawn
@@ -411,6 +426,7 @@ def render_edge(self, edge, color_key: str = "color_key", last_pos=None, tol: fl
                 self.set_source_color((0.5, 0.5, 1.0))
 
             if edge.attributes.get("delete", False):
+                assert isinstance(self.line_width, (int, float))
                 dc.set_dash([self.line_width * 2, self.line_width * 3])
                 dc.stroke()
                 dc.set_dash([])
@@ -419,14 +435,15 @@ def render_edge(self, edge, color_key: str = "color_key", last_pos=None, tol: fl
                 # dc.stroke()
         return self.surface if last_pos is None else (self.surface, edge.dest[self.position_key])
 
-    def render_vertex(self, vertex, color_key: str = "color_key") -> None:
+    def render_vertex(self, vertex: Vertex, color_key: str = "color_key") -> None:
         """Draw a vertex as a small disc, picking colour from ``vertex[color_key]`` or attribute flags."""
         dc = self.dc
         dc.set_line_width(0)
         radius = vertex.get("vertex_radius", self.vertex_radius)
         if radius == 0:
             return
-        dc.arc(*vertex[self.position_key], vertex.get("vertex_radius", self.vertex_radius), 0, 2 * np.pi)
+        pos = vertex[self.position_key]
+        dc.arc(pos[0], pos[1], vertex.get("vertex_radius", self.vertex_radius), 0, 2 * np.pi)
         if color_key in vertex.attributes:
             self.set_source_color(vertex[color_key])
         elif vertex.attributes.get("join", False):
@@ -438,9 +455,10 @@ def render_vertex(self, vertex, color_key: str = "color_key") -> None:
         dc.fill_preserve()
         dc.set_source_rgb(0.0, 0.0, 0.0)
         dc.stroke()
+        assert isinstance(self.line_width, (int, float))
         dc.set_line_width(self.line_width)
 
-    def autoscale(self, graph):
+    def autoscale(self, graph: HalfEdgeGraph) -> CairoRenderer:
         """Pick an isotropic scale to fit *graph* in the surface (origin-centred)."""
         positions = np.array([v[self.position_key] for v in graph.vertices])
         max_abs_pos = np.max(np.abs(positions), axis=0)
@@ -448,7 +466,7 @@ def autoscale(self, graph):
         self.dc.scale(scale, scale)
         return self
 
-    def autocenterscale(self, graph):
+    def autocenterscale(self, graph: HalfEdgeGraph) -> CairoRenderer:
         """Pick scale and translation to fit *graph* in the surface, centred on its bounding box."""
         positions = np.array([v[self.position_key] for v in graph.vertices])
         offset = (np.max(positions, axis=0) + np.min(positions, axis=0)) / 2
@@ -460,7 +478,7 @@ def autocenterscale(self, graph):
         return self
 
     @staticmethod
-    def auto_line_width(graph) -> float:
+    def auto_line_width(graph: HalfEdgeGraph) -> float:
         """Default line width for *graph*: ``min(min_edge / 2, mean_edge / 10)``."""
         lengths = np.array(
             [
@@ -472,12 +490,12 @@ def auto_line_width(graph) -> float:
 
     def render_graph(
         self,
-        graph: "HalfEdgeGraph",
+        graph: HalfEdgeGraph,
         render_vertices: bool = True,
         render_faces: bool = True,
         render_edges: bool = True,
         for_cutting: bool = False,
-    ) -> "Rendering":
+    ) -> Rendering:
         """Render the whole *graph* (faces, edges, vertices) onto the cairo surface.
 
         Args:
@@ -506,6 +524,7 @@ def render_graph(
         if self.scale == "auto":
             self.autocenterscale(graph)
         else:
+            assert isinstance(self.scale, (int, float))
             self.dc.scale(self.scale, self.scale)
         # self.dc.set_font_size(18.0 / self.scale)
 
@@ -513,8 +532,11 @@ def render_graph(
             self.line_width = self.auto_line_width(graph)
         elif isinstance(self.line_width, str) and self.line_width.endswith("%"):
             self.line_width = float(self.line_width[:-1]) / 100 * self.auto_line_width(graph)
-        self.vertex_radius = self.vertex_radius if self.vertex_radius is not None else self.line_width
-        self.face_inset = self.face_inset if self.face_inset is not None else self.line_width
+        assert isinstance(self.line_width, (int, float))
+        if self.vertex_radius is None:
+            self.vertex_radius = self.line_width
+        if self.face_inset is None:
+            self.face_inset = self.line_width
 
         if render_faces:
             for f in graph.faces:
@@ -582,7 +604,7 @@ def __init__(self, position_key: str = "pos", curve_position_key: str = "curve_p
         self.position_key = position_key
         self.curve_position_key = curve_position_key
 
-    def _render_halfedges(self, halfedges, dwg, bbox, scale):
+    def _render_halfedges(self, halfedges: Iterable[HalfEdge], dwg: Any, bbox: NDArray, scale: float) -> None:
         # TODO: adapt this to get an even better path: https://stackoverflow.com/a/44080908
         edges = list(halfedges)
         edge_to_index = {e: i for i, e in enumerate(edges)}
@@ -613,14 +635,14 @@ def _render_halfedges(self, halfedges, dwg, bbox, scale):
         # print(f'Rendering {len(polylines)} polylines, with {[len(line) for line in polylines]} line segments.')
 
         for pts in polylines:
-            pts = np.array(pts)
-            pts -= bbox[:, 0][None]
-            pts *= scale
+            pts_arr = np.array(pts)
+            pts_arr -= bbox[:, 0][None]
+            pts_arr *= scale
             pth = dwg.path(fill_opacity=0, stroke_width="0.05", stroke="black")
-            pth.push("M", *pts)
+            pth.push("M", *pts_arr)
             dwg.add(pth)
 
-    def create_drawing(self, filename: str, pts: np.ndarray, height: float, unit=svgwrite.cm) -> None:
+    def create_drawing(self, filename: str, pts: NDArray, height: float, unit: Any = svgwrite.cm) -> None:
         """Initialise the underlying ``svgwrite.Drawing`` based on the bounding box of *pts*."""
         self.bbox = np.array([[f(pts[:, i]) for f in [np.min, np.max]] for i in [0, 1]])
         aspect_ratio = (self.bbox[0, 1] - self.bbox[0, 0]) / (self.bbox[1, 1] - self.bbox[1, 0])
@@ -634,13 +656,13 @@ def create_drawing(self, filename: str, pts: np.ndarray, height: float, unit=svg
     def render_graph(
         self,
         filename: str,
-        graph,
+        graph: Any,
         render_vertices: bool = False,
         render_faces: bool = False,
         render_edges: bool = True,
         for_cutting: bool = True,
         height: float = 30,
-        unit=svgwrite.cm,
+        unit: Any = svgwrite.cm,
         render_interior_and_borders: bool = True,
         extra_render_keys: tuple[str, ...] = ("drawing_edge",),
     ) -> None:
diff --git a/eucare/shrink_rotate/__init__.py b/eucare/shrink_rotate/__init__.py
index 24fc6fe..f140b5f 100644
--- a/eucare/shrink_rotate/__init__.py
+++ b/eucare/shrink_rotate/__init__.py
@@ -28,6 +28,8 @@
 
 from __future__ import annotations
 
+from typing import Any
+
 from .crease_orientation import (
     THIS_WAY,
     assign_this_way_by_distance,
@@ -61,7 +63,7 @@
 ]
 
 
-def __getattr__(name: str):
+def __getattr__(name: str) -> Any:
     """Lazily expose :class:`ShrinkRotateExplorer` (avoids importing matplotlib eagerly)."""
     if name == "ShrinkRotateExplorer":
         from .widgets import ShrinkRotateExplorer
diff --git a/eucare/shrink_rotate/crease_orientation.py b/eucare/shrink_rotate/crease_orientation.py
index 6f188e1..2c8b30c 100644
--- a/eucare/shrink_rotate/crease_orientation.py
+++ b/eucare/shrink_rotate/crease_orientation.py
@@ -50,6 +50,7 @@
 from typing import Iterable
 
 import numpy as np
+from numpy.typing import NDArray
 
 from .. import base
 from ..half import Face, HalfEdgeGraph, Vertex
@@ -145,10 +146,13 @@ def assign_this_way_by_bfs(G: HalfEdgeGraph, source: Vertex | Face | set[Vertex]
 
     Skips edges that already have THIS_WAY assigned on either side.
     """
+    source_set: set[Vertex | Face]
     if not isinstance(source, set):
-        source = {source}
-    source_faces = {f for s in source for f in (s.true_face_iter() if isinstance(s, Vertex) else [s])}
-    source_vertices = {v for s in source for v in (s.vertex_iter() if isinstance(s, Face) else [s])}
+        source_set = {source}
+    else:
+        source_set = set(source)
+    source_faces: set[Face] = {f for s in source_set for f in (s.true_face_iter() if isinstance(s, Vertex) else [s])}
+    source_vertices: set[Vertex] = {v for s in source_set for v in (s.vertex_iter() if isinstance(s, Face) else [s])}
 
     assign_this_way_by_face_bfs(G, source_faces)  # primary: face BFS
     assign_this_way_by_vertex_bfs(G, source_vertices)  # tiebreaker: vertex BFS
@@ -221,7 +225,7 @@ def assign_this_way_from_center(G: HalfEdgeGraph) -> None:
     assign_this_way_by_bfs(G, src)
 
 
-def assign_this_way_by_distance(G: HalfEdgeGraph, point=None) -> None:
+def assign_this_way_by_distance(G: HalfEdgeGraph, point: NDArray[np.floating] | None = None) -> None:
     """Orient interior edges by distance of face midpoints from *point*.
 
     The face whose midpoint is *farther* from *point* lies below (gets
diff --git a/eucare/shrink_rotate/pipeline.py b/eucare/shrink_rotate/pipeline.py
index 0cf1011..549ee6c 100644
--- a/eucare/shrink_rotate/pipeline.py
+++ b/eucare/shrink_rotate/pipeline.py
@@ -15,13 +15,14 @@
 from __future__ import annotations
 
 import logging
+from typing import Any
 
 import numpy as np
 
 from ..base import rotation_matrix
 from ..conway import shrink_rotate_graph
 from ..flat_foldable import max_kawasaki_sum
-from ..half import EuclideanPositionHEG, GeometricHEG, HalfEdgeGraph
+from ..half import EuclideanPositionHEG, Face, GeometricHEG, HalfEdge, HalfEdgeGraph, Vertex
 from ..overlap import BORDER, CREASE_ASSIGNMENT, MOUNTAIN, VALLEY
 from ..rendering import BORDER_COLOR, MOUNTAIN_COLOR, VALLEY_COLOR
 from ..utils import invert_mapping
@@ -38,8 +39,8 @@ def shrink_rotate_pattern(
     *,
     assign_creases: bool = True,
     simplify_boundary: bool = True,
-    **reciprocal_figure_kwargs,
-) -> EuclideanPositionHEG:
+    **reciprocal_figure_kwargs: Any,
+) -> GeometricHEG:
     """Build a shrink-rotate crease pattern from tiling *G*.
 
     Each face of *G* is subdivided by the shrink-rotate Conway operator;
@@ -155,9 +156,11 @@ def assign_shrink_rotate_creases(SRG: HalfEdgeGraph) -> None:
         while e_twist.rev.on_border():
             e_twist = e_twist.nex
         # find the original-graph edge corresponding to e_twist.
-        e = None
+        e: HalfEdge | None = None
+        neighbor_face = e_twist.rev.nex.nex.rev.face
+        assert neighbor_face is not None
         for e_orig in f["pre_conway"].halfedge_iter():
-            if e_orig.rev in e_twist.rev.nex.nex.rev.face["pre_conway"].halfedge_iter():
+            if e_orig.rev in neighbor_face["pre_conway"].halfedge_iter():
                 e = e_orig
                 break
         assert e is not None
diff --git a/eucare/shrink_rotate/reciprocal_figures.py b/eucare/shrink_rotate/reciprocal_figures.py
index a059031..db406c8 100644
--- a/eucare/shrink_rotate/reciprocal_figures.py
+++ b/eucare/shrink_rotate/reciprocal_figures.py
@@ -16,12 +16,14 @@
 import logging
 from copy import copy
 
+
 import numpy as np
+from numpy.typing import NDArray
 import scipy as sc
 
 from ..base import rotation_matrix
 from ..conway import dual_graph
-from ..half import GeometricHEG
+from ..half import Face, GeometricHEG, HalfEdge, HalfEdgeGraph, Vertex
 from ..utils import invert_mapping, random_directed_set
 
 logger = logging.getLogger(__name__)
@@ -31,7 +33,7 @@ def reciprocal_figure(
     G: GeometricHEG,
     reciprocal_pos_key: str = "reciprocal_pos",
     rcond: float = 1e-7,
-):
+) -> GeometricHEG:
     """Compute the reciprocal figure of *G* and return it as a face graph.
 
     Stores reciprocal positions on faces of *G* under
@@ -73,13 +75,15 @@ def reciprocal_figure(
     # vertices approximate primal face centroids.
     to_process = set(G.faces)
     anchor = to_process.pop()
-    coefficients = {anchor: np.zeros(n_edges, dtype=np.float32)}
-    border = {anchor}
+    coefficients: dict[Face, NDArray[np.float32]] = {anchor: np.zeros(n_edges, dtype=np.float32)}
+    border: set[Face] = {anchor}
     while border:
-        new_border = set()
+        new_border: set[Face] = set()
         for f in border:
             for e in f.halfedge_iter():
                 f2 = e.rev.face
+                if f2 is None:
+                    continue
                 if f2 not in coefficients:
                     if e in directed_edges:
                         coefficients[f2] = copy(coefficients[f])
diff --git a/eucare/shrink_rotate/widgets.py b/eucare/shrink_rotate/widgets.py
index 324f120..30076fc 100644
--- a/eucare/shrink_rotate/widgets.py
+++ b/eucare/shrink_rotate/widgets.py
@@ -6,13 +6,17 @@
 
 from __future__ import annotations
 
+from typing import Any
+
 import numpy as np
 
+from ..half import HalfEdgeGraph
+
 from .. import base
 from ..utils import random_directed_set
 
 
-def _require_notebook_deps():
+def _require_notebook_deps() -> None:
     try:
         import ipywidgets as widgets  # noqa: F401
         import matplotlib.pyplot as plt  # noqa: F401
@@ -56,7 +60,7 @@ class ShrinkRotateExplorer:
 
     def __init__(
         self,
-        SRG,
+        SRG: HalfEdgeGraph,
         *,
         alpha0: float = 1 / 6,
         factor0: float = 0.58,
@@ -79,8 +83,8 @@ def __init__(
         # the plot reliably (and exactly once).
         with plt.ioff():
             fig = plt.figure(num=figure_id, figsize=figsize, clear=True)
-        fig.canvas.header_visible = False
-        fig.canvas.footer_visible = False
+        fig.canvas.header_visible = False  # type: ignore[attr-defined]
+        fig.canvas.footer_visible = False  # type: ignore[attr-defined]
         ax = fig.add_subplot(1, 1, 1)
         lc = LineCollection(self._segments(), antialiased=True, color="k", linewidth=1)
         pc = PolyCollection(self._polys(), antialiased=True, color="k", alpha=0.1)
@@ -131,7 +135,7 @@ def __init__(
     # ------------------------------------------------------------------ state
 
     @staticmethod
-    def _build_state(SRG) -> dict:
+    def _build_state(SRG: HalfEdgeGraph) -> dict[str, Any]:
         """Precompute index arrays so each tick is one batched matmul."""
         vertex_list = list(SRG.vertices)
         vidx = {id(v): i for i, v in enumerate(vertex_list)}
@@ -171,18 +175,27 @@ def _reshrinkrotate(self, alpha: float, factor: float, global_scale: float = 1.0
         for i, v in enumerate(s["vertex_list"]):
             v["pos"] = pos[i]
 
-    def _segments(self):
+    def _segments(self) -> Any:
         s = self._state
         return s["positions"][s["edge_idx"]]
 
-    def _polys(self):
+    def _polys(self) -> list[Any]:
         s = self._state
         pos = s["positions"]
         return [pos[idx] for idx in s["face_idx"]]
 
     # ----------------------------------------------------------------- update
 
-    def _update(self, alpha, factor, folded, reparametrized, scale_folded, show_lines, show_polys):
+    def _update(
+        self,
+        alpha: float,
+        factor: float,
+        folded: bool,
+        reparametrized: bool,
+        scale_folded: bool,
+        show_lines: bool,
+        show_polys: bool,
+    ) -> None:
         if self._in_update:
             return
         if factor <= 0:
@@ -238,7 +251,7 @@ def _update(self, alpha, factor, folded, reparametrized, scale_folded, show_line
 
     # ----------------------------------------------------------------- public
 
-    def display(self):
+    def display(self) -> None:
         """Display the widget UI in a Jupyter notebook.
 
         Re-displays the figure canvas explicitly, so re-executing a cell
diff --git a/eucare/svg.py b/eucare/svg.py
index e4a83cf..84d3cf3 100644
--- a/eucare/svg.py
+++ b/eucare/svg.py
@@ -6,10 +6,12 @@
 import operator
 import re
 from collections import defaultdict
+from typing import Any
 
 import networkx as nx
 import numpy as np
 import svgpathtools as spt
+from numpy.typing import NDArray
 
 import eucare as ec
 from eucare.overlap import CREASE_ASSIGNMENT, MOUNTAIN, VALLEY
@@ -37,28 +39,29 @@ def load_svg(filepath: str) -> ec.half.EuclideanPositionHEG:
         crease assignments stored on each half-edge under :data:`CREASE_ASSIGNMENT`.
     """
 
-    def get_stroke(attrs: dict) -> str | None:
+    def get_stroke(attrs: dict[str, Any]) -> str | None:
         """Extract the ``stroke:...;`` colour from an SVG ``style`` attribute string."""
         # hits = re.search('stroke:(#.{6})', attrs['style'])
-        hits = re.search("stroke:(.*?);", attrs["style"])
-        if hits is not None:
-            hits = hits.groups()
-        if not hits:
+        match = re.search("stroke:(.*?);", attrs["style"])
+        if match is None:
             return None
-        elif len(hits) == 1:
-            return hits[0]
+        groups = match.groups()
+        if not groups:
+            return None
+        elif len(groups) == 1:
+            return groups[0]
         else:
-            raise ValueError(f"Found multiple strokes: {list(hits)}")
+            raise ValueError(f"Found multiple strokes: {list(groups)}")
 
     paths, attributes = spt.svg2paths(filepath, convert_rectangles_to_paths=False)
 
-    counts_by_stroke = defaultdict(int)
-    counts_by_style = defaultdict(int)
+    counts_by_stroke: dict[str | None, int] = defaultdict(int)
+    counts_by_style: dict[str, int] = defaultdict(int)
 
     # step 1:
 
-    points = []
-    edges = []
+    points_list: list[NDArray[np.float32]] = []
+    edges: list[tuple[int, int, dict[str, Any]]] = []
 
     for path, attrs in zip(paths, attributes):
         for line in path:
@@ -67,7 +70,7 @@ def get_stroke(attrs: dict) -> str | None:
             start = np.array([line.start.real, line.start.imag], dtype=np.float32)
             end = np.array([line.end.real, line.end.imag], dtype=np.float32)
             # print(start, end)
-            crease_type = None
+            crease_type: int | None = None
             # this works for cps exported from oripa
             if "style" in attrs:
                 if "red" in attrs["style"]:
@@ -76,7 +79,7 @@ def get_stroke(attrs: dict) -> str | None:
                     crease_type = VALLEY
                 elif "gray" in attrs["style"]:
                     continue
-            edge_attrs = dict() if crease_type is None else {CREASE_ASSIGNMENT: crease_type}
+            edge_attrs: dict[str, Any] = dict() if crease_type is None else {CREASE_ASSIGNMENT: crease_type}
 
             # this is e.g. for robert langs cps
             if len(path) == 1:
@@ -89,9 +92,9 @@ def get_stroke(attrs: dict) -> str | None:
                     counts_by_stroke[attrs["stroke"]] += 1
                     edge_attrs["svg_stroke"] = attrs["stroke"]
 
-            edges.append((len(points), len(points) + 1, edge_attrs))
-            points.extend([start, end])
-    points = np.stack(points)
+            edges.append((len(points_list), len(points_list) + 1, edge_attrs))
+            points_list.extend([start, end])
+    points: NDArray[np.float32] = np.stack(points_list)
     points -= np.mean(points, axis=0)
     points /= 2 * np.max(np.abs(points))
     points += [[0.5, 0.5]]
@@ -100,12 +103,14 @@ def get_stroke(attrs: dict) -> str | None:
     first_occurences = np.argmax(clustering[None] == np.arange(np.max(clustering) + 1)[:, None], axis=1)
     merged_points = points[first_occurences]
 
-    G = nx.Graph()
-    G.add_edges_from(
+    nxG: nx.Graph[Any] = nx.Graph()
+    nxG.add_edges_from(
         [(tuple(merged_points[clustering[i]]), tuple(merged_points[clustering[j]]), attrs) for i, j, attrs in edges]
     )
 
-    G = ec.conversions.EHEG_from_nx(G)
+    result = ec.conversions.EHEG_from_nx(nxG)
+    assert isinstance(result, ec.half.EuclideanPositionHEG)
+    G: ec.half.EuclideanPositionHEG = result
 
     if len(counts_by_stroke) not in (2, 3):
         logger.warning(
@@ -113,17 +118,18 @@ def get_stroke(attrs: dict) -> str | None:
         )
         pass
     else:
+        crease_strokes: list[str]
         if len(counts_by_stroke) == 3:  # assume one is the border stroke
-            on_border_counts = {key: 0 for key in list(counts_by_stroke.keys()) + [None]}
+            on_border_counts: dict[str | None, int] = {key: 0 for key in list(counts_by_stroke.keys()) + [None]}
             for e in G.border_edges():
                 on_border_counts[e.attributes.get("svg_stroke", None)] += 1
                 on_border_counts[e.rev.attributes.get("svg_stroke", None)] += 1
             del on_border_counts[None]
             border_stroke = max(on_border_counts.items(), key=operator.itemgetter(1))[0]
             del on_border_counts[border_stroke]
-            crease_strokes = sorted(on_border_counts)
+            crease_strokes = sorted(k for k in on_border_counts if k is not None)
         else:  # 2 strokes
-            crease_strokes = sorted(counts_by_stroke)
+            crease_strokes = sorted(k for k in counts_by_stroke if k is not None)
         for e in G.halfedges:
             stroke = e.attributes.get("svg_stroke", None)
             if stroke in crease_strokes:
diff --git a/eucare/utils.py b/eucare/utils.py
index 1950b0f..0764730 100755
--- a/eucare/utils.py
+++ b/eucare/utils.py
@@ -36,10 +36,9 @@ def random_directed_set(edges: HalfEdgeGraph | Iterable[HalfEdge]) -> set[HalfEd
         order of the input determines which side of each pair is kept; the
         result is therefore deterministic only for ordered inputs.
     """
-    if isinstance(edges, HalfEdgeGraph):
-        edges = edges.halfedges
+    halfedges: Iterable[HalfEdge] = edges.halfedges if isinstance(edges, HalfEdgeGraph) else edges
     directed_edges: set[HalfEdge] = set()
-    for e in edges:
+    for e in halfedges:
         if e.rev not in directed_edges:
             directed_edges.add(e)
     return directed_edges
@@ -84,9 +83,10 @@ def print_attribute_info(objs: HalfEdgeGraph | Iterable[AttributeObject]) -> Non
         print_attribute_info(objs.faces)
         return
 
-    counter = defaultdict(int)
-    attribute_dict = defaultdict(set)
-    for obj in objs:
+    obj_list = list(objs)
+    counter: dict[str, int] = defaultdict(int)
+    attribute_dict: dict[str, set[object]] = defaultdict(set)
+    for obj in obj_list:
         assert isinstance(obj, AttributeObject)
         for key, val in obj.attributes.items():
             try:
@@ -94,6 +94,6 @@ def print_attribute_info(objs: HalfEdgeGraph | Iterable[AttributeObject]) -> Non
             except TypeError:
                 pass
             counter[key] += 1
-    logger.info("%d Objects", len(objs))
+    logger.info("%d Objects", len(obj_list))
     for key, count in sorted(counter.items()):
         logger.info("Key '%s': %d objects (%d distinct hashable values)", key, count, len(attribute_dict[key]))
diff --git a/pyproject.toml b/pyproject.toml
index fe4359c..34a0517 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -26,7 +26,20 @@ threed = ["meshio"]
 torch = ["torch", "einops"]
 image = ["scikit-image", "mahotas"]
 intersecting_cylinders = ["rdp", "plotly"]
-dev = ["pytest", "pytest-cov", "black", "mypy", "pre-commit", "nbstripout", "rdp", "plotly"]
+dev = [
+    "pytest",
+    "pytest-cov",
+    "black",
+    "mypy",
+    "pre-commit",
+    "nbstripout",
+    "rdp",
+    "plotly",
+    "types-PyYAML",
+    "types-networkx",
+    "types-tqdm",
+    "scipy-stubs",
+]
 docs = [
     "mkdocs-materialx",
     "mkdocstrings[python]",
@@ -111,16 +124,16 @@ extend-exclude = '''
 python_version = "3.10"
 warn_unused_ignores = true
 warn_redundant_casts = true
-# Start permissive; tighten per module via overrides as typing lands.
-ignore_missing_imports = false
-disallow_untyped_defs = false
+ignore_missing_imports = true
+disallow_untyped_defs = true
 check_untyped_defs = true
 no_implicit_optional = true
+exclude = ["tests/", "docs/"]
 
 [[tool.mypy.overrides]]
-# Modules that have been fully annotated — promote to strict here.
-module = []
-disallow_untyped_defs = true
+# tifffile uses Python 3.12+ `type` statement syntax; skip parsing it under py3.10.
+module = "tifffile.*"
+follow_imports = "skip"
 
 # to make it so torch cpu is installed