Source code for nanomesh.image2mesh._mesher2d._mesher

from __future__ import annotations

from typing import TYPE_CHECKING, List

import matplotlib.pyplot as plt
import numpy as np
from scipy.spatial.distance import cdist
from skimage import measure

from nanomesh._constants import BACKGROUND, FEATURE
from nanomesh._doc import doc
from nanomesh.region_markers import RegionMarker, RegionMarkerList

from .._mesher import Mesher
from ._helpers import append_to_segment_markers, generate_segment_markers, pad
from ._polygon import Polygon

    from ..image import Plane
    from ..mesh import LineMesh
    from ..mesh_container import MeshContainer

def _polygons_to_line_mesh(polygons: List[Polygon],
                           bbox: np.ndarray) -> LineMesh:
    """Generate line-mesh from polygons and surrounding bbox. The polygons are
    stacked and missing corners are obtained from the bounding box coordinates.

    polygons : List[Polygon]
        List of polygons.
    bbox : (n, 2) numpy.ndarray
        Coordinates for the bounding box. These define the convex hull
        of the meshing area.

    points : (m,2) numpy.ndarray
        List of points.
    segments : (n,2) numpy.ndarray
        List of segments.
    from nanomesh import LineMesh

    segments = _generate_segments(polygons)

    all_points = np.vstack([polygon.points for polygon in polygons])

    corner_idx = np.argwhere(cdist(bbox, all_points) == 0)

    if len(corner_idx) < len(bbox):
        # Add missing corners and add them where necessary
        missing_corners = np.delete(bbox, corner_idx[:, 0], axis=0)
        all_points = np.vstack([all_points, missing_corners])
        corner_idx = np.argwhere(cdist(bbox, all_points) == 0)

    R = corner_idx[:, 1].tolist()
    additional_segments = list(zip(R, R[1:] + R[:1]))
    segments = np.vstack([segments, additional_segments])

    segment_markers = generate_segment_markers(polygons)

    fields = {}
    for i in np.unique(segment_markers):
        fields[f'L{i}'] = i

    segment_markers = append_to_segment_markers(segment_markers,

    mesh = LineMesh(points=all_points,

    return mesh

def _generate_background_region(polygons: List[Polygon],
                                bbox: np.ndarray) -> RegionMarker:
    """Generate marker for background. This point is inside the bbox, but
    outside the given polygons.

    polygons : List[Polygon]
        List of polygons.
    bbox : (n, 2) numpy.ndarray
        Coordinates for the bounding box. These define the convex hull
        of the meshing area.

    region : RegionMarker
        Region marker to describe the background feature
    point = bbox.mean(axis=0)

    xmin, ymin = bbox.min(axis=0)
    xmax, ymax = bbox.max(axis=0)

    while any(polygon.contains_point(point) for polygon in polygons):
        point = np.random.uniform(xmin, xmax), np.random.uniform(ymin, ymax)

    return RegionMarker(label=BACKGROUND, point=point, name='background')

def _generate_regions(polygons: List[Polygon]) -> RegionMarkerList:
    """Generate regions for triangle.

    polygons : List[Polygon]
        List of polygons.

    regions : RegionMarkerList
        List of region markers describing each feature
    regions = RegionMarkerList()

    for i, polygon in enumerate(polygons):
        point = polygon.find_point()

        regions.append(RegionMarker(label=FEATURE, point=point, name='X'))

    return regions

def _generate_segments(polygons: List[Polygon]) -> np.ndarray:
    """Generate segments for triangle.

    polygons : List[Polygon]
        List of polygons

    segments : numpy.ndarray
        Segment connectivity array
    i = 0
    segments = []

    for polygon in polygons:
        n_points = len(polygon)
        rng = np.arange(i, i + n_points)

        # generate segment connectivity matrix
        segment = np.vstack([rng, np.roll(rng, shift=-1)]).T

        i += n_points

    return np.vstack(segments)

[docs]@doc(Mesher, prefix='triangular mesh from 2D image data') class Mesher2D(Mesher, ndim=2): def __init__(self, image: np.ndarray | Plane): super().__init__(image) self.contour: LineMesh self._bbox = None
[docs] def generate_contour( self, level: float = None, precision: int = 1, max_edge_dist: int = 5, group_regions: bool = True, ): """Generate contours using marching cubes algorithm. Contours are approximated by a polygon, where the maximum distance between points is decided by `max_edge_dist`. Parameters ---------- level : float, optional Contour value to search for isosurfaces (i.e. the threshold value). By default takes the average of the min and max value. Can be ignored if a binary image is passed to :class:`Mesher2D`. precision : int, optional Maximum distance from original points in polygon approximation routine. max_edge_dist : int, optional Divide long edges so that maximum distance between points does not exceed this value. group_regions : bool, optional If True, assign the same label to all features If False, label regions sequentially """ polygons = [ Polygon(points) for points in measure.find_contours(self.image, level=level) ] polygons = [polygon.approximate(precision) for polygon in polygons] polygons = [ polygon.subdivide(max_dist=max_edge_dist) for polygon in polygons ] polygons = [ polygon.close_corner(self.image.shape) for polygon in polygons ] polygons = [polygon.remove_duplicate_points() for polygon in polygons] # remove polygons with no area, # fixes polygons = [polygon for polygon in polygons if len(polygon) > 2] regions = _generate_regions(polygons) regions.append(_generate_background_region(polygons, self.bbox)) if not group_regions: regions = regions.label_sequentially(FEATURE, fmt_name='X{}') contour = _polygons_to_line_mesh(polygons, self.bbox) contour.region_markers = regions self.polygons = polygons self.contour = contour
@property def image_bbox(self) -> np.ndarray: """Return bbox from image shape. Returns ------- bbox : (4,2) np.array Coordinates of bounding box contour. """ x, y = self.image.shape return np.array(( (0, 0), (x - 1, 0), (x - 1, y - 1), (0, y - 1), )) @property def bbox(self) -> np.ndarray: """Return bbox attribute. If not explicity set, returns :attr:`Mesher2D.image_bbox`. Sequence: x0, y0 x1, y0 x1, y1 x0, y0 Returns ------- bbox : np.ndarray Bounding box set for output mesh. """ if self._bbox is None: return self.image_bbox else: return self._bbox @bbox.setter def bbox(self, bbox: np.ndarray): """Set bounding box attribute. Parameters ---------- bbox : np.ndarray List of coordinates for bounding box corners: x0, y0 x1, y0 x1, y1 x0, y0 Raises ------ ValueError Raised if `bbox` has the wrong shape. """ bbox = np.array(bbox) if bbox.shape != (4, 2): raise ValueError('Bounding box must be an array with shape (4,2).') self._bbox = bbox
[docs] def triangulate(self, opts='pAq30a10', **kwargs) -> MeshContainer: """Triangulate contours. Mandatory switches for `opts`: - `e`: ensure edges get returned - `p`: planar straight line graph - `A`: assign regional attribute to each triangle If missing, these will be added. Parameters ---------- opts : str, optional Options passed to :func:`triangulate`. For more info, see: Returns ------- mesh : MeshContainer Triangulated 2D mesh with domain labels """ default_opts = {'p': True, 'A': True, 'e': True} mesh = self.contour.triangulate(opts=opts, default_opts=default_opts, **kwargs) return mesh
[docs] @doc(pad, prefix='Pad the contour using :func:`image2mesh._mesher2d.pad`') def pad_contour(self, **kwargs): self.contour = pad(self.contour, **kwargs)
[docs] def plot_contour(self, ax: plt.Axes = None, cmap: str = None, **kwargs): """Plot contours on image. Parameters ---------- ax : matplotlib.axes.Axes Axes to use for plotting. cmap : str Matplotlib color map for :func:`matplotlib.pyplot.imshow` **kwargs These parameters are passed to :func:`plotting.linemeshplot` Returns ------- ax : matplotlib.axes.Axes """ if not ax: fig, ax = plt.subplots() ax.set_title('Contours') self.contour.plot_mpl(ax=ax, **kwargs) ax.imshow(self.image, cmap=cmap) ax.axis('image') ax.set_xticks([]) ax.set_yticks([]) return ax
[docs]def plane2mesh(image: np.ndarray | Plane, *, level: float = None, max_edge_dist: int = 5, opts: str = 'q30a10', plot: bool = False) -> 'MeshContainer': """Generate a triangular mesh from a 2D segmented image. Parameters ---------- image : (i,j) numpy.ndarray or Plane Input image to mesh. level : float, optional Level to generate contours at from image max_edge_dist : int, optional Maximum distance between neighbouring pixels in contours. opts : str, optional Options passed to :func:`triangulate`. For more info, see: Returns ------- mesh : MeshContainer Triangulated 2D mesh with domain labels. """ mesher = Mesher2D(image) mesher.generate_contour(max_edge_dist=5, level=level) return mesher.triangulate(opts=opts)