feat(meshing): disable polygon cleanup in corner finder

Corner finder: reduce the number of vertices for smoothly curved cross section

Mesher and gap refinement: performance improvement in parse_structures

Author: weiliangjin2021

tests/test_components/test_geometry.py

@@ -236,6 +236,35 @@ def test_intersections_plane_inf():
     assert len(c.intersections_plane(y=0)) == 1
 
 
+@pytest.mark.parametrize("component", [BOX, CYLINDER, SPHERE, POLYSLAB, UNION, GROUP])
+@pytest.mark.parametrize("cleanup", [True, False])
+def test_intersections_plane_cleanup_param(component, cleanup):
+    """Test that cleanup parameter is accepted by all geometry types."""
+    shapes = component.intersections_plane(z=0, cleanup=cleanup)
+    assert isinstance(shapes, list)
+    shapes_default = component.intersections_plane(z=0)
+    # Both should return valid shapely objects
+    for shape in shapes:
+        assert hasattr(shape, "is_valid")
+    for shape in shapes_default:
+        assert hasattr(shape, "is_valid")
+
+
+@pytest.mark.parametrize("component", [SPHERE, CYLINDER])
+@pytest.mark.parametrize("quad_segs", [8, 50, 200])
+def test_intersections_plane_quad_segs(component, quad_segs):
+    """Test that quad_segs parameter controls discretization of circular shapes."""
+    shapes = component.intersections_plane(z=0, quad_segs=quad_segs)
+    assert len(shapes) > 0
+    # For circular shapes, quad_segs should affect the number of vertices
+    # Higher quad_segs should generally give more vertices
+    shape = shapes[0]
+    num_coords = len(shape.exterior.coords)
+    assert num_coords > 4 * quad_segs, (
+        f"Expected more than {4 * quad_segs} coords, got {num_coords}"
+    )
+
+
 def test_center_not_inf_validate():
     with pytest.raises(pydantic.ValidationError):
         _ = td.Box(center=(td.inf, 0, 0))

tidy3d/components/geometry/base.py

@@ -530,7 +530,12 @@ def bounds(self) -> Bound:
         return self.trimesh.bounds
 
     def intersections_tilted_plane(
-        self, normal: Coordinate, origin: Coordinate, to_2D: MatrixReal4x4
+        self,
+        normal: Coordinate,
+        origin: Coordinate,
+        to_2D: MatrixReal4x4,
+        cleanup: bool = True,
+        quad_segs: Optional[int] = None,
     ) -> list[Shapely]:
         """Return a list of shapely geometries at the plane specified by normal and origin.
 
@@ -542,6 +547,10 @@ def intersections_tilted_plane(
             Vector defining the plane origin.
         to_2D : MatrixReal4x4
             Transformation matrix to apply to resulting shapes.
+        cleanup : bool = True
+            If True, removes extremely small features from each polygon's boundary.
+        quad_segs : Optional[int] = None
+            Number of segments used to discretize circular shapes. Not used for TriangleMesh.
 
         Returns
         -------
@@ -557,7 +566,12 @@ def intersections_tilted_plane(
         return path.polygons_full
 
     def intersections_plane(
-        self, x: Optional[float] = None, y: Optional[float] = None, z: Optional[float] = None
+        self,
+        x: Optional[float] = None,
+        y: Optional[float] = None,
+        z: Optional[float] = None,
+        cleanup: bool = True,
+        quad_segs: Optional[int] = None,
     ) -> list[Shapely]:
         """Returns list of shapely geometries at plane specified by one non-None value of x,y,z.
 
@@ -569,6 +583,10 @@ def intersections_plane(
             Position of plane in y direction, only one of x,y,z can be specified to define plane.
         z : float = None
             Position of plane in z direction, only one of x,y,z can be specified to define plane.
+        cleanup : bool = True
+            If True, removes extremely small features from each polygon's boundary.
+        quad_segs : Optional[int] = None
+            Number of segments used to discretize circular shapes. Not used for TriangleMesh.
 
         Returns
         -------
@@ -623,7 +641,7 @@ def intersections_plane(
                     "Using bounding box instead."
                 )
             log.warning(f"Error encountered: {e}")
-            return self.bounding_box.intersections_plane(x=x, y=y, z=z)
+            return self.bounding_box.intersections_plane(x=x, y=y, z=z, cleanup=cleanup)
 
     def inside(
         self, x: np.ndarray[float], y: np.ndarray[float], z: np.ndarray[float]

tidy3d/components/geometry/mesh.py

@@ -587,7 +587,11 @@ def _move_axis_reverse(arr: NDArray) -> NDArray:
 
     @verify_packages_import(["trimesh"])
     def _do_intersections_tilted_plane(
-        self, normal: Coordinate, origin: Coordinate, to_2D: MatrixReal4x4
+        self,
+        normal: Coordinate,
+        origin: Coordinate,
+        to_2D: MatrixReal4x4,
+        quad_segs: Optional[int] = None,
     ) -> list[Shapely]:
         """Return a list of shapely geometries at the plane specified by normal and origin.
 
@@ -599,6 +603,8 @@ def _do_intersections_tilted_plane(
             Vector defining the plane origin.
         to_2D : MatrixReal4x4
             Transformation matrix to apply to resulting shapes.
+        quad_segs : Optional[int] = None
+            Number of segments used to discretize circular shapes. Not used for PolySlab geometry.
 
         Returns
         -------
@@ -641,7 +647,7 @@ def _do_intersections_tilted_plane(
         path, _ = section.to_2D(to_2D=to_2D)
         return path.polygons_full
 
-    def _intersections_normal(self, z: float) -> list[Shapely]:
+    def _intersections_normal(self, z: float, quad_segs: Optional[int] = None) -> list[Shapely]:
         """Find shapely geometries intersecting planar geometry with axis normal to slab.
 
         Parameters
@@ -2567,7 +2573,12 @@ def _dilation_value_at_reference_to_coord(self, dilation: float) -> float:
         return z_coord
 
     def intersections_tilted_plane(
-        self, normal: Coordinate, origin: Coordinate, to_2D: MatrixReal4x4
+        self,
+        normal: Coordinate,
+        origin: Coordinate,
+        to_2D: MatrixReal4x4,
+        cleanup: bool = True,
+        quad_segs: Optional[int] = None,
     ) -> list[Shapely]:
         """Return a list of shapely geometries at the plane specified by normal and origin.
 
@@ -2579,6 +2590,10 @@ def intersections_tilted_plane(
             Vector defining the plane origin.
         to_2D : MatrixReal4x4
             Transformation matrix to apply to resulting shapes.
+        cleanup : bool = True
+            If True, removes extremely small features from each polygon's boundary.
+        quad_segs : Optional[int] = None
+            Number of segments used to discretize circular shapes. Not used for PolySlab.
 
         Returns
         -------
@@ -2592,7 +2607,9 @@ def intersections_tilted_plane(
                 [
                     base.Geometry.evaluate_inf_shape(shape)
                     for polyslab in self.sub_polyslabs
-                    for shape in polyslab.intersections_tilted_plane(normal, origin, to_2D)
+                    for shape in polyslab.intersections_tilted_plane(
+                        normal, origin, to_2D, cleanup=cleanup, quad_segs=quad_segs
+                    )
                 ]
             )
         ]

tidy3d/components/geometry/polyslab.py

@@ -28,7 +28,7 @@
 _N_SAMPLE_CURVE_SHAPELY = 40
 
 # for shapely circular shapes discretization in visualization
-_N_SHAPELY_QUAD_SEGS = 200
+_N_SHAPELY_QUAD_SEGS_VISUALIZATION = 200
 
 # Default number of points to discretize polyslab in `Cylinder.to_polyslab()`
 _N_PTS_CYLINDER_POLYSLAB = 51
@@ -71,7 +71,12 @@ def inside(
         return (dist_x**2 + dist_y**2 + dist_z**2) <= (self.radius**2)
 
     def intersections_tilted_plane(
-        self, normal: Coordinate, origin: Coordinate, to_2D: MatrixReal4x4
+        self,
+        normal: Coordinate,
+        origin: Coordinate,
+        to_2D: MatrixReal4x4,
+        cleanup: bool = True,
+        quad_segs: Optional[int] = None,
     ) -> list[Shapely]:
         """Return a list of shapely geometries at the plane specified by normal and origin.
 
@@ -83,6 +88,11 @@ def intersections_tilted_plane(
             Vector defining the plane origin.
         to_2D : MatrixReal4x4
             Transformation matrix to apply to resulting shapes.
+        cleanup : bool = True
+            If True, removes extremely small features from each polygon's boundary.
+        quad_segs : Optional[int] = None
+            Number of segments used to discretize circular shapes. If ``None``, uses
+            ``_N_SHAPELY_QUAD_SEGS_VISUALIZATION`` for high-quality visualization.
 
         Returns
         -------
@@ -91,6 +101,9 @@ def intersections_tilted_plane(
             For more details refer to
             `Shapely's Documentation <https://shapely.readthedocs.io/en/stable/project.html>`_.
         """
+        if quad_segs is None:
+            quad_segs = _N_SHAPELY_QUAD_SEGS_VISUALIZATION
+
         normal = np.array(normal)
         unit_normal = normal / (np.sum(normal**2) ** 0.5)
         projection = np.dot(np.array(origin) - np.array(self.center), unit_normal)
@@ -106,13 +119,18 @@ def intersections_tilted_plane(
         u /= np.sum(u**2) ** 0.5
         v = np.cross(unit_normal, u)
 
-        angles = np.linspace(0, 2 * np.pi, _N_SHAPELY_QUAD_SEGS * 4 + 1)[:-1]
+        angles = np.linspace(0, 2 * np.pi, quad_segs * 4 + 1)[:-1]
         circ = center + np.outer(np.cos(angles), radius * u) + np.outer(np.sin(angles), radius * v)
         vertices = np.dot(np.hstack((circ, np.ones((angles.size, 1)))), to_2D.T)
         return [shapely.Polygon(vertices[:, :2])]
 
     def intersections_plane(
-        self, x: Optional[float] = None, y: Optional[float] = None, z: Optional[float] = None
+        self,
+        x: Optional[float] = None,
+        y: Optional[float] = None,
+        z: Optional[float] = None,
+        cleanup: bool = True,
+        quad_segs: Optional[int] = None,
     ) -> list[BaseGeometry]:
         """Returns shapely geometry at plane specified by one non None value of x,y,z.
 
@@ -124,22 +142,30 @@ def intersections_plane(
             Position of plane in x direction, only one of x,y,z can be specified to define plane.
         z : float = None
             Position of plane in x direction, only one of x,y,z can be specified to define plane.
+        cleanup : bool = True
+            If True, removes extremely small features from each polygon's boundary.
+        quad_segs : Optional[int] = None
+            Number of segments used to discretize circular shapes. If ``None``, uses
+            ``_N_SHAPELY_QUAD_SEGS_VISUALIZATION`` for high-quality visualization.
 
         Returns
         -------
         List[shapely.geometry.base.BaseGeometry]
             List of 2D shapes that intersect plane.
             For more details refer to
-            `Shapely's Documentation <https://shapely.readthedocs.io/en/stable/project.html>`_.
+            `Shapely's Documentation <https://shapely.readthedocs.io/en/stable/project.html>``.
         """
+        if quad_segs is None:
+            quad_segs = _N_SHAPELY_QUAD_SEGS_VISUALIZATION
+
         axis, position = self.parse_xyz_kwargs(x=x, y=y, z=z)
         if not self.intersects_axis_position(axis, position):
             return []
         z0, (x0, y0) = self.pop_axis(self.center, axis=axis)
         intersect_dist = self._intersect_dist(position, z0)
         if not intersect_dist:
             return []
-        return [shapely.Point(x0, y0).buffer(0.5 * intersect_dist, quad_segs=_N_SHAPELY_QUAD_SEGS)]
+        return [shapely.Point(x0, y0).buffer(0.5 * intersect_dist, quad_segs=quad_segs)]
 
     @cached_property
     def bounds(self) -> Bound:
@@ -430,7 +456,11 @@ def _update_from_bounds(self, bounds: tuple[float, float], axis: Axis) -> Cylind
 
     @verify_packages_import(["trimesh"])
     def _do_intersections_tilted_plane(
-        self, normal: Coordinate, origin: Coordinate, to_2D: MatrixReal4x4
+        self,
+        normal: Coordinate,
+        origin: Coordinate,
+        to_2D: MatrixReal4x4,
+        quad_segs: Optional[int] = None,
     ) -> list[Shapely]:
         """Return a list of shapely geometries at the plane specified by normal and origin.
 
@@ -442,6 +472,9 @@ def _do_intersections_tilted_plane(
             Vector defining the plane origin.
         to_2D : MatrixReal4x4
             Transformation matrix to apply to resulting shapes.
+        quad_segs : Optional[int] = None
+            Number of segments used to discretize circular shapes. If ``None``, uses
+            ``_N_SHAPELY_QUAD_SEGS_VISUALIZATION`` for high-quality visualization.
 
         Returns
         -------
@@ -452,6 +485,9 @@ def _do_intersections_tilted_plane(
         """
         import trimesh
 
+        if quad_segs is None:
+            quad_segs = _N_SHAPELY_QUAD_SEGS_VISUALIZATION
+
         z0, (x0, y0) = self.pop_axis(self.center, self.axis)
         half_length = self.finite_length_axis / 2
 
@@ -471,7 +507,7 @@ def _do_intersections_tilted_plane(
                 r_bot = 0
                 z_bot = z0 + self._radius_z(z0) / self._tanq
 
-        angles = np.linspace(0, 2 * np.pi, _N_SHAPELY_QUAD_SEGS * 4 + 1)
+        angles = np.linspace(0, 2 * np.pi, quad_segs * 4 + 1)
 
         if r_bot > 0:
             x_bot = x0 + r_bot * np.cos(angles)
@@ -525,13 +561,18 @@ def _do_intersections_tilted_plane(
         path, _ = section.to_2D(to_2D=to_2D)
         return path.polygons_full
 
-    def _intersections_normal(self, z: float) -> list[BaseGeometry]:
+    def _intersections_normal(
+        self, z: float, quad_segs: Optional[int] = None
+    ) -> list[BaseGeometry]:
         """Find shapely geometries intersecting cylindrical geometry with axis normal to slab.
 
         Parameters
         ----------
         z : float
             Position along the axis normal to slab
+        quad_segs : Optional[int] = None
+            Number of segments used to discretize circular shapes. If ``None``, uses
+            ``_N_SHAPELY_QUAD_SEGS_VISUALIZATION`` for high-quality visualization.
 
         Returns
         -------
@@ -540,6 +581,8 @@ def _intersections_normal(self, z: float) -> list[BaseGeometry]:
             For more details refer to
             `Shapely's Documentation <https://shapely.readthedocs.io/en/stable/project.html>`_.
         """
+        if quad_segs is None:
+            quad_segs = _N_SHAPELY_QUAD_SEGS_VISUALIZATION
 
         static_self = self.to_static()
 
@@ -550,7 +593,7 @@ def _intersections_normal(self, z: float) -> list[BaseGeometry]:
             return []
 
         _, (x0, y0) = self.pop_axis(static_self.center, axis=self.axis)
-        return [shapely.Point(x0, y0).buffer(radius_offset, quad_segs=_N_SHAPELY_QUAD_SEGS)]
+        return [shapely.Point(x0, y0).buffer(radius_offset, quad_segs=quad_segs)]
 
     def _intersections_side(self, position: float, axis: int) -> list[BaseGeometry]:
         """Find shapely geometries intersecting cylindrical geometry with axis orthogonal to length.

tidy3d/components/geometry/primitives.py

@@ -93,6 +93,8 @@ def merging_geometries_on_plane(
     plane: Box,
     property_list: list[Any],
     interior_disjoint_geometries: bool = False,
+    cleanup: bool = True,
+    quad_segs: Optional[int] = None,
 ) -> list[tuple[Any, Shapely]]:
     """Compute list of shapes on plane. Overlaps are removed or merged depending on
     provided property_list.
@@ -107,6 +109,11 @@ def merging_geometries_on_plane(
         Property value for each structure.
     interior_disjoint_geometries: bool = False
         If ``True``, geometries of different properties on the plane must not be overlapping.
+    cleanup : bool = True
+        If True, removes extremely small features from each polygon's boundary.
+    quad_segs : Optional[int] = None
+        Number of segments used to discretize circular shapes. If ``None``, uses
+        high-quality visualization settings.
 
     Returns
     -------
@@ -122,7 +129,7 @@ def merging_geometries_on_plane(
     shapes = []
     for geo, prop in zip(geometries, property_list):
         # get list of Shapely shapes that intersect at the plane
-        shapes_plane = plane.intersections_with(geo)
+        shapes_plane = plane.intersections_with(geo, cleanup=cleanup, quad_segs=quad_segs)
 
         # Append each of them and their property information to the list of shapes
         for shape in shapes_plane:

tidy3d/components/geometry/utils.py

@@ -16,6 +16,10 @@
 from tidy3d.constants import inf
 
 CORNER_ANGLE_THRESOLD = 0.25 * np.pi
+# For shapely circular shapes discretization.
+N_SHAPELY_QUAD_SEGS = 8
+# whether to clean tiny features that sometimes occurs in shapely operations
+SHAPELY_CLEANUP = False
 
 
 class CornerFinderSpec(Tidy3dBaseModel):
@@ -132,7 +136,12 @@ def _merged_pec_on_plane(
             medium_list = [PEC for _ in geometry_list]
         # merge geometries
         merged_geos = merging_geometries_on_plane(
-            geometry_list, plane, medium_list, interior_disjoint_geometries
+            geometry_list,
+            plane,
+            medium_list,
+            interior_disjoint_geometries,
+            cleanup=SHAPELY_CLEANUP,
+            quad_segs=N_SHAPELY_QUAD_SEGS,
         )
 
         return merged_geos