diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index f684e2a..cc4fc11 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -32,7 +32,7 @@ repos:
     rev: 6.1.0  # Use the latest stable version
     hooks:
       - id: flake8
-        args: [--max-line-length=88, '--ignore=D205,D400,D105,E731,W503,D401,E203,D209']
+        args: [--max-line-length=88, '--ignore=D205,D400,D105,E731,W503,D401,E203,D209,D202']
         additional_dependencies: [flake8-docstrings]
 
   - repo: https://github.com/pre-commit/mirrors-mypy
diff --git a/pyalphashape/AlphaShape.py b/pyalphashape/AlphaShape.py
index 061c446..f7736be 100644
--- a/pyalphashape/AlphaShape.py
+++ b/pyalphashape/AlphaShape.py
@@ -20,7 +20,7 @@
   α‑shape in *any* dimension, supports strict/relaxed connectivity rules,
   optional hole‑patching, incremental point insertion, inside/point‑to‑surface
   queries, centroid computation, and access to perimeter vertices, edges and
-  faces.
+  facets.
 
 In practice you construct an α‑shape from an ``(N,d)`` array of points,
 tune α until the boundary is as tight as required, and then use the resulting
@@ -33,13 +33,12 @@
 from collections import defaultdict
 from scipy.spatial import Delaunay
 import numpy as np
-from typing import Tuple, Set, List, Literal, Optional
+from typing import Tuple, Set, List, Literal, Optional, Generator, Union, Dict, Any
 from pyalphashape.GraphClosure import GraphClosureTracker
 
 
 def circumcenter(points: np.ndarray) -> np.ndarray:
-    """
-    Compute the circumcenter of a simplex in arbitrary dimensions.
+    """Compute the circumcenter of a simplex in arbitrary dimensions.
 
     Parameters
     ----------
@@ -51,23 +50,22 @@ def circumcenter(points: np.ndarray) -> np.ndarray:
     -------
     np.ndarray
         The barycentric coordinates of the circumcenter of the simplex.
+
     """
 
     n, _ = points.shape
 
     # build the (d+1) × (d+1) system with plain ndarrays
-    A = np.block([
-        [2 * points @ points.T, np.ones((n, 1))],
-        [np.ones((1, n)), np.zeros((1, 1))]
-    ])
+    A = np.block(
+        [[2 * points @ points.T, np.ones((n, 1))], [np.ones((1, n)), np.zeros((1, 1))]]
+    )
     b = np.concatenate([np.sum(points * points, axis=1), np.array([1.0])])
 
     return np.linalg.solve(A, b)[:-1]
 
 
 def circumradius(points: np.ndarray) -> float:
-    """
-    Compute the circumradius of a simplex in arbitrary dimensions.
+    """Compute the circumradius of a simplex in arbitrary dimensions.
 
     Parameters
     ----------
@@ -78,15 +76,17 @@ def circumradius(points: np.ndarray) -> float:
     -------
     float
         The circumradius of the simplex.
+
     """
 
     return np.linalg.norm(points[0, :] - np.dot(circumcenter(points), points))
 
 
-
-def alphasimplices(points: np.ndarray) -> Tuple[np.ndarray, float, np.ndarray]:
-    """
-    Generate all simplices in the Delaunay triangulation along with their circumradii.
+def alphasimplices(
+    points: np.ndarray,
+) -> Generator[Tuple[np.ndarray, float, np.ndarray], None, None]:
+    """Generate all simplices in the Delaunay triangulation along with their
+    circumradii.
 
     Parameters
     ----------
@@ -100,6 +100,7 @@ def alphasimplices(points: np.ndarray) -> Tuple[np.ndarray, float, np.ndarray]:
         - The simplex as an array of indices,
         - Its circumradius,
         - The coordinates of the simplex vertices.
+
     """
 
     coords = np.asarray(points)
@@ -110,33 +111,82 @@ def alphasimplices(points: np.ndarray) -> Tuple[np.ndarray, float, np.ndarray]:
         try:
             yield simplex, circumradius(simplex_points), simplex_points
         except np.linalg.LinAlgError:
-            logging.warn('Singular matrix. Likely caused by all points '
-                         'lying in an N-1 space.')
+            logging.warn(
+                "Singular matrix. Likely caused by all points " "lying in an N-1 space."
+            )
 
 
-class AlphaShape:
-    """
-    Compute the α-shape (concave hull) of a point cloud in arbitrary dimensions.
+def point_to_simplex_distance(
+    p: np.ndarray, simplex: np.ndarray, tol: float = 1e-9
+) -> float:
+    """Compute the shortest Euclidean distance from a point to a k‑simplex in
+    n‑dimensional space.
 
     Parameters
     ----------
-    points : np.ndarray
-        An (N, d) array of points.
-    alpha : float, optional
-        The α parameter controlling the "tightness" of the shape. Default is 0.
-    connectivity : {"strict", "relaxed"}, optional
-        Connectivity rule for filtering simplices. Default is "strict".
-    ensure_closure : bool, default True
-        If True, triangles that are otherwise too large but are fully enclosed by
-        accepted triangles will be included. This prevents holes
+    p : np.ndarray of shape (n,)
+        A query point in the ambient n‑dimensional space.
+    simplex : np.ndarray of shape (k+1, n)
+        An array of k+1 vertices defining a k‑dimensional simplex embedded in n‑space.
+    tol : float, optional
+        Tolerance for the barycentric coordinate “inside” test. Defaults to 1e-9.
+
+    Returns
+    -------
+    float
+        The Euclidean distance from `p` to the simplex. Returns 0.0 if `p` projects
+        inside or on the simplex (within the given tolerance).
+
     """
+    k = len(simplex) - 1
+    if k == 0:
+        # 0-simplex: just a point
+        return np.linalg.norm(p - simplex[0])
+
+    base = simplex[0]
+    A = (simplex[1:] - base).T  # shape (d, k)
+    λ, *_ = np.linalg.lstsq(A, p - base, rcond=None)
+    proj = base + A @ λ
+    bary = np.empty(k + 1)
+    bary[0] = 1 - λ.sum()
+    bary[1:] = λ
+
+    # inside (or on) the simplex?
+    if np.all(bary >= -tol):
+        return np.linalg.norm(p - proj)
+
+    # otherwise, recurse on all k facets (remove one vertex at a time)
+    return min(
+        point_to_simplex_distance(p, np.delete(simplex, i, axis=0), tol)
+        for i in range(k + 1)
+    )
+
 
-    def __init__(self,
-                 points: np.ndarray,
-                 alpha: float = 0.,
-                 connectivity: Literal["strict", "relaxed"] = "strict",
-                 ensure_closure: bool = True
+class AlphaShape:
+    """Compute the α-shape (concave hull) of a point cloud in arbitrary dimensions."""
+
+    def __init__(
+        self,
+        points: np.ndarray,
+        alpha: float = 0.0,
+        connectivity: Literal["strict", "relaxed"] = "strict",
+        ensure_closure: bool = True,
     ):
+        """Compute the α-shape (concave hull) of a point cloud in arbitrary dimensions.
+
+        Parameters
+        ----------
+        points : np.ndarray
+            An (N, d) array of points.
+        alpha : float, optional
+            The α parameter controlling the "tightness" of the shape. Default is 0.
+        connectivity : {"strict", "relaxed"}, optional
+            Connectivity rule for filtering simplices. Default is "strict".
+        ensure_closure : bool, default True
+            If True, triangles that are otherwise too large but are fully enclosed by
+            accepted triangles will be included. This prevents holes
+
+        """
         self._dim = points.shape[1]
         if self._dim < 2:
             raise ValueError("dimension must be ≥ 2")
@@ -151,7 +201,7 @@ def __init__(self,
 
         self.simplices: Set[Tuple[int, ...]] = set()
         self.perimeter_edges: List[Tuple[np.ndarray, np.ndarray]] = []
-        self.perimeter_points: np.ndarray | None = None
+        self.perimeter_points: Union[np.ndarray, None] = None
         self.GCT = GraphClosureTracker(len(points))
 
         # build once
@@ -159,20 +209,19 @@ def __init__(self,
 
     @property
     def vertices(self) -> Optional[np.ndarray]:
-        """
-        Get the perimeter vertices of the alpha shape.
+        """Get the perimeter vertices of the alpha shape.
 
         Returns
         -------
         np.ndarray or None
             Array of perimeter points, or None if not computed.
+
         """
 
         return self.perimeter_points
 
     def contains_point(self, pt: np.ndarray, tol: float = 1e-8) -> bool:
-        """
-        Check whether a given point lies inside or on the alpha shape.
+        """Check whether a given point lies inside or on the alpha shape.
 
         Parameters
         ----------
@@ -185,9 +234,10 @@ def contains_point(self, pt: np.ndarray, tol: float = 1e-8) -> bool:
         -------
         bool
             True if the point lies inside or on the alpha shape; False otherwise.
+
         """
 
-        if len(self.perimeter_points) == 0:
+        if self.perimeter_points is None or len(self.perimeter_points) == 0:
             return False
 
         # 1. Close to any perimeter vertex?
@@ -214,7 +264,8 @@ def contains_point(self, pt: np.ndarray, tol: float = 1e-8) -> bool:
                 A = np.vstack([verts.T, np.ones(len(verts))])  # (d+1, d+1)
                 b = np.append(pt, 1.0)
 
-                # Full‑rank simplex → solve directly; otherwise fall back to least‑squares
+                # Full‑rank simplex → solve directly;
+                # otherwise fall back to least‑squares
                 if A.shape[0] == A.shape[1]:
                     bary = np.linalg.solve(A, b)
                 else:
@@ -228,8 +279,7 @@ def contains_point(self, pt: np.ndarray, tol: float = 1e-8) -> bool:
         return False
 
     def add_points(self, new_pts: np.ndarray, perimeter_only: bool = False) -> None:
-        """
-        Add new points to the alpha shape (batch rebuild).
+        """Add new points to the alpha shape (batch rebuild).
 
         Parameters
         ----------
@@ -237,59 +287,60 @@ def add_points(self, new_pts: np.ndarray, perimeter_only: bool = False) -> None:
             A (N, d) array of new points to add. The alpha shape is rebuilt.
         perimeter_only: bool
             If True, only pass perimeter points to new shape. Otherwise, pass all points
+
         """
 
         if perimeter_only:
             pts = np.vstack([self.points, new_pts])
         else:
             pts = np.vstack([self.perimeter_points, new_pts])
-        self.__init__(pts, alpha=self.alpha, connectivity=self.connectivity,
-                      ensure_closure=self.ensure_closure)
-
-    def _get_boundary_faces(self) -> Set[Tuple[int, ...]]:
-        """
-        Identify and return the boundary (d-1)-faces of the alpha shape.
-
-        Returns
-        -------
-        Set[Tuple[int, ...]]
-            A set of index tuples representing the boundary faces.
-        """
-
-        if hasattr(self, "_boundary_faces"):
-            return self._boundary_faces
+        self.__init__(  # type: ignore[misc]
+            pts,
+            alpha=self.alpha,
+            connectivity=self.connectivity,
+            ensure_closure=self.ensure_closure,
+        )
+
+    def _get_boundary_facets(self) -> Set[Tuple[int, ...]]:
+        """Identify and return the boundary (d-1)-facets of the alpha shape, and cache
+        per-facet filter data for fast distance queries."""
+        if hasattr(self, "_boundary_facets"):
+            return self._boundary_facets
 
         dim = self._dim
-        faces: Set[Tuple[int, ...]] = set()
+        facets: Set[Tuple[int, ...]] = set()
         for s in self.simplices:
             for f in itertools.combinations(s, dim):
                 f = tuple(sorted(f))
-                if f in faces:
-                    faces.remove(f)
+                if f in facets:
+                    facets.remove(f)
                 else:
-                    faces.add(f)
-        # cache
-        self._boundary_faces = faces
-        return faces
+                    facets.add(f)
+        # cache boundary facets
+        self._boundary_facets: Set[Tuple[int, ...]] = facets
+
+        # cache filter data: (normal, base_point, sphere_center, sphere_radius)
+        self._facet_filter_data: Dict[
+            Tuple[int, ...], Tuple[np.ndarray, np.ndarray, np.ndarray, float]
+        ] = {}
+        for facet in facets:
+            verts = self.points[list(facet)]  # (d, d)
+            base = verts[0]
+            A = (verts[1:] - base).T  # (d, d-1)
+            # unit normal via nullspace of A^T
+            _, _, Vt = np.linalg.svd(A.T)
+            normal = Vt[-1]
+            normal /= np.linalg.norm(normal)
+            # bounding sphere: centroid + max vertex distance
+            center = verts.mean(axis=0)
+            radius = np.max(np.linalg.norm(verts - center, axis=1))
+            self._facet_filter_data[facet] = (normal, base, center, radius)
+
+        return facets
 
     def distance_to_surface(self, point: np.ndarray, tol: float = 1e-9) -> float:
-        """
-        Compute the shortest Euclidean distance from a point to the alpha shape surface.
-
-        Parameters
-        ----------
-        point : np.ndarray
-            A point of shape (d,) in the same ambient space as the alpha shape.
-        tol : float, optional
-            Tolerance for barycentric coordinate test. Default is 1e-9.
-
-        Returns
-        -------
-        float
-            Distance from the point to the alpha shape surface. Returns 0 if inside
-            or on surface.
-        """
-
+        """Compute the shortest Euclidean distance from a point to the alpha shape
+        surface, using plane & sphere filters for speed."""
         p = np.asarray(point, dtype=float)
         if p.shape[-1] != self._dim:
             raise ValueError("point dimensionality mismatch")
@@ -298,41 +349,34 @@ def distance_to_surface(self, point: np.ndarray, tol: float = 1e-9) -> float:
         if self.contains_point(p):
             return 0.0
 
-        # 2. gather boundary faces and vertices
-        faces = self._get_boundary_faces()
-        if not faces:
-            # degenerate case (e.g. only 1–2 input points)
-            # fall back to nearest perimeter vertex
-            return np.min(np.linalg.norm(self.perimeter_points - p, axis=1))
-
-        dists = []
-
-        for f in faces:
-            verts = self.points[list(f)]  # shape (d, d)
-            base = verts[0]
-            A = verts[1:] - base  # (d‑1, d)
+        # 2. boundary facets
+        facets = self._get_boundary_facets()
+        if not facets:
+            return float(np.min(np.linalg.norm(self.perimeter_points - p, axis=1)))
 
-            # orthogonal projection of p onto the face’s affine span
-            # Solve A x = (p - base)  →  least‑squares because A is tall
-            x_hat, *_ = np.linalg.lstsq(A.T, (p - base), rcond=None)
-            proj = base + A.T @ x_hat
-
-            # barycentric coordinates to test if proj is inside the simplex
-            # coords = [1 - sum(x_hat), *x_hat]
-            bary = np.concatenate(([1.0 - x_hat.sum()], x_hat))
-            if np.all(bary >= -tol):  # inside (or on) the face
-                dists.append(np.linalg.norm(p - proj))
-            else:
-                # outside → distance to nearest vertex of this face
-                dists.extend(np.linalg.norm(verts - p, axis=1))
+        best = np.inf
+        for facet in facets:
+            normal, base, center, radius = self._facet_filter_data[facet]
+            # plane-distance filter
+            if abs(np.dot(p - base, normal)) >= best:
+                continue
+            # bounding-sphere filter
+            if np.linalg.norm(p - center) - radius >= best:
+                continue
+            # exact distance on this facet
+            verts = self.points[list(facet)]
+            d = point_to_simplex_distance(p, verts, tol)
+            if d < best:
+                best = d
 
-        return float(min(dists))
+        return float(best)
 
     def _build_batch(self) -> None:
-        """
-        Construct the alpha shape using Delaunay triangulation and filtering by alpha.
+        """Construct the alpha shape using Delaunay triangulation and filtering by
+        alpha.
 
         This method is automatically called upon initialization.
+
         """
 
         dim, pts = self._dim, self.points
@@ -356,8 +400,9 @@ def _build_batch(self) -> None:
 
         for simp, r in simplices:
             root_set = {uf.find(v) for v in simp}
-            keep = (r <= r_filter) or \
-                   (self.connectivity == "relaxed" and len(root_set) > 1)
+            keep = (r <= r_filter) or (
+                self.connectivity == "relaxed" and len(root_set) > 1
+            )
             if not keep:
                 continue
             uf.add_fully_connected_subgraph(list(simp))
@@ -375,15 +420,15 @@ def _build_batch(self) -> None:
 
             # Identify the largest connected component
             comp_sizes = {root: len(nodes) for root, nodes in gct.components.items()}
-            main_root = max(comp_sizes, key=comp_sizes.get)
+            main_root = max(comp_sizes, key=lambda k: comp_sizes[k])
             main_verts = gct.components[main_root]
 
             # Build edge-to-triangle map
             edge_to_triangles = defaultdict(list)
             for simp in all_passed:
                 for edge in itertools.combinations(simp, 2):
-                    edge = tuple(sorted(edge))
-                    edge_to_triangles[edge].append(simp)
+                    edge_ = tuple(sorted(edge))
+                    edge_to_triangles[edge_].append(simp)
 
             # Keep only triangles in the main component that share an edge with another
             kept = []
@@ -415,13 +460,13 @@ def _build_batch(self) -> None:
         # ---------- 3.  rebuild perimeter from *kept* simplices ----------
         self.simplices = set(kept)
         self.GCT = GraphClosureTracker(n)  # final tracker
-        edge_counts = defaultdict(int)
+        edge_counts: Dict[Tuple[Any, ...], int] = defaultdict(int)
 
         for s in self.simplices:
             self.GCT.add_fully_connected_subgraph(list(s))
             for edge in itertools.combinations(s, 2):  # triangle edges
-                edge = tuple(sorted(edge))
-                edge_counts[edge] += 1
+                edge_ = tuple(sorted(edge))
+                edge_counts[edge_] += 1
 
         # ---------- 4.  store perimeter ----------------------------------
         perimeter_edges_idx = [e for e, count in edge_counts.items() if count == 1]
@@ -432,41 +477,41 @@ def _build_batch(self) -> None:
 
     @property
     def is_empty(self) -> bool:
-        """
-        Check whether the alpha shape is empty (i.e., no perimeter points).
+        """Check whether the alpha shape is empty (i.e., no perimeter points).
 
         Returns
         -------
         bool
             True if empty; False otherwise.
-        """
 
+        """
+        if self.perimeter_points is None:
+            return True
         return len(self.perimeter_points) == 0
 
     @property
-    def triangle_faces(self) -> List[np.ndarray]:
-        """
-        Get the triangle faces (simplices) that make up the alpha shape.
+    def triangle_facets(self) -> List[np.ndarray]:
+        """Get the triangle facets (simplices) that make up the alpha shape.
 
         Returns
         -------
         List[np.ndarray]
             List of arrays containing vertex coordinates of each simplex.
+
         """
 
         return [self.points[list(s)] for s in self.simplices]
 
-
     @property
     def centroid(self) -> np.ndarray:
-        """
-        Compute the hyper-volumetric centroid of the Euclidean alpha shape
-        using direct determinant-based simplex volume.
+        """Compute the hyper-volumetric centroid of the Euclidean alpha shape using
+        direct determinant-based simplex volume.
 
         Returns
         -------
         np.ndarray
             A (d,) array representing the centroid in Euclidean space.
+
         """
         if len(self.simplices) == 0:
             return np.full(self._dim, np.nan)
@@ -495,4 +540,3 @@ def centroid(self) -> np.ndarray:
             return np.mean(self.points, axis=0)
 
         return weighted_sum / total_volume
-
diff --git a/pyalphashape/GraphClosure.py b/pyalphashape/GraphClosure.py
index e5e1bc9..d6ce8b9 100644
--- a/pyalphashape/GraphClosure.py
+++ b/pyalphashape/GraphClosure.py
@@ -35,23 +35,25 @@
 number or composition of connected components in real time.
 """
 
-from typing import List
+from typing import List, Union, Tuple, Set
+
+
 class GraphClosureTracker:
-    """
-    A dynamic Union-Find (Disjoint Set) data structure with explicit tracking
-    of connected components and support for dynamic resizing.
+    """A dynamic Union-Find (Disjoint Set) data structure with explicit tracking of
+    connected components and support for dynamic resizing.
 
     Useful for managing dynamic connectivity in undirected graphs.
+
     """
 
     def __init__(self, num_nodes: int):
-        """
-        Initialize the tracker with a specified number of nodes.
+        """Initialize the tracker with a specified number of nodes.
 
         Parameters
         ----------
         num_nodes : int
             The initial number of nodes in the graph.
+
         """
 
         self.parent = list(range(num_nodes))
@@ -59,15 +61,14 @@ def __init__(self, num_nodes: int):
         self.num_nodes = num_nodes
         self.components = {i: {i} for i in range(num_nodes)}  # Initial components
 
-
     def _ensure_capacity(self, node: int) -> None:
-        """
-        Dynamically expand internal arrays to include a node with the given index.
+        """Dynamically expand internal arrays to include a node with the given index.
 
         Parameters
         ----------
         node : int
             The node index to ensure capacity for.
+
         """
 
         if node >= self.num_nodes:
@@ -80,8 +81,7 @@ def _ensure_capacity(self, node: int) -> None:
 
     # modify the public methods to call it
     def find(self, node: int) -> int:
-        """
-        Find the root representative of the set containing the node.
+        """Find the root representative of the set containing the node.
 
         Parameters
         ----------
@@ -92,6 +92,7 @@ def find(self, node: int) -> int:
         -------
         int
             The root node of the component.
+
         """
 
         self._ensure_capacity(node)
@@ -100,8 +101,7 @@ def find(self, node: int) -> int:
         return self.parent[node]
 
     def union(self, node1: int, node2: int) -> None:
-        """
-        Merge the components containing node1 and node2.
+        """Merge the components containing node1 and node2.
 
         Parameters
         ----------
@@ -109,6 +109,7 @@ def union(self, node1: int, node2: int) -> None:
             First node.
         node2 : int
             Second node.
+
         """
 
         root1 = self.find(node1)
@@ -132,8 +133,7 @@ def union(self, node1: int, node2: int) -> None:
                 del self.components[root2]
 
     def add_edge(self, node1: int, node2: int) -> None:
-        """
-        Add an undirected edge between two nodes by merging their components.
+        """Add an undirected edge between two nodes by merging their components.
 
         Parameters
         ----------
@@ -141,18 +141,21 @@ def add_edge(self, node1: int, node2: int) -> None:
             First node.
         node2 : int
             Second node.
+
         """
 
         self.union(node1, node2)
 
-    def add_fully_connected_subgraph(self, nodes: List[int]) -> None:
-        """
-        Fully connect a list of nodes by merging all pairs into one component.
+    def add_fully_connected_subgraph(
+        self, nodes: Union[List[int], Tuple[int, ...]]
+    ) -> None:
+        """Fully connect a list of nodes by merging all pairs into one component.
 
         Parameters
         ----------
         nodes : List[int]
             A list of node indices to be fully connected.
+
         """
 
         for i in range(len(nodes)):
@@ -160,8 +163,7 @@ def add_fully_connected_subgraph(self, nodes: List[int]) -> None:
                 self.union(nodes[i], nodes[j])
 
     def subgraph_is_already_connected(self, nodes: List[int]) -> bool:
-        """
-        Check whether all nodes in the list belong to the same connected component.
+        """Check whether all nodes in the list belong to the same connected component.
 
         Parameters
         ----------
@@ -172,6 +174,7 @@ def subgraph_is_already_connected(self, nodes: List[int]) -> bool:
         -------
         bool
             True if all nodes are connected, False otherwise.
+
         """
 
         if not nodes:
@@ -182,8 +185,7 @@ def subgraph_is_already_connected(self, nodes: List[int]) -> bool:
         return all(self.find(node) == root for node in nodes)
 
     def is_connected(self, node1: int, node2: int) -> bool:
-        """
-        Check whether two nodes are in the same connected component.
+        """Check whether two nodes are in the same connected component.
 
         Parameters
         ----------
@@ -196,25 +198,25 @@ def is_connected(self, node1: int, node2: int) -> bool:
         -------
         bool
             True if node1 and node2 are connected, False otherwise.
+
         """
 
         return self.find(node1) == self.find(node2)
 
     def __iter__(self):
-        """
-        Iterate over the current connected components.
+        """Iterate over the current connected components.
 
         Returns
         -------
         Iterator[Set[int]]
             An iterator over sets of node indices.
+
         """
 
         return iter(self.components.values())
 
-    def __getitem__(self, index: int) -> List[int]:
-        """
-        Index into the list of connected components.
+    def __getitem__(self, index: int) -> Set[int]:
+        """Index into the list of connected components.
 
         Parameters
         ----------
@@ -225,18 +227,19 @@ def __getitem__(self, index: int) -> List[int]:
         -------
         List[int]
             List of nodes in the selected connected component.
+
         """
 
         return list(self.components.values())[index]
 
     def __len__(self) -> int:
-        """
-        Return the number of connected components.
+        """Return the number of connected components.
 
         Returns
         -------
         int
             The number of components currently being tracked.
+
         """
 
-        return len(self.components)
\ No newline at end of file
+        return len(self.components)
diff --git a/pyalphashape/SphericalAlphaShape.py b/pyalphashape/SphericalAlphaShape.py
index d1e6c77..a9889e2 100644
--- a/pyalphashape/SphericalAlphaShape.py
+++ b/pyalphashape/SphericalAlphaShape.py
@@ -18,7 +18,7 @@
 - Provides geometric queries:
   * `contains_point` – point‑in‑shape test.
   * `distance_to_surface` – shortest arc distance to the boundary.
-  * `triangle_faces`, `triangle_faces_latlon` – access to all retained
+  * `triangle_facets`, `triangle_facets_latlon` – access to all retained
     triangles in either coordinate system.
   * `centroid` – area‑weighted centre of the shape.
 - Accepts incremental point insertion via `add_points`.
@@ -48,41 +48,46 @@
 import numpy as np
 import itertools
 from collections import defaultdict
-from typing import Literal, Set, Tuple, List, Optional
+from typing import Literal, Set, Tuple, List, Optional, Union, Any, Dict
 from pyalphashape.SphericalDelaunay import SphericalDelaunay
 from pyalphashape.sphere_utils import (
     latlon_to_unit_vectors,
     unit_vectors_to_latlon,
     arc_distance,
     spherical_circumradius,
-    spherical_triangle_area
+    spherical_triangle_area,
 )
 from pyalphashape.GraphClosure import GraphClosureTracker
 
 
 class SphericalAlphaShape:
-    """
-    Compute the α-shape (concave hull) of points defined on the surface of a unit sphere.
-
-    Parameters
-    ----------
-    points : np.ndarray
-        An (N, 2) array of latitude and longitude coordinates in degrees.
-    alpha : float, optional
-        The α parameter controlling shape detail. Smaller values yield tighter shapes.
-    connectivity : {"strict", "relaxed"}, optional
-        Rule for keeping connected components during filtering.
-    ensure_closure : bool, default True
-        If True, triangles that are otherwise too large but are fully enclosed by
-        accepted triangles will be included. This prevents holes
-    """
-
-    def __init__(self,
-                 points: np.ndarray,
-                 alpha: float = 0.,
-                 connectivity: Literal["strict", "relaxed"] = "strict",
-                 ensure_closure: bool = True
+    """Compute the α-shape (concave hull) of points defined on the surface of a unit
+    sphere."""
+
+    def __init__(
+        self,
+        points: np.ndarray,
+        alpha: float = 0.0,
+        connectivity: Literal["strict", "relaxed"] = "strict",
+        ensure_closure: bool = True,
     ):
+        """Compute the α-shape (concave hull) of points defined on the surface of a unit
+        sphere.
+
+        Parameters
+        ----------
+        points : np.ndarray
+            An (N, 2) array of latitude and longitude coordinates in degrees.
+        alpha : float, optional
+            The α parameter controlling shape detail. Smaller values yield tighter
+            shapes.
+        connectivity : {"strict", "relaxed"}, optional
+            Rule for keeping connected components during filtering.
+        ensure_closure : bool, default True
+            If True, triangles that are otherwise too large but are fully enclosed by
+            accepted triangles will be included. This prevents holes
+
+        """
         self._dim = points.shape[1]
         if self._dim < 2:
             raise ValueError("dimension must be ≥ 2")
@@ -98,8 +103,8 @@ def __init__(self,
 
         self.simplices: Set[Tuple[int, ...]] = set()
         self.perimeter_edges: List[Tuple[np.ndarray, np.ndarray]] = []
-        self.perimeter_points: np.ndarray | None = None
-        self.perimeter_points_latlon: np.ndarray | None = None
+        self.perimeter_points: Union[np.ndarray, None] = None
+        self.perimeter_points_latlon: Union[np.ndarray, None] = None
         self.GCT = GraphClosureTracker(len(points))
 
         # build once
@@ -107,21 +112,21 @@ def __init__(self,
 
     @property
     def vertices(self) -> Optional[np.ndarray]:
-        """
-        Return the perimeter vertices of the spherical alpha shape (in 3D unit vector form).
+        """Return the perimeter vertices of the spherical alpha shape (in 3D unit vector
+        form).
 
         Returns
         -------
         np.ndarray or None
             The perimeter points of the alpha shape, or None if uninitialized.
+
         """
 
         return self.perimeter_points
 
     def contains_point(self, pt_latlon: np.ndarray, tol: float = 1e-8) -> bool:
-        """
-        Return ``True`` iff the latitude/longitude point lies inside *or on*
-        the spherical α‑shape.
+        """Return ``True`` iff the latitude/longitude point lies inside *or on* the
+        spherical α‑shape.
 
         Parameters
         ----------
@@ -129,20 +134,19 @@ def contains_point(self, pt_latlon: np.ndarray, tol: float = 1e-8) -> bool:
             [latitude, longitude] **degrees**.
         tol : float, default 1e‑8
             Angular tolerance (radians) for vertex proximity and half‑space test.
+
         """
         # ── 0. quick outs ──────────────────────────────────────────────────
-        if len(self.perimeter_points) == 0:
+        if self.perimeter_points is None or len(self.perimeter_points) == 0:
             return False
 
         P = latlon_to_unit_vectors(pt_latlon[None, :])[0]  # unit vector
 
         # close to a perimeter vertex?
-        if np.any(
-                np.arccos(np.clip(self.perimeter_points @ P, -1.0, 1.0)) < tol
-        ):
+        if np.any(np.arccos(np.clip(self.perimeter_points @ P, -1.0, 1.0)) < tol):
             return True
 
-        # no faces yet  → no enclosed area
+        # no facets yet  → no enclosed area
         if len(self.simplices) == 0:
             return False
 
@@ -178,8 +182,7 @@ def contains_point(self, pt_latlon: np.ndarray, tol: float = 1e-8) -> bool:
         return False
 
     def add_points(self, new_pts: np.ndarray, perimeter_only: bool = False) -> None:
-        """
-        Add new latitude-longitude points and rebuild the spherical alpha shape.
+        """Add new latitude-longitude points and rebuild the spherical alpha shape.
 
         Parameters
         ----------
@@ -187,42 +190,47 @@ def add_points(self, new_pts: np.ndarray, perimeter_only: bool = False) -> None:
             An (M, 2) array of new points in degrees [lat, lon].
         perimeter_only: bool
             If True, only pass perimeter points to new shape. Otherwise, pass all points
+
         """
         if perimeter_only:
             pts = np.vstack([self.perimeter_points_latlon, new_pts])
         else:
             pts = np.vstack([self.points_latlon, new_pts])
-        self.__init__(pts, alpha=self.alpha, connectivity=self.connectivity,
-                      ensure_closure=self.ensure_closure)
+        self.__init__(  # type: ignore[misc]
+            pts,
+            alpha=self.alpha,
+            connectivity=self.connectivity,
+            ensure_closure=self.ensure_closure,
+        )
 
-    def _get_boundary_faces(self) -> Set[Tuple[int, ...]]:
-        """
-        Return the set of (d-1)-faces that form the boundary of the alpha shape.
+    def _get_boundary_facets(self) -> Set[Tuple[int, ...]]:
+        """Return the set of (d-1)-facets that form the boundary of the alpha shape.
 
         Returns
         -------
         Set[Tuple[int, ...]]
-            Set of index tuples representing the boundary faces.
+            Set of index tuples representing the boundary facets.
+
         """
 
-        if hasattr(self, "_boundary_faces"):
-            return self._boundary_faces
+        if hasattr(self, "_boundary_facets"):
+            return self._boundary_facets
 
-        faces: Set[Tuple[int, ...]] = set()
+        facets: Set[Tuple[int, ...]] = set()
         for s in self.simplices:
             for f in itertools.combinations(s, 2):
-                f = tuple(sorted(f))
-                if f in faces:
-                    faces.remove(f)
+                f_: tuple[int, ...] = tuple(sorted(f))
+                if f_ in facets:
+                    facets.remove(f_)
                 else:
-                    faces.add(f)
+                    facets.add(f_)
         # cache
-        self._boundary_faces = faces
-        return faces
+        self._boundary_facets: Set[Tuple[int, ...]] = facets
+        return facets
 
     def distance_to_surface(self, point: np.ndarray) -> float:
-        """
-        Compute the angular arc distance from a (lat, lon) point to the alpha shape surface.
+        """Compute the angular arc distance from a (lat, lon) point to the alpha shape
+        surface.
 
         Parameters
         ----------
@@ -234,28 +242,36 @@ def distance_to_surface(self, point: np.ndarray) -> float:
         float
             The arc distance in radians from the point to the alpha shape surface.
             Returns 0 if the point lies inside or on the surface.
+
         """
 
         if point.shape[-1] != 2:
             raise ValueError("Input point must be (lat, lon) in degrees")
 
+        if self.perimeter_points is None:
+            raise ValueError("Cannot find distance to surface with no perimeter points")
+
         if self.contains_point(point):
             return 0.0
 
         # Convert (lat, lon) to 3D unit vector
         p = latlon_to_unit_vectors(point[None, :])[0]
 
-        faces = self._get_boundary_faces()
-        if not faces:
+        facets = self._get_boundary_facets()
+        if not facets:
             # Fallback to nearest perimeter point
-            return float(np.min([
-                np.arccos(np.clip(np.dot(p, q), -1.0, 1.0))
-                for q in self.perimeter_points
-            ]))
+            return float(
+                np.min(
+                    [
+                        np.arccos(np.clip(np.dot(p, q), -1.0, 1.0))
+                        for q in self.perimeter_points
+                    ]
+                )
+            )
 
         # Compute minimum arc distance to all boundary edges
         dists = []
-        for f in faces:
+        for f in facets:
             idx = list(f)
             if len(idx) != 2:
                 raise ValueError(
@@ -267,10 +283,11 @@ def distance_to_surface(self, point: np.ndarray) -> float:
         return float(min(dists))
 
     def _build_batch(self) -> None:
-        """
-        Construct the spherical alpha shape by computing Delaunay triangles
-        and filtering them by circumradius and connectivity.
+        """Construct the spherical alpha shape by computing Delaunay triangles and
+        filtering them by circumradius and connectivity.
+
         This method is automatically called on initialization.
+
         """
 
         pts, pts_latlon = self.points, self.points_latlon
@@ -299,8 +316,9 @@ def _build_batch(self) -> None:
 
         for simp, r in simplices:
             root_set = {uf.find(v) for v in simp}
-            keep = (r <= r_filter) or \
-                   (self.connectivity == "relaxed" and len(root_set) > 1)
+            keep = (r <= r_filter) or (
+                self.connectivity == "relaxed" and len(root_set) > 1
+            )
             if not keep:
                 continue
             uf.add_fully_connected_subgraph(list(simp))
@@ -318,16 +336,15 @@ def _build_batch(self) -> None:
 
             # Identify the largest connected component
             comp_sizes = {root: len(nodes) for root, nodes in gct.components.items()}
-            main_root = max(comp_sizes, key=comp_sizes.get)
+            main_root = max(comp_sizes, key=lambda k: comp_sizes[k])
             main_verts = gct.components[main_root]
 
             # Build edge-to-triangle map
             edge_to_triangles = defaultdict(list)
             for simp in all_passed:
                 for edge in itertools.combinations(simp, 2):
-                    edge = tuple(sorted(edge))
-                    edge_to_triangles[edge].append(simp)
-
+                    edge_: Tuple[Any, ...] = tuple(sorted(edge))
+                    edge_to_triangles[edge_].append(simp)
 
             # Keep only triangles that:
             # (a) have all vertices in the main component, and
@@ -364,12 +381,12 @@ def _build_batch(self) -> None:
         self.simplices = set(kept)
         self.GCT = GraphClosureTracker(n)  # final tracker
 
-        edge_counts = defaultdict(int)
+        edge_counts: Dict[Tuple[Any, ...], int] = defaultdict(int)
         for s in self.simplices:
             self.GCT.add_fully_connected_subgraph(list(s))
             for edge in itertools.combinations(s, 2):  # triangle edges
-                edge = tuple(sorted(edge))
-                edge_counts[edge] += 1
+                edge_ = tuple(sorted(edge))
+                edge_counts[edge_] += 1
 
         # ---------- 4.  store perimeter ----------------------------------
         # Only edges that appear once are on the perimeter
@@ -379,57 +396,61 @@ def _build_batch(self) -> None:
         self.perimeter_points = pts[list(sorted(perim_idx))]
         self.perimeter_points_latlon = pts_latlon[list(sorted(perim_idx))]
         self.perimeter_edges = [(pts[i], pts[j]) for i, j in perimeter_edges_idx]
-        self.perimeter_edges_latlon = [(pts_latlon[i], pts_latlon[j]) for i, j in
-                                       perimeter_edges_idx]
+        self.perimeter_edges_latlon = [
+            (pts_latlon[i], pts_latlon[j]) for i, j in perimeter_edges_idx
+        ]
 
     @property
     def is_empty(self) -> bool:
-        """
-        Check whether the alpha shape has any perimeter points.
+        """Check whether the alpha shape has any perimeter points.
 
         Returns
         -------
         bool
             True if no perimeter has been constructed, False otherwise.
-        """
 
+        """
+        if self.perimeter_points is None:
+            return True
         return len(self.perimeter_points) == 0
 
     @property
-    def triangle_faces(self) -> List[np.ndarray]:
-        """
-        Return all triangles (simplices) of the alpha shape in unit vector coordinates.
+    def triangle_facets(self) -> List[np.ndarray]:
+        """Return all triangles (simplices) of the alpha shape in unit vector
+        coordinates.
 
         Returns
         -------
         List[np.ndarray]
             List of (3, 3) arrays representing triangle vertices as 3D unit vectors.
+
         """
 
         return [self.points[list(s)] for s in self.simplices]
 
     @property
-    def triangle_faces_latlon(self) -> List[np.ndarray]:
-        """
-        Return all triangles (simplices) of the alpha shape in (latitude, longitude) degrees.
+    def triangle_facets_latlon(self) -> List[np.ndarray]:
+        """Return all triangles (simplices) of the alpha shape in (latitude, longitude)
+        degrees.
 
         Returns
         -------
         List[np.ndarray]
             List of (3, 2) arrays representing triangle vertices in (lat, lon).
+
         """
 
         return [self.points_latlon[list(s)] for s in self.simplices]
 
     @property
     def centroid(self) -> np.ndarray:
-        """
-        Compute the center of area of the spherical alpha shape.
+        """Compute the center of area of the spherical alpha shape.
 
         Returns
         -------
         np.ndarray
             A (2,) array representing the centroid in [latitude, longitude] degrees.
+
         """
         if len(self.simplices) == 0:
             return np.array([np.nan, np.nan])
@@ -449,5 +470,3 @@ def centroid(self) -> np.ndarray:
 
         centroid_vec /= np.linalg.norm(centroid_vec)
         return unit_vectors_to_latlon(centroid_vec[None])[0]
-
-
diff --git a/pyalphashape/plotting.py b/pyalphashape/plotting.py
index a2191ed..4d7f181 100644
--- a/pyalphashape/plotting.py
+++ b/pyalphashape/plotting.py
@@ -36,37 +36,37 @@
 ```
 """
 
-
 import matplotlib.pyplot as plt
 from pyalphashape.AlphaShape import AlphaShape
 from pyalphashape.SphericalAlphaShape import SphericalAlphaShape
 from pyalphashape.SphericalDelaunay import SphericalDelaunay
-from typing import Any, Optional, Tuple, List, Dict
+from typing import Any, Optional, Tuple, List, Dict, cast
 from matplotlib.axes import Axes
 from mpl_toolkits.mplot3d import Axes3D  # noqa
 import numpy as np
 import plotly.graph_objects as go
 
+
 def plot_polygon_edges(
     edges: np.ndarray,
     ax: Axes,
     line_color: str = "b",
     line_width: float = 1.0,
 ):
-    """
-    Plot a set of polygon edges on a 2D or 3D Matplotlib axis.
+    """Plot a set of polygon edges on a 2D or 3D Matplotlib axis.
 
     Parameters
     ----------
     edges : np.ndarray
-        An array of shape (n, 2) where each entry is a pair of points (each a coordinate array)
-        representing an edge of the polygon.
+        An array of shape (n, 2) where each entry is a pair of points
+        (each a coordinate array) representing an edge of the polygon.
     ax : matplotlib.axes.Axes
         A Matplotlib Axes object (2D or 3D) to plot on.
     line_color : str, optional
         Color of the edges, by default 'b'.
     line_width : float, optional
         Width of the edge lines, by default 1.0.
+
     """
     for p1, p2 in edges:
         x = [p1[0], p2[0]]
@@ -91,9 +91,11 @@ def plot_polygon_edges(
             )
 
 
-def interpolate_great_arc(A: np.ndarray, B: np.ndarray, num_points: int = 100) -> np.ndarray:
-    """
-    Compute evenly spaced points along the great arc connecting two points on the unit sphere.
+def interpolate_great_arc(
+    A: np.ndarray, B: np.ndarray, num_points: int = 100
+) -> np.ndarray:
+    """Compute evenly spaced points along the great arc connecting two points on the
+    unit sphere.
 
     Parameters
     ----------
@@ -107,7 +109,9 @@ def interpolate_great_arc(A: np.ndarray, B: np.ndarray, num_points: int = 100) -
     Returns
     -------
     np.ndarray
-        Array of shape (num_points, 3) containing interpolated points on the unit sphere.
+        Array of shape (num_points, 3) containing interpolated points on the unit
+        sphere.
+
     """
 
     A = A / np.linalg.norm(A)
@@ -120,8 +124,9 @@ def interpolate_great_arc(A: np.ndarray, B: np.ndarray, num_points: int = 100) -
 
     sin_theta = np.sin(theta)
     t_vals = np.linspace(0, 1, num_points)
-    arc_points = (np.sin((1 - t_vals) * theta)[:, None] * A +
-                  np.sin(t_vals * theta)[:, None] * B) / sin_theta
+    arc_points = (
+        np.sin((1 - t_vals) * theta)[:, None] * A + np.sin(t_vals * theta)[:, None] * B
+    ) / sin_theta
     return arc_points
 
 
@@ -135,10 +140,9 @@ def plot_spherical_triangulation(
     marker_symbol: Optional[str] = "circle",
     line_width: float = 1.5,
     line_color: Any = "blue",
-    line_style: str = "-"
+    line_style: str = "-",
 ):
-    """
-    Visualize a spherical Delaunay triangulation using either Matplotlib or Plotly.
+    """Visualize a spherical Delaunay triangulation using either Matplotlib or Plotly.
 
     Parameters
     ----------
@@ -167,6 +171,7 @@ def plot_spherical_triangulation(
     -------
     plotly.graph_objects.Figure or matplotlib.axes.Axes
         The figure or axis object used for plotting.
+
     """
 
     if ax is not None:
@@ -175,7 +180,7 @@ def plot_spherical_triangulation(
         ax.grid(False)
 
         # Sphere surface
-        u, v = np.mgrid[0:2*np.pi:60j, 0:np.pi:30j]
+        u, v = np.mgrid[0 : 2 * np.pi : 60j, 0 : np.pi : 30j]  # type: ignore[misc]
         xs = np.cos(u) * np.sin(v)
         ys = np.sin(u) * np.sin(v)
         zs = np.cos(v)
@@ -183,7 +188,9 @@ def plot_spherical_triangulation(
 
         # Points
         pts = triangulation.points_xyz
-        ax.scatter(pts[:, 0], pts[:, 1], pts[:, 2], color=marker_color, s=marker_size**2)
+        ax.scatter(
+            pts[:, 0], pts[:, 1], pts[:, 2], color=marker_color, s=marker_size**2
+        )
 
         # Arcs
         for tri in triangulation.triangles:
@@ -192,15 +199,26 @@ def plot_spherical_triangulation(
                 A = pts[tri[indices[i]]]
                 B = pts[tri[indices[i + 1]]]
                 arc_pts = interpolate_great_arc(A, B)
-                ax.plot(arc_pts[:, 0], arc_pts[:, 1], arc_pts[:, 2],
-                        color=line_color, linewidth=line_width, linestyle=line_style)
+                ax.plot(
+                    arc_pts[:, 0],
+                    arc_pts[:, 1],
+                    arc_pts[:, 2],
+                    color=line_color,
+                    linewidth=line_width,
+                    linestyle=line_style,
+                )
 
         return ax
 
     return _plot_spherical_triangulation_plotly(
-        triangulation, title, fig,
-        marker_size, marker_color, marker_symbol,
-        line_width, line_color
+        triangulation,
+        title,
+        fig,
+        marker_size,
+        marker_color,
+        marker_symbol,
+        line_width,
+        line_color,
     )
 
 
@@ -212,10 +230,9 @@ def _plot_spherical_triangulation_plotly(
     marker_color: Any,
     marker_symbol: Optional[str],
     line_width: float,
-    line_color: Any
+    line_color: Any,
 ):
-    """
-    Internal helper to render a spherical triangulation using Plotly.
+    """Internal helper to render a spherical triangulation using Plotly.
 
     Parameters
     ----------
@@ -240,30 +257,39 @@ def _plot_spherical_triangulation_plotly(
     -------
     plotly.graph_objects.Figure
         The updated or newly created Plotly figure.
+
     """
 
     if fig is None:
         fig = go.Figure()
 
-        u, v = np.mgrid[0:2*np.pi:60j, 0:np.pi:30j]
+        u, v = np.mgrid[0 : 2 * np.pi : 60j, 0 : np.pi : 30j]  # type: ignore[misc]
         xs = np.cos(u) * np.sin(v)
         ys = np.sin(u) * np.sin(v)
         zs = np.cos(v)
-        fig.add_trace(go.Surface(
-            x=xs, y=ys, z=zs,
-            opacity=0.15,
-            showscale=False,
-            colorscale='Greys',
-            name='Sphere'
-        ))
-
-    pts = triangulation.perimeter_points
-    fig.add_trace(go.Scatter3d(
-        x=pts[:, 0], y=pts[:, 1], z=pts[:, 2],
-        mode='markers',
-        marker=dict(size=marker_size, color=marker_color, symbol=marker_symbol),
-        name='Points'
-    ))
+        fig.add_trace(
+            go.Surface(
+                x=xs,
+                y=ys,
+                z=zs,
+                opacity=0.15,
+                showscale=False,
+                colorscale="Greys",
+                name="Sphere",
+            )
+        )
+
+    pts = triangulation.points_xyz
+    fig.add_trace(
+        go.Scatter3d(
+            x=pts[:, 0],
+            y=pts[:, 1],
+            z=pts[:, 2],
+            mode="markers",
+            marker=dict(size=marker_size, color=marker_color, symbol=marker_symbol),
+            name="Points",
+        )
+    )
 
     for tri in triangulation.triangles:
         indices = [0, 1, 2, 0]
@@ -271,20 +297,26 @@ def _plot_spherical_triangulation_plotly(
             A = pts[tri[indices[i]]]
             B = pts[tri[indices[i + 1]]]
             arc_pts = interpolate_great_arc(A, B)
-            fig.add_trace(go.Scatter3d(
-                x=arc_pts[:, 0], y=arc_pts[:, 1], z=arc_pts[:, 2],
-                mode='lines',
-                line=dict(color=line_color, width=line_width),
-                showlegend=False
-            ))
+            fig.add_trace(
+                go.Scatter3d(
+                    x=arc_pts[:, 0],
+                    y=arc_pts[:, 1],
+                    z=arc_pts[:, 2],
+                    mode="lines",
+                    line=dict(color=line_color, width=line_width),
+                    showlegend=False,
+                )
+            )
 
     fig.update_layout(
         title=title,
-        scene=dict(xaxis=dict(showgrid=False),
-                   yaxis=dict(showgrid=False),
-                   zaxis=dict(showgrid=False),
-                   aspectmode='data'),
-        margin=dict(l=0, r=0, b=0, t=30)
+        scene=dict(
+            xaxis=dict(showgrid=False),
+            yaxis=dict(showgrid=False),
+            zaxis=dict(showgrid=False),
+            aspectmode="data",
+        ),
+        margin=dict(l=0, r=0, b=0, t=30),
     )
     return fig
 
@@ -293,10 +325,9 @@ def plot_alpha_shape(
     shape: AlphaShape,
     ax: Optional[Axes] = None,
     line_width: int = 1,
-    line_color: Any = "r"
+    line_color: Any = "r",
 ):
-    """
-    Plot a 2D alpha shape perimeter using Matplotlib.
+    """Plot a 2D alpha shape perimeter using Matplotlib.
 
     Parameters
     ----------
@@ -308,32 +339,27 @@ def plot_alpha_shape(
         Width of the perimeter lines, by default 1.
     line_color : Any, optional
         Color of the perimeter lines, by default 'r'.
+
     """
 
     edges = [(e1, e2) for e1, e2 in shape.perimeter_edges]
 
     if ax is None:
         fig, ax = plt.subplots()
-    plot_polygon_edges(
-        edges, ax, line_width=line_width, line_color=line_color
-    )
+    plot_polygon_edges(edges, ax, line_width=line_width, line_color=line_color)
 
 
-import numpy as np
-from typing import Tuple, Dict, List
-
 def generate_geodesic_fill_mesh(
-    shape: 'SphericalAlphaShape',
-    resolution: int = 10
+    shape: "SphericalAlphaShape", resolution: int = 10
 ) -> Tuple[np.ndarray, np.ndarray]:
-    """
-    Generate a unified, watertight geodesic mesh to fill the interior of a spherical alpha shape.
-    All triangle interiors are subdivided and projected back to the sphere to prevent gaps.
+    """Generate a unified, watertight geodesic mesh to fill the interior of a spherical
+    alpha shape. All triangle interiors are subdivided and projected back to the sphere
+    to prevent gaps.
 
     Parameters
     ----------
     shape : SphericalAlphaShape
-        The alpha shape instance with triangle faces.
+        The alpha shape instance with triangle facets.
     resolution : int, optional
         Number of subdivisions per triangle edge. Higher means smoother (default: 10).
 
@@ -341,8 +367,9 @@ def generate_geodesic_fill_mesh(
     -------
     vertices : np.ndarray
         Array of shape (N, 3) of points on the unit sphere.
-    faces : np.ndarray
+    facets : np.ndarray
         Array of shape (M, 3) of triangle indices into `vertices`.
+
     """
 
     def slerp(a: np.ndarray, b: np.ndarray, t: float) -> np.ndarray:
@@ -355,17 +382,17 @@ def slerp(a: np.ndarray, b: np.ndarray, t: float) -> np.ndarray:
 
     vertices: List[np.ndarray] = []
     vertex_cache: Dict[Tuple[float, float, float], int] = {}
-    faces: List[Tuple[int, int, int]] = []
+    facets: List[Tuple[int, int, int]] = []
 
     def add_vertex(v: np.ndarray) -> int:
         v /= np.linalg.norm(v)
-        key = tuple(np.round(v, 8))
+        key = cast(tuple[float, float, float], tuple(np.round(v, 8)))
         if key not in vertex_cache:
             vertex_cache[key] = len(vertices)
             vertices.append(v)
         return vertex_cache[key]
 
-    for triangle in shape.triangle_faces:
+    for triangle in shape.triangle_facets:
         A, B, C = triangle
         index_grid = []
 
@@ -386,12 +413,12 @@ def add_vertex(v: np.ndarray) -> int:
                 v0 = index_grid[i][j]
                 v1 = index_grid[i + 1][j]
                 v2 = index_grid[i + 1][j + 1]
-                faces.append((v0, v1, v2))
+                facets.append((v0, v1, v2))
                 if j < i:
                     v3 = index_grid[i][j + 1]
-                    faces.append((v0, v2, v3))
+                    facets.append((v0, v2, v3))
 
-    return np.array(vertices), np.array(faces)
+    return np.array(vertices), np.array(facets)
 
 
 def plot_spherical_alpha_shape(
@@ -407,10 +434,10 @@ def plot_spherical_alpha_shape(
     ax: Optional[Axes] = None,
     fill: bool = True,
     fill_color: Any = "skyblue",
-    fill_alpha: float = 0.4
+    fill_alpha: float = 0.4,
 ):
-    """
-    Visualize the perimeter of a spherical alpha shape using either Matplotlib or Plotly.
+    """Visualize the perimeter of a spherical alpha shape using either Matplotlib or
+    Plotly.
 
     Parameters
     ----------
@@ -445,6 +472,7 @@ def plot_spherical_alpha_shape(
     -------
     plotly.graph_objects.Figure or matplotlib.axes.Axes
         The figure or axis used for plotting.
+
     """
 
     if ax is not None:
@@ -453,7 +481,7 @@ def plot_spherical_alpha_shape(
         ax.grid(False)
 
         # Sphere surface
-        u, v = np.mgrid[0:2*np.pi:60j, 0:np.pi:30j]
+        u, v = np.mgrid[0 : 2 * np.pi : 60j, 0 : np.pi : 30j]  # type: ignore[misc]
         xs = np.cos(u) * np.sin(v)
         ys = np.sin(u) * np.sin(v)
         zs = np.cos(v)
@@ -463,26 +491,50 @@ def plot_spherical_alpha_shape(
         if fill:
             verts, tris = generate_geodesic_fill_mesh(shape, resolution=20)
             ax.plot_trisurf(
-                verts[:, 0], verts[:, 1], verts[:, 2],
+                verts[:, 0],
+                verts[:, 1],
+                verts[:, 2],
                 triangles=tris,
-                color=fill_color, alpha=fill_alpha, edgecolor="none"
+                color=fill_color,
+                alpha=fill_alpha,
+                edgecolor="none",
             )
 
         # Arcs
         for A, B in shape.perimeter_edges:
             arc_pts = interpolate_great_arc(A, B)
-            ax.plot(arc_pts[:, 0], arc_pts[:, 1], arc_pts[:, 2],
-                    color=line_color, linewidth=line_width, linestyle=line_style)
+            ax.plot(
+                arc_pts[:, 0],
+                arc_pts[:, 1],
+                arc_pts[:, 2],
+                color=line_color,
+                linewidth=line_width,
+                linestyle=line_style,
+            )
 
         # Points
+        assert (
+            shape.perimeter_points is not None
+        ), "Cannot plot without any perimeter points"
         pts = shape.perimeter_points
-        ax.scatter(pts[:, 0], pts[:, 1], pts[:, 2], color=marker_color, s=marker_size**2)
+        ax.scatter(
+            pts[:, 0], pts[:, 1], pts[:, 2], color=marker_color, s=marker_size**2
+        )
 
         return ax
 
     return _plot_spherical_alpha_shape_plotly(
-        shape, line_width, line_color, marker_size, marker_color, marker_symbol,
-        title, fig, fill, fill_color, fill_alpha
+        shape,
+        line_width,
+        line_color,
+        marker_size,
+        marker_color,
+        marker_symbol,
+        title,
+        fig,
+        fill,
+        fill_color,
+        fill_alpha,
     )
 
 
@@ -497,15 +549,15 @@ def _plot_spherical_alpha_shape_plotly(
     fig: Optional[go.Figure],
     fill: bool,
     fill_color: Any,
-    fill_alpha: float
+    fill_alpha: float,
 ):
-    """
-    Internal helper to render a spherical alpha shape using Plotly in 3D.
+    """Internal helper to render a spherical alpha shape using Plotly in 3D.
 
     Parameters
     ----------
     shape : SphericalAlphaShape
-        The spherical alpha shape object containing perimeter edges, simplices, and point cloud.
+        The spherical alpha shape object containing perimeter edges, simplices, and
+        point cloud.
     line_width : float
         Width of the perimeter arcs.
     line_color : Any
@@ -521,7 +573,8 @@ def _plot_spherical_alpha_shape_plotly(
     fig : plotly.graph_objects.Figure or None
         Existing Plotly figure to modify. If None, a new one is created.
     fill : bool
-        Whether to shade the interior area of the alpha shape using the spherical triangles.
+        Whether to shade the interior area of the alpha shape using the spherical
+        triangles.
     fill_color : Any
         Fill color for the triangle mesh.
     fill_alpha : float
@@ -531,56 +584,79 @@ def _plot_spherical_alpha_shape_plotly(
     -------
     plotly.graph_objects.Figure
         The Plotly figure object containing the spherical alpha shape.
+
     """
 
     if fig is None:
         fig = go.Figure()
-        u, v = np.mgrid[0:2*np.pi:60j, 0:np.pi:30j]
+        u, v = np.mgrid[0 : 2 * np.pi : 60j, 0 : np.pi : 30j]  # type: ignore[misc]
         xs = np.cos(u) * np.sin(v)
         ys = np.sin(u) * np.sin(v)
         zs = np.cos(v)
-        fig.add_trace(go.Surface(
-            x=xs, y=ys, z=zs,
-            opacity=0.15,
-            showscale=False,
-            colorscale='Greys',
-            name='Sphere'
-        ))
+        fig.add_trace(
+            go.Surface(
+                x=xs,
+                y=ys,
+                z=zs,
+                opacity=0.15,
+                showscale=False,
+                colorscale="Greys",
+                name="Sphere",
+            )
+        )
 
     # Optional filled triangles
     if fill:
         verts, tris = generate_geodesic_fill_mesh(shape, resolution=20)
-        fig.add_trace(go.Mesh3d(
-            x=verts[:, 0], y=verts[:, 1], z=verts[:, 2],
-            i=tris[:, 0], j=tris[:, 1], k=tris[:, 2],
-            color=fill_color,
-            opacity=fill_alpha,
-            showscale=False
-        ))
+        fig.add_trace(
+            go.Mesh3d(
+                x=verts[:, 0],
+                y=verts[:, 1],
+                z=verts[:, 2],
+                i=tris[:, 0],
+                j=tris[:, 1],
+                k=tris[:, 2],
+                color=fill_color,
+                opacity=fill_alpha,
+                showscale=False,
+            )
+        )
 
     for A, B in shape.perimeter_edges:
         arc_pts = interpolate_great_arc(A, B)
-        fig.add_trace(go.Scatter3d(
-            x=arc_pts[:, 0], y=arc_pts[:, 1], z=arc_pts[:, 2],
-            mode='lines',
-            line=dict(color=line_color, width=line_width),
-            showlegend=False
-        ))
-
+        fig.add_trace(
+            go.Scatter3d(
+                x=arc_pts[:, 0],
+                y=arc_pts[:, 1],
+                z=arc_pts[:, 2],
+                mode="lines",
+                line=dict(color=line_color, width=line_width),
+                showlegend=False,
+            )
+        )
+    assert (
+        shape.perimeter_points is not None
+    ), "Cannot plot without any perimeter points"
     pts = shape.perimeter_points
-    fig.add_trace(go.Scatter3d(
-        x=pts[:, 0], y=pts[:, 1], z=pts[:, 2],
-        mode='markers',
-        marker=dict(size=marker_size, color=marker_color, symbol=marker_symbol),
-        name="Points"
-    ))
+    fig.add_trace(
+        go.Scatter3d(
+            x=pts[:, 0],
+            y=pts[:, 1],
+            z=pts[:, 2],
+            mode="markers",
+            marker=dict(size=marker_size, color=marker_color, symbol=marker_symbol),
+            name="Points",
+        )
+    )
 
     fig.update_layout(
         title=title,
-        scene=dict(xaxis=dict(showgrid=False),
-                   yaxis=dict(showgrid=False),
-                   zaxis=dict(showgrid=False),
-                   aspectmode='data'),
-        margin=dict(l=0, r=0, b=0, t=30)
+        scene=dict(
+            xaxis=dict(showgrid=False),
+            yaxis=dict(showgrid=False),
+            zaxis=dict(showgrid=False),
+            aspectmode="data",
+        ),
+        margin=dict(l=0, r=0, b=0, t=30),
     )
     return fig
diff --git a/unit_tests/test_AlphaShape.py b/unit_tests/test_AlphaShape.py
index 01152d2..fdf4084 100644
--- a/unit_tests/test_AlphaShape.py
+++ b/unit_tests/test_AlphaShape.py
@@ -33,20 +33,56 @@ def test_alpha_shape_basic_perimeter():
     assert isinstance(alpha_shape.perimeter_edges, list)
 
 
-def test_contains_point_inside_and_outside():
+def test_contains_point_inside():
     points = np.array([[0, 0], [1, 0], [0.5, 1]])
     shape = AlphaShape(points, alpha=0.0)
     assert shape.contains_point(np.array([0.5, 0.5]))
-    assert not shape.contains_point(np.array([2, 2]))
-
 
-def test_distance_to_surface_inside_and_outside():
+def test_contains_point_outside():
     points = np.array([[0, 0], [1, 0], [0.5, 1]])
     shape = AlphaShape(points, alpha=0.0)
-    inside_point = np.array([0.5, 0.5])
-    outside_point = np.array([2.0, 2.0])
-    assert np.isclose(shape.distance_to_surface(inside_point), 0.0)
-    assert shape.distance_to_surface(outside_point) > 0.0
+    assert not shape.contains_point(np.array([2, 2]))
+
+def test_distance_to_surface_various_cases():
+    # Define a simple 2D triangle
+    points = np.array([
+        [0.0, 0.0],
+        [1.0, 0.0],
+        [0.5, 1.0],
+    ])
+    shape = AlphaShape(points, alpha=0.0)
+
+    # 1) Interior point ⇒ distance = 0
+    inside = np.array([0.5, 0.5])
+    assert np.isclose(shape.distance_to_surface(inside), 0.0, atol=1e-9)
+
+    # 2) Far above ⇒ closest to vertex (0.5,1.0):
+    #    distance = sqrt((2-0.5)^2 + (2-1)^2) = sqrt(3.25)
+    above = np.array([2.0, 2.0])
+    expected_above = np.sqrt(3.25)
+    got_above = shape.distance_to_surface(above)
+    assert np.isclose(got_above, expected_above, atol=1e-9), (
+        f"Above: expected {expected_above:.9f}, got {got_above:.9f}"
+    )
+
+    # 3) Directly below base edge ⇒ closest to edge y=0 between x=0 and x=1:
+    #    point = (0.5, -1) ⇒ distance = 1.0
+    below = np.array([0.5, -1.0])
+    expected_below = 1.0
+    got_below = shape.distance_to_surface(below)
+    assert np.isclose(got_below, expected_below, atol=1e-9), (
+        f"Below: expected {expected_below:.9f}, got {got_below:.9f}"
+    )
+
+    # 4) Down‑left corner ⇒ closest to vertex (0,0):
+    #    point = (-1, -1) ⇒ distance = sqrt(2)
+    down_left = np.array([-1.0, -1.0])
+    expected_dl = np.sqrt(2)
+    got_dl = shape.distance_to_surface(down_left)
+    assert np.isclose(got_dl, expected_dl, atol=1e-9), (
+        f"Down-left: expected {expected_dl:.9f}, got {got_dl:.9f}"
+    )
+
 
 
 def test_add_points_grows_shape():
@@ -58,12 +94,12 @@ def test_add_points_grows_shape():
     assert n_perim_after > n_perim_before
 
 
-def test_get_boundary_faces_returns_faces():
+def test_get_boundary_facets_returns_facets():
     points = np.random.rand(10, 3)
     shape = AlphaShape(points, alpha=1.0)
-    faces = shape._get_boundary_faces()
-    assert isinstance(faces, set)
-    for f in faces:
+    facets = shape._get_boundary_facets()
+    assert isinstance(facets, set)
+    for f in facets:
         assert isinstance(f, tuple)
         assert len(f) == shape._dim
 
diff --git a/unit_tests/test_SphericalAlphaShape.py b/unit_tests/test_SphericalAlphaShape.py
index 1bce091..33716c2 100644
--- a/unit_tests/test_SphericalAlphaShape.py
+++ b/unit_tests/test_SphericalAlphaShape.py
@@ -56,7 +56,7 @@ def test_add_points_behavior():
     assert n_after > n_before
 
 
-def test_triangle_faces_properties():
+def test_triangle_facets_properties():
     points = np.array([
         [0, 0],
         [0, 90],
@@ -65,11 +65,11 @@ def test_triangle_faces_properties():
         [15, 60]
     ])
     shape = SphericalAlphaShape(points, alpha=1.0)
-    faces_vec = shape.triangle_faces
-    faces_ll = shape.triangle_faces_latlon
-    for face in faces_vec:
+    facets_vec = shape.triangle_facets
+    facets_ll = shape.triangle_facets_latlon
+    for face in facets_vec:
         assert face.shape == (3, 3)  # 3 unit vectors
-    for face in faces_ll:
+    for face in facets_ll:
         assert face.shape == (3, 2)  # 3 lat/lon pairs
 
 
@@ -94,7 +94,7 @@ def test_degenerate_shape_returns_nan_centroid():
     assert np.isnan(centroid).all()
 
 
-def test_get_boundary_faces_valid_format():
+def test_get_boundary_facets_valid_format():
     points = np.array([
         [0, 0],
         [0, 90],
@@ -103,7 +103,7 @@ def test_get_boundary_faces_valid_format():
         [15, 60]
     ])
     shape = SphericalAlphaShape(points, alpha=1.0)
-    faces = shape._get_boundary_faces()
-    for f in faces:
+    facets = shape._get_boundary_facets()
+    for f in facets:
         assert isinstance(f, tuple)
         assert len(f) == 2  # spherical edges are between 2 points