plots of objects defined by shape intersection logic no longer display thin line artifacts. VERSION 2 (issue #2560)

Author: jewettaijfc

CHANGELOG.md

@@ -11,17 +11,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Add support for `np.unwrap` in `tidy3d.plugins.autograd`.
 - Add Nunley variant to germanium material library based on Nunley et al. 2016 data.
 
-### Changed
-- Switched to an analytical gradient calculation for spatially-varying pole-residue models (`CustomPoleResidue`).
-
-### Changed
-- Significantly improved performance of the `tidy3d.plugins.autograd.grey_dilation` morphological operation and its gradient calculation. The new implementation is orders of magnitude faster, especially for large arrays and kernel sizes.
-
 ### Fixed
 - Arrow lengths are now scaled consistently in the X and Y directions, and their lengths no longer exceed the height of the plot window.
 - Bug in `PlaneWave` defined with a negative `angle_theta` which would lead to wrong injection.
+- Plots of objects defined by shape intersection logic will no longer display thin line artifacts.
 
 ### Changed
+- Switched to an analytical gradient calculation for spatially-varying pole-residue models (`CustomPoleResidue`).
+- Significantly improved performance of the `tidy3d.plugins.autograd.grey_dilation` morphological operation and its gradient calculation. The new implementation is orders of magnitude faster, especially for large arrays and kernel sizes.
 - `GaussianBeam` and `AstigmaticGaussianBeam` default `num_freqs` reset to 1 (it was set to 3 in v2.8.0) and a warning is issued for a broadband, angled beam for which `num_freqs` may not be sufficiently large.
 - Set the maximum `num_freqs` to 20 for all broadband sources (we have been warning about the introduction of this hard limit for a while).
 

tests/test_components/test_geometry.py

@@ -14,6 +14,8 @@
 import trimesh
 
 import tidy3d as td
+from tidy3d.compat import _shapely_is_older_than
+from tidy3d.components.geometry.base import cleanup_shapely_object
 from tidy3d.components.geometry.mesh import AREA_SIZE_THRESHOLD
 from tidy3d.components.geometry.utils import (
     SnapBehavior,
@@ -1249,3 +1251,50 @@ def wrong_shape_height_func(x, y):
     error_message = str(excinfo.value)
     assert f"shape {expected_shape}" in error_message
     assert "shape (3, 3)" in error_message
+
+
+def test_cleanup_shapely_object():
+    if _shapely_is_older_than("2.1"):
+        # (Old versions of shapely don't support `shapely.make_valid()` with the correct arguments.
+        # However older alternatives like `.buffer(0)` are not as robust.  `.buffer(0)` is likely
+        # to generate polygons which look correct, but have extra vertices, causing test to fail.
+        # So if `shapely.make_valid()` is not supported, the safest thing to do is skip this test.)
+        pytest.skip("This test requires `shapely` version 2.1 or later")
+
+    # Test 1: A square containing a triangular hole, and an infinitley thin hole.
+    square_with_spikes_5x5 = np.array(
+        (
+            (0, 0),
+            (5, 0),
+            (5, 0),
+            (5, 10),  # this vertex creates an outward spike and should be removed
+            (5, 5),
+            (0, 5 - 1e-13),  # this vertex should be rounded to (0, 5)
+            (3, 3),  # this vertex creates an inward spike and should be removed
+            (0, 5 + 1e-13),  # this vertex should be removed because it duplicates (0, 5 - 1e-13)
+        )
+    )
+    triangle_empty_tails = np.array(((1, 1), (3, 1), (2, 2), (2.5, 2.5), (0.5, 0.5)))
+    triangle_collinear = np.array(((4, 2), (3, 3), (2, 4), (4, 2)))  # has zero area
+    # NOTE: Test will fail for self intersecting polys like: ((0,0), (1,1), (2,-1), (3,1), (4,0))
+    # Now build a shapely polygon with the 4 small polygons enclosed by big_square_5x5
+    exterior_coords = square_with_spikes_5x5
+    interior_coords_list = [
+        triangle_empty_tails,
+        triangle_collinear,  # this triangle should be eliminated
+    ]
+    # Test using a non-empty exterior polygon (big_square_5x5)
+    orig_polygon = shapely.Polygon(exterior_coords, interior_coords_list)
+    new_polygon = cleanup_shapely_object(orig_polygon, tolerance_ratio=1e-12)
+    # Delete any nearby or overlapping vertices (cleanup_shapely_object() does not do this).
+    new_polygon = shapely.simplify(new_polygon, tolerance=1e-10)
+    # Now `new_polygon` should only contain the coordinates of the square (with a duplicate at end).
+    assert len(new_polygon.exterior.coords) == 5  # squares have 4 vertices but shapely adds 1
+    assert len(new_polygon.interiors) == 1  # only the "triangle_empty_tails" interior hole survives
+    assert len(new_polygon.interiors[0].coords) == 4  # triangles have 3 vertices but shapely adds 1
+    # Test 2: An infinitely thin triangle exterior with some holes (which should be deleted)
+    exterior_coords = triangle_collinear  # has zero area
+    orig_polygon = shapely.Polygon(exterior_coords)
+    new_polygon = cleanup_shapely_object(orig_polygon, tolerance_ratio=1e-12)
+    new_polygon = shapely.simplify(new_polygon, tolerance=1e-10)
+    assert len(new_polygon.exterior.coords) == 0  # empty / collinear polygons should get deleted

tidy3d/compat.py

@@ -2,9 +2,23 @@
 
 from __future__ import annotations
 
+import importlib
+
+from packaging.version import parse as parse_version
+
 try:
     from xarray.structure import alignment
 except ImportError:
     from xarray.core import alignment
 
+
+_SHAPELY_VERSION = parse_version(importlib.metadata.version("shapely"))
+
+
+def _shapely_is_older_than(version: str) -> bool:
+    if _SHAPELY_VERSION < parse_version(version):
+        return True
+    return False
+
+
 __all__ = ["alignment"]

tidy3d/components/geometry/base.py

@@ -17,6 +17,7 @@
 except ImportError:
     pass
 
+from tidy3d.compat import _shapely_is_older_than
 from tidy3d.components.autograd import AutogradFieldMap, TracedCoordinate, TracedSize, get_static
 from tidy3d.components.autograd.derivative_utils import DerivativeInfo, integrate_within_bounds
 from tidy3d.components.base import Tidy3dBaseModel, cached_property
@@ -63,6 +64,7 @@
 from .bound_ops import bounds_intersection, bounds_union
 
 POLY_GRID_SIZE = 1e-12
+POLY_TOLERANCE_RATIO = 1e-12
 
 
 _shapely_operations = {
@@ -2862,27 +2864,45 @@ def _geometries_untraced(cls, val):
         return val
 
     @staticmethod
-    def to_polygon_list(base_geometry: Shapely) -> list[Shapely]:
+    def to_polygon_list(base_geometry: Shapely, cleanup: bool = False) -> list[Shapely]:
         """Return a list of valid polygons from a shapely geometry, discarding points, lines, and
-        empty polygons.
+        empty polygons, and empty triangles within polygons.
 
         Parameters
         ----------
         base_geometry : shapely.geometry.base.BaseGeometry
             Base geometry for inspection.
+        cleanup: bool = False
+            If True, removes extremely small features from each polygon's boundary.
+            This is useful for removing artifacts from 2D plots displayed to the user.
 
         Returns
         -------
         List[shapely.geometry.base.BaseGeometry]
             Valid polygons retrieved from ``base geometry``.
         """
+        unfiltered_geoms = []
         if base_geometry.geom_type == "GeometryCollection":
-            return [p for geom in base_geometry.geoms for p in ClipOperation.to_polygon_list(geom)]
+            unfiltered_geoms = [
+                p for geom in base_geometry.geoms for p in ClipOperation.to_polygon_list(geom)
+            ]
         if base_geometry.geom_type == "MultiPolygon":
-            return [p for p in base_geometry.geoms if not p.is_empty]
+            unfiltered_geoms = [p for p in base_geometry.geoms if not p.is_empty]
         if base_geometry.geom_type == "Polygon" and not base_geometry.is_empty:
-            return [base_geometry]
-        return []
+            unfiltered_geoms = [base_geometry]
+        geoms = []
+        if cleanup:
+            # Optional: "clean" each of the polygons (by removing extremely small or thin features).
+            for geom in unfiltered_geoms:
+                geom_clean = cleanup_shapely_object(geom)
+                if geom_clean.geom_type == "Polygon":
+                    geoms.append(geom_clean)
+                if geom_clean.geom_type == "MultiPolygon":
+                    geoms += [p for p in geom_clean.geoms if not p.is_empty]
+                # Ignore other types of shapely objects (points and lines)
+        else:
+            geoms = unfiltered_geoms
+        return geoms
 
     @property
     def _shapely_operation(self) -> Callable[[Shapely, Shapely], Shapely]:
@@ -2931,7 +2951,10 @@ def intersections_tilted_plane(
         b = self.geometry_b.intersections_tilted_plane(normal, origin, to_2D)
         geom_a = shapely.unary_union([Geometry.evaluate_inf_shape(g) for g in a])
         geom_b = shapely.unary_union([Geometry.evaluate_inf_shape(g) for g in b])
-        return ClipOperation.to_polygon_list(self._shapely_operation(geom_a, geom_b))
+        return ClipOperation.to_polygon_list(
+            self._shapely_operation(geom_a, geom_b),
+            cleanup=True,
+        )
 
     def intersections_plane(
         self, x: Optional[float] = None, y: Optional[float] = None, z: Optional[float] = None
@@ -2958,7 +2981,10 @@ def intersections_plane(
         b = self.geometry_b.intersections_plane(x, y, z)
         geom_a = shapely.unary_union([Geometry.evaluate_inf_shape(g) for g in a])
         geom_b = shapely.unary_union([Geometry.evaluate_inf_shape(g) for g in b])
-        return ClipOperation.to_polygon_list(self._shapely_operation(geom_a, geom_b))
+        return ClipOperation.to_polygon_list(
+            self._shapely_operation(geom_a, geom_b),
+            cleanup=True,
+        )
 
     @cached_property
     def bounds(self) -> Bound:
@@ -3293,4 +3319,92 @@ def _compute_derivatives(self, derivative_info: DerivativeInfo) -> AutogradField
         return grad_vjps
 
 
+def cleanup_shapely_object(obj: Shapely, tolerance_ratio: float = POLY_TOLERANCE_RATIO) -> Shapely:
+    """Remove small geometric features from the boundaries of a shapely object including
+    inward and outward spikes, thin holes, and thin connections between larger regions.
+
+    Parameters
+    ----------
+    obj : shapely
+        a shapely object (typically a ``Polygon`` or a ``MultiPolygon``)
+    tolerance_ratio : float = ``POLY_TOLERANCE_RATIO``
+        Features on the boundaries of polygons will be discarded if they are smaller
+        or narrower than ``tolerance_ratio`` multiplied by the size of the object.
+
+    Returns
+    -------
+    Shapely
+        A new shapely object whose small features (eg. thin spikes or holes) are removed.
+
+    Notes
+    -----
+    This function does not attempt to delete overlapping, nearby, or collinear vertices.
+    To solve that problem, use ``shapely.simplify()`` afterwards.
+    """
+    if _shapely_is_older_than("2.1"):
+        log.warning(
+            "Using old versions of the shapely library (prior to v2.1) may cause "
+            "plot errors.  This can be solved by upgrading to Python 3.10 "
+            "(or later) and reinstalling Tidy3d.",
+            log_once=True,
+        )
+        return obj
+    if obj.is_empty:
+        return obj
+    centroid = obj.centroid
+    object_size = min(obj.bounds[2] - obj.bounds[0], obj.bounds[3] - obj.bounds[1])
+    if object_size == 0.0:
+        return shapely.Polygon([])
+    # In order to prevent numerical overflow or underflow errors, we first subtract
+    # the centroid and divide by (rescale) the size of the object so it is not too big.
+    normalized_obj = shapely.affinity.affine_transform(
+        # https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations
+        obj,
+        matrix=[
+            1 / object_size,
+            0.0,
+            0.0,
+            1 / object_size,
+            -centroid.x / object_size,
+            -centroid.y / object_size,
+        ],
+    )
+    # Important: Remove any self intersections beforehand using `shapely.make_valid()`.
+    valid_obj = shapely.make_valid(normalized_obj, method="structure", keep_collapsed=False)
+    # To get rid of small thin features, erode(shrink), dilate(expand), and erode again.
+    eroded_obj = shapely.buffer(  # This removes outward spikes
+        valid_obj,
+        distance=-tolerance_ratio,
+        cap_style="square",  # (optional parameter to reduce computation time)
+        quad_segs=3,  # (optional parameter to reduce computation time)
+    )
+    dilated_obj = shapely.buffer(  # This removes inward spikes and tiny holes
+        eroded_obj,
+        distance=2 * tolerance_ratio,
+        cap_style="square",
+        quad_segs=3,
+    )
+    cleaned_obj = dilated_obj
+    # Optional: Now shrink the polygon back to the original size.
+    cleaned_obj = shapely.buffer(
+        cleaned_obj,
+        distance=-tolerance_ratio,
+        cap_style="square",
+        quad_segs=3,
+    )
+    # Revert to the original scale and position.
+    rescaled_clean_obj = shapely.affinity.affine_transform(
+        cleaned_obj,
+        matrix=[
+            object_size,
+            0.0,
+            0.0,
+            object_size,
+            centroid.x,
+            centroid.y,
+        ],
+    )
+    return rescaled_clean_obj
+
+
 from .utils import GeometryType, from_shapely, vertices_from_shapely  # noqa: E402

tidy3d/components/viz/descartes.py

@@ -35,7 +35,10 @@ def __init__(self, context):
     @property
     def exterior(self):
         """Get polygon exterior."""
-        return getattr(self.context, "exterior", None) or self.context[0]
+        value = getattr(self.context, "exterior", None)
+        if value is None:
+            value = self.context[0]
+        return value
 
     @property
     def interiors(self):
@@ -53,9 +56,13 @@ def polygon_path(polygon):
     def coding(obj):
         # The codes will be all "LINETO" commands, except for "MOVETO"s at the
         # beginning of each subpath
-        n = len(getattr(obj, "coords", None) or obj)
+        crds = getattr(obj, "coords", None)
+        if crds is None:
+            crds = obj
+        n = len(crds)
         vals = ones(n, dtype=Path.code_type) * Path.LINETO
-        vals[0] = Path.MOVETO
+        if len(vals) > 0:
+            vals[0] = Path.MOVETO
         return vals
 
     ptype = polygon.geom_type