Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
29 changes: 29 additions & 0 deletions tests/test_components/test_geometry.py
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down
200 changes: 171 additions & 29 deletions tidy3d/components/geometry/base.py

Large diffs are not rendered by default.

24 changes: 21 additions & 3 deletions tidy3d/components/geometry/mesh.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand All @@ -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
-------
Expand All @@ -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.

Expand All @@ -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
-------
Expand Down Expand Up @@ -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]
Expand Down
25 changes: 21 additions & 4 deletions tidy3d/components/geometry/polyslab.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand All @@ -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
-------
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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.

Expand All @@ -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
-------
Expand All @@ -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
)
]
)
]
63 changes: 53 additions & 10 deletions tidy3d/components/geometry/primitives.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.

Expand All @@ -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
-------
Expand All @@ -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)
Expand All @@ -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.

Expand All @@ -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:
Expand Down Expand Up @@ -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.

Expand All @@ -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
-------
Expand All @@ -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

Expand All @@ -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)
Expand Down Expand Up @@ -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
-------
Expand All @@ -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()

Expand All @@ -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.
Expand Down
9 changes: 8 additions & 1 deletion tidy3d/components/geometry/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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
-------
Expand All @@ -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:
Expand Down
11 changes: 10 additions & 1 deletion tidy3d/components/grid/corner_finder.py
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down Expand Up @@ -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
Expand Down
Loading