neurospatial.layout.engines.triangular_mesh

triangular_mesh

Classes

TriangularMeshLayout

TriangularMeshLayout()

Bases: _NeighborPairsMixin

A LayoutEngine that builds a triangular mesh over interior points (auto-generated) clipped to a boundary polygon. Each triangle whose centroid lies inside the polygon is kept as an active bin. Connectivity by shared faces.

Source code in src/neurospatial/layout/engines/triangular_mesh.py

def __init__(self) -> None:
    self.bin_centers = np.empty((0, 2), dtype=float)
    self.dimension_ranges = None
    self.grid_edges = ()  # For non-grid layouts, this is typically empty
    self.grid_shape = None
    self.active_mask = None  # Will be updated in build
    self._build_params_used = {}

    self._full_delaunay_tri = None
    self._active_original_simplex_indices = None
    self._boundary_polygon_stored = None

Functions

build

build(boundary_polygon: Polygon, *, point_spacing: float) -> None

Build the triangular mesh layout.

Parameters:

Name	Type	Description	Default
boundary_polygon	Polygon	The polygon defining the boundary. Triangles with centroids outside this polygon are discarded.	required
point_spacing	float	Spacing of the square lattice of interior sample points, in the units of the polygon coordinates. Every ring of the polygon is also sampled at most this far apart.	required

Raises:

Type	Description
TypeError	If boundary_polygon is not a Shapely Polygon.
ValueError	If boundary_polygon is empty; if point_spacing is not a single positive, finite number; if fewer than 3 lattice points fall inside the polygon (checked before the boundary points are added, so a polygon narrower than point_spacing is refused); if the points cannot be triangulated; or if no triangle's centroid lies inside the polygon.

Source code in src/neurospatial/layout/engines/triangular_mesh.py

@capture_build_params(data=("boundary_polygon",))
def build(self, boundary_polygon: Polygon, *, point_spacing: float) -> None:
    """Build the triangular mesh layout.

    Parameters
    ----------
    boundary_polygon : shapely.geometry.Polygon
        The polygon defining the boundary. Triangles with centroids
        outside this polygon are discarded.
    point_spacing : float
        Spacing of the square lattice of interior sample points, in the
        units of the polygon coordinates. Every ring of the polygon is
        also sampled at most this far apart.

    Raises
    ------
    TypeError
        If `boundary_polygon` is not a Shapely Polygon.
    ValueError
        If `boundary_polygon` is empty; if `point_spacing` is not a
        single positive, finite number; if fewer than 3 lattice points
        fall inside the polygon (checked before the boundary points are
        added, so a polygon narrower than `point_spacing` is refused);
        if the points cannot be triangulated; or if no triangle's
        centroid lies inside the polygon.

    """
    self._begin_build()
    if not isinstance(boundary_polygon, Polygon):
        raise TypeError("boundary_polygon must be a Shapely Polygon.")
    if boundary_polygon.is_empty:
        raise ValueError("boundary_polygon cannot be empty.")
    point_spacing = float(
        validate_cell_size(point_spacing, name="point_spacing", scalar_only=True)[0]
    )

    self._boundary_polygon_stored = boundary_polygon

    # 1. Generate sample points for triangulation
    sample_points = _generate_interior_points_for_mesh(
        boundary_polygon,
        point_spacing,
    )
    if sample_points.shape[0] < 3:  # Delaunay needs at least N+1 points in N-D
        raise ValueError(
            f"Not enough interior sample points ({sample_points.shape[0]}) generated "
            "to form any triangle. Try decreasing point_spacing or ensuring "
            "the polygon is large enough relative to the spacing.",
        )
    # 2. Triangulate the lattice with every ring of the polygon as sides
    self._full_delaunay_tri = _conforming_triangulation(
        sample_points,
        boundary_polygon,
        point_spacing,
        triangulate=_triangulate_points,
    )

    # 3. Filter active simplices (triangles)
    active_original_indices, all_centroids = _filter_active_simplices_by_centroid(
        self._full_delaunay_tri,
        boundary_polygon,
    )
    n_total_delaunay_triangles = self._full_delaunay_tri.simplices.shape[0]

    if active_original_indices.size == 0:
        raise ValueError(
            "No triangles found with centroids inside the boundary polygon. "
            "Check boundary_polygon shape, point_spacing, or point generation strategy.",
        )

    self._active_original_simplex_indices = active_original_indices

    # 4. Populate core attributes for active triangles
    self.bin_centers = all_centroids[active_original_indices]

    # 5. Pair active triangles that share a side
    pairs = mesh_pairs(active_original_indices, self._full_delaunay_tri.neighbors)
    self._set_neighbor_pairs(
        pairs, euclidean_pair_distances(self.bin_centers, pairs)
    )

    # 6. Cell extent: the mesh boundary's bounds
    self.dimension_ranges = component_extent(
        outline=self.outline(), shapes=self.bin_shapes()
    )

    # 7. The full "grid" is the list of all Delaunay triangles; the
    # active mask marks the kept ones.
    self.grid_shape = (n_total_delaunay_triangles,)
    self.active_mask = np.zeros(n_total_delaunay_triangles, dtype=bool)
    self.active_mask[active_original_indices] = True
    self._own_arrays("bin_centers", "active_mask")
    self._locator = MeshLocator(
        self._full_delaunay_tri.points,
        self._full_delaunay_tri.simplices,
        mask_active_map(self.active_mask),
        boundary_polygon,
    )

binning

binning() -> Binning

Return owned triangles, lookup and lazy shared-edge faces.

Raises:

Type	Description
LayoutNotBuiltError	If the layout has not been built.

Source code in src/neurospatial/layout/engines/triangular_mesh.py

def binning(self) -> Binning:
    """Return owned triangles, lookup and lazy shared-edge faces.

    Raises
    ------
    LayoutNotBuiltError
        If the layout has not been built.
    """
    tri = self._full_delaunay_tri
    if tri is None or self._active_original_simplex_indices is None:
        raise LayoutNotBuiltError(type(self).__name__, "binning")
    pairs, distances = self._neighbor_pairs()
    return assemble_binning(
        pairs,
        distances,
        kind=self.layout_type,
        centers=self.bin_centers,
        sizes=self.bin_sizes(),
        measure_dim=self.bin_centers.shape[1],
        locator=self.locator(),
        faces="mesh",
        face_arrays=(
            tri.points,
            tri.simplices,
            np.asarray(self._active_original_simplex_indices),
        ),
        shapes=self.bin_shapes(),
        outline=self.outline(),
        dimension_ranges=self.dimension_ranges,
    )

bin_shapes

bin_shapes() -> TriangleShapes

Return mesh points and the active triangles in bin order.

Source code in src/neurospatial/layout/engines/triangular_mesh.py

def bin_shapes(self) -> TriangleShapes:
    """Return mesh points and the active triangles in bin order."""
    if (
        self._full_delaunay_tri is None
        or self._active_original_simplex_indices is None
    ):
        raise LayoutNotBuiltError(type(self).__name__, "bin_shapes")
    return TriangleShapes(
        self._full_delaunay_tri.points,
        self._full_delaunay_tri.simplices[self._active_original_simplex_indices],
    )

outline

outline() -> Outline | None

Return the mesh boundary, including its interior rings when present.

Source code in src/neurospatial/layout/engines/triangular_mesh.py

def outline(self) -> Outline | None:
    """Return the mesh boundary, including its interior rings when present."""
    return (
        Outline(self._boundary_polygon_stored)
        if self._boundary_polygon_stored is not None
        else None
    )

bin_sizes

bin_sizes() -> NDArray[np.float64]

Return the N-dimensional volume of each active N-simplex (bin).

For 2D, this is area. For 3D, this is volume, and so on. The volume of an N-simplex with vertices v0, v1, ..., vn is calculated as: V = (1 / n!) * |det([v1-v0, v2-v0, ..., vn-v0])| where n is the number of dimensions.

Returns:

Type	Description
NDArray[float64]	Array of N-dimensional volumes, shape (n_active_simplices,).

Raises:

Type	Description
LayoutNotBuiltError	If the layout has not been built.

Source code in src/neurospatial/layout/engines/triangular_mesh.py

def bin_sizes(self) -> NDArray[np.float64]:
    """Return the N-dimensional volume of each active N-simplex (bin).

    For 2D, this is area. For 3D, this is volume, and so on.
    The volume of an N-simplex with vertices v0, v1, ..., vn is calculated as:
    V = (1 / n!) * |det([v1-v0, v2-v0, ..., vn-v0])|
    where n is the number of dimensions.

    Returns
    -------
    NDArray[np.float64]
        Array of N-dimensional volumes, shape (n_active_simplices,).

    Raises
    ------
    LayoutNotBuiltError
        If the layout has not been built.

    """
    if (
        self._full_delaunay_tri is None
        or self._active_original_simplex_indices is None
    ):
        raise LayoutNotBuiltError(type(self).__name__, "bin_sizes")

    # Get the vertices of all active N-simplices
    # .points has shape (total_points, n_dim)
    # .simplices has shape (total_simplices, n_dim + 1)
    all_mesh_points = self._full_delaunay_tri.points
    active_simplices_vertex_indices = self._full_delaunay_tri.simplices[
        self._active_original_simplex_indices
    ]

    # nsimplex_vertices will have shape:
    # (num_active_simplices, num_vertices_per_simplex, n_dim)
    nsimplex_vertices = all_mesh_points[active_simplices_vertex_indices]

    n_dim = all_mesh_points.shape[1]

    # Select one vertex from each simplex as the origin (v0)
    # v0 will have shape (num_active_simplices, n_dim)
    v0 = nsimplex_vertices[:, 0, :]

    # Create vectors from v0 to all other vertices (v1-v0, v2-v0, ..., vn-v0)
    # These vectors form the rows (or columns) of the matrix for the determinant.
    # nsimplex_vertices[:, 1:, :] has shape (n_active, n_dim, n_dim)
    # v0[:, np.newaxis, :] broadcasts v0 to match for subtraction
    # matrix_for_determinant will have shape (n_active, n_dim, n_dim)
    matrix_for_determinant = nsimplex_vertices[:, 1:, :] - v0[:, np.newaxis, :]

    # Calculate the determinant for each simplex's matrix
    # np.linalg.det operates on the last two axes by default.
    determinants = np.linalg.det(matrix_for_determinant)

    # Volume = abs(determinant) / n!
    n_factorial = float(math.factorial(n_dim))
    volumes = np.abs(determinants) / n_factorial

    return np.asarray(volumes, dtype=np.float64)

Functions