feat(io): added pixel_exact option to gds export

Author: mahlau-flex

CHANGELOG.md

@@ -13,6 +13,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Added autograd support for `TriangleMesh`, allowing gradient computation with respect to mesh vertices for inverse design.
 - Added `EMESimulation.store_coeffs` to store coefficients from the EME solver, including mode overlaps, interface S matrices, and effective propagation indices.
 - Get cell-related information from violation markers in `DRCResults` and `DRCViolation` to the klayout plugin: Use for example `DRCResults.violations_by_cell` to group them.
+- `pixel_exact` option for `to_gds` export of simulations and structures. If this option is set, any custom medium will be exported as rectangular pixels according to the specification of the medium coordinates.
 
 ### Changed
 - Removed validator that would warn if `PerturbationMedium` values could become numerically unstable, since an error will anyway be raised if this actually happens when the medium is converted using actual perturbation data.

tests/test_components/test_structure.py

@@ -264,3 +264,27 @@ def f(x):
 
     grad = ag.grad(f)(1.0)
     assert not np.isclose(grad, 0.0)
+
+
+def test_to_gdstk_pixel_exact(tmp_path):
+    box = td.Box(center=(0, 0, 0), size=(2, 2, 2))
+    nx, ny = 50, 50
+    arr = np.ones((nx, ny, 1))
+    arr[20:30, 20:30] = 2
+    x = np.linspace(-1, 1, nx)
+    y = np.linspace(-1, 1, ny)
+    z = np.asarray([0])
+    coords = {"x": x, "y": y, "z": z}
+    permittivity = td.SpatialDataArray(arr, coords=coords)
+
+    medium = td.CustomMedium(permittivity=permittivity)
+    structure = td.Structure(geometry=box, medium=medium)
+    polygons = structure.to_gdstk(z=0, frequency=3e14, permittivity_threshold=1.5, pixel_exact=True)
+    assert polygons
+
+    fname = str(tmp_path / "structure-exact.gds")
+    structure.to_gds_file(fname, z=0, permittivity_threshold=1.5, frequency=3e14, pixel_exact=True)
+    cells = gdstk.read_gds(fname).cells
+    cell = cells[0]
+    assert np.allclose(cell.bounding_box(), ((-0.2, -0.2), (0.2, 0.2)), atol=0.01)
+    assert gdstk.inside([(0.1, 0.9), (0.5, 2.5), (0, 0)], cell.polygons) == (False, False, True)

tidy3d/components/simulation.py

@@ -5145,6 +5145,7 @@ def to_gdstk(
         gds_layer_dtype_map: Optional[
             dict[AbstractMedium, tuple[pydantic.NonNegativeInt, pydantic.NonNegativeInt]]
         ] = None,
+        pixel_exact: bool = False,
     ) -> list:
         """Convert a simulation's planar slice to a .gds type polygon list.
 
@@ -5163,6 +5164,8 @@ def to_gdstk(
             Frequency for permittivity evaluation in case of custom medium (Hz).
         gds_layer_dtype_map : Dict
             Dictionary mapping mediums to GDSII layer and data type tuples.
+        pixel_exact : bool = False
+            If true export gds as pixel exact rectangles instead of gdstk contour if a custom medium is provided.
 
         Return
         ------
@@ -5194,6 +5197,7 @@ def to_gdstk(
                 frequency=frequency,
                 gds_layer=gds_layer,
                 gds_dtype=gds_dtype,
+                pixel_exact=pixel_exact,
             ):
                 pmin, pmax = polygon.bounding_box()
                 if pmin[0] < bmin[0] or pmin[1] < bmin[1] or pmax[0] > bmax[0] or pmax[1] > bmax[1]:
@@ -5216,6 +5220,7 @@ def to_gds(
         gds_layer_dtype_map: Optional[
             dict[AbstractMedium, tuple[pydantic.NonNegativeInt, pydantic.NonNegativeInt]]
         ] = None,
+        pixel_exact: bool = False,
     ) -> None:
         """Append the simulation structures to a .gds cell.
 
@@ -5236,6 +5241,8 @@ def to_gds(
             Frequency for permittivity evaluation in case of custom medium (Hz).
         gds_layer_dtype_map : Dict
             Dictionary mapping mediums to GDSII layer and data type tuples.
+        pixel_exact : bool = False
+            If true export gds as pixel exact rectangles instead of gdstk contour if a custom medium is provided.
         """
         if gds_layer_dtype_map is None:
             gds_layer_dtype_map = {}
@@ -5248,6 +5255,7 @@ def to_gds(
                 permittivity_threshold=permittivity_threshold,
                 frequency=frequency,
                 gds_layer_dtype_map=gds_layer_dtype_map,
+                pixel_exact=pixel_exact,
             )
             if len(polygons) > 0:
                 cell.add(*polygons)
@@ -5271,6 +5279,7 @@ def to_gds_file(
             dict[AbstractMedium, tuple[pydantic.NonNegativeInt, pydantic.NonNegativeInt]]
         ] = None,
         gds_cell_name: str = "MAIN",
+        pixel_exact: bool = False,
     ) -> None:
         """Append the simulation structures to a .gds cell.
 
@@ -5293,6 +5302,8 @@ def to_gds_file(
             Dictionary mapping mediums to GDSII layer and data type tuples.
         gds_cell_name : str = 'MAIN'
             Name of the cell created in the .gds file to store the geometry.
+        pixel_exact : bool = False
+            If true export gds as pixel exact rectangles instead of gdstk contour if a custom medium is provided.
         """
         if gdstk_available:
             library = gdstk.Library()
@@ -5326,6 +5337,7 @@ def to_gds_file(
             permittivity_threshold=permittivity_threshold,
             frequency=frequency,
             gds_layer_dtype_map=gds_layer_dtype_map,
+            pixel_exact=pixel_exact,
         )
         fname = pathlib.Path(fname)
         fname.parent.mkdir(parents=True, exist_ok=True)

tidy3d/components/structure.py

@@ -424,6 +424,7 @@ def to_gdstk(
         frequency: pydantic.PositiveFloat = 0,
         gds_layer: pydantic.NonNegativeInt = 0,
         gds_dtype: pydantic.NonNegativeInt = 0,
+        pixel_exact: bool = False,
     ) -> None:
         """Convert a structure's planar slice to a .gds type polygon.
 
@@ -435,15 +436,16 @@ def to_gdstk(
             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.
-        permittivity_threshold : float = 1.1
-            Permitivitty value used to define the shape boundaries for structures with custom
-            medim
+        permittivity_threshold : float = 1
+            Permitivitty value used to define the shape boundaries for structures with custom medium
         frequency : float = 0
             Frequency for permittivity evaluaiton in case of custom medium (Hz).
         gds_layer : int = 0
             Layer index to use for the shapes stored in the .gds file.
         gds_dtype : int = 0
             Data-type index to use for the shapes stored in the .gds file.
+        pixel_exact : bool = False
+            If true export gds as pixel exact rectangles instead of gdstk contour if a custom medium is provided.
 
         Return
         ------
@@ -457,27 +459,72 @@ def to_gdstk(
             axis, _ = self.geometry.parse_xyz_kwargs(x=x, y=y, z=z)
             bb_min, bb_max = self.geometry.bounds
 
-            # Set the contour scale to be the minimal cooridante step size w.r.t. the 3 main axes,
-            # skipping those with a single coordniate. In case all axes have only a single coordinate,
-            # use the largest bounding box dimension.
             eps, _, _ = self.medium.eps_dataarray_freq(frequency=frequency)
-            scale = max(b - a for a, b in zip(bb_min, bb_max))
-            for coord in (eps.x, eps.y, eps.z):
-                if len(coord) > 1:
-                    scale = min(scale, np.diff(coord).min())
-
-            coords = Coords(
-                x=np.arange(bb_min[0], bb_max[0] + scale * 0.9, scale) if x is None else x,
-                y=np.arange(bb_min[1], bb_max[1] + scale * 0.9, scale) if y is None else y,
-                z=np.arange(bb_min[2], bb_max[2] + scale * 0.9, scale) if z is None else z,
-            )
+            if pixel_exact:
+                coords = Coords(
+                    x=eps.x if x is None else x,
+                    y=eps.y if y is None else y,
+                    z=eps.z if z is None else z,
+                )
+            else:
+                # Set the contour scale to be the minimal cooridante step size w.r.t. the 3 main axes,
+                # skipping those with a single coordniate. In case all axes have only a single coordinate,
+                # use the largest bounding box dimension.
+                scale = max(b - a for a, b in zip(bb_min, bb_max))
+                for coord in (eps.x, eps.y, eps.z):
+                    if len(coord) > 1:
+                        scale = min(scale, np.diff(coord).min())
+                coords = Coords(
+                    x=np.arange(bb_min[0], bb_max[0] + scale * 0.9, scale) if x is None else x,
+                    y=np.arange(bb_min[1], bb_max[1] + scale * 0.9, scale) if y is None else y,
+                    z=np.arange(bb_min[2], bb_max[2] + scale * 0.9, scale) if z is None else z,
+                )
+
             eps = self.medium.eps_diagonal_on_grid(frequency=frequency, coords=coords)
-            eps = np.stack((eps[0].real, eps[1].real, eps[2].real), axis=3).max(axis=3).squeeze()
-            contours = gdstk.contour(eps.T, permittivity_threshold, scale, precision=scale * 1e-3)
+            eps = (
+                np.stack((eps[0].real, eps[1].real, eps[2].real), axis=3)
+                .max(axis=3)
+                .squeeze(axis=axis)
+            )
+
+            if pixel_exact:
+                # Convert coordinates to numpy arrays for efficient processing
+                _, (w, h) = self.geometry.pop_axis((coords.x, coords.y, coords.z), axis)
+                w, h = np.asarray(w), np.asarray(h)
+                _, (wmin, hmin) = self.geometry.pop_axis(bb_min, axis)
+                _, (wmax, hmax) = self.geometry.pop_axis(bb_max, axis)
+
+                # Determine boundaries by taking the midpoint between adjacent coordinates
+                if w.size > 1:
+                    dw = np.diff(w) * 0.5
+                    wb = np.concatenate(([wmin], w[:-1] + dw, [wmax]))
+                else:
+                    wb = np.array([wmin, wmax])
+
+                if h.size > 1:
+                    dh = np.diff(h) * 0.5
+                    hb = np.concatenate(([hmin], h[:-1] + dh, [hmax]))
+                else:
+                    hb = np.array([hmin, hmax])
+
+                # Create boolean mask where permittivity exceeds threshold
+                mask = eps > permittivity_threshold
+                w_idxs, h_idxs = np.where(mask)
+
+                # Generate list of gdstk.Polygon (rectangles)
+                contours = [
+                    gdstk.rectangle((wb[wi], hb[hi]), (wb[wi + 1], hb[hi + 1]))
+                    for wi, hi in zip(w_idxs, h_idxs)
+                ]
+
+            else:
+                contours = gdstk.contour(
+                    eps.T, permittivity_threshold, scale, precision=scale * 1e-3
+                )
 
-            _, (dx, dy) = self.geometry.pop_axis(bb_min, axis)
-            for polygon in contours:
-                polygon.translate(dx, dy)
+                _, (dx, dy) = self.geometry.pop_axis(bb_min, axis)
+                for polygon in contours:
+                    polygon.translate(dx, dy)
 
             polygons = gdstk.boolean(polygons, contours, "and", layer=gds_layer, datatype=gds_dtype)
 
@@ -493,6 +540,7 @@ def to_gds(
         frequency: pydantic.PositiveFloat = 0,
         gds_layer: pydantic.NonNegativeInt = 0,
         gds_dtype: pydantic.NonNegativeInt = 0,
+        pixel_exact: bool = False,
     ) -> None:
         """Append a structure's planar slice to a .gds cell.
 
@@ -515,6 +563,8 @@ def to_gds(
             Layer index to use for the shapes stored in the .gds file.
         gds_dtype : int = 0
             Data-type index to use for the shapes stored in the .gds file.
+        pixel_exact : bool = False
+            If true export gds as pixel exact rectangles instead of gdstk contour if a custom medium is provided.
         """
         if not isinstance(cell, gdstk.Cell):
             if "gdstk" in cell.__class__.__name__.lower() and not gdstk_available:
@@ -531,6 +581,7 @@ def to_gds(
             frequency=frequency,
             gds_layer=gds_layer,
             gds_dtype=gds_dtype,
+            pixel_exact=pixel_exact,
         )
         if polygons:
             cell.add(*polygons)
@@ -546,6 +597,7 @@ def to_gds_file(
         gds_layer: pydantic.NonNegativeInt = 0,
         gds_dtype: pydantic.NonNegativeInt = 0,
         gds_cell_name: str = "MAIN",
+        pixel_exact: bool = False,
     ) -> None:
         """Export a structure's planar slice to a .gds file.
 
@@ -570,6 +622,8 @@ def to_gds_file(
             Data-type index to use for the shapes stored in the .gds file.
         gds_cell_name : str = 'MAIN'
             Name of the cell created in the .gds file to store the geometry.
+        pixel_exact : bool = False
+            If true export gds as pixel exact rectangles instead of gdstk contour if a custom medium is provided.
         """
         try:
             import gdstk
@@ -590,6 +644,7 @@ def to_gds_file(
             frequency=frequency,
             gds_layer=gds_layer,
             gds_dtype=gds_dtype,
+            pixel_exact=pixel_exact,
         )
         fname = pathlib.Path(fname)
         fname.parent.mkdir(parents=True, exist_ok=True)