From a97af623cc2caa8d0ebdb3d7e1b8e1fe158caeba Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Mon, 1 Sep 2025 21:47:51 -0700 Subject: [PATCH 1/3] Update SSP with erosions/dilations Here we update the SSP algorithm to handle erosion and dilations functions. --- python/adjoint/filters.py | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/python/adjoint/filters.py b/python/adjoint/filters.py index 09d2bef01..10b70238d 100644 --- a/python/adjoint/filters.py +++ b/python/adjoint/filters.py @@ -5,7 +5,7 @@ """ import sys -from typing import List, Tuple, Union +from typing import List, Optional, Tuple, Union import numpy as np from autograd import numpy as npa @@ -712,6 +712,9 @@ def smoothed_projection( beta: float, eta: float, resolution: float, + resolution_simulation: Optional[float] = None, + erosion_dilation: float = 0.0, + ): """Project using subpixel smoothing, which allows for β→∞. @@ -747,6 +750,10 @@ def smoothed_projection( eta: The threshold point in the range [0, 1]. resolution: resolution of the design grid (not the Meep grid resolution). + resolution_simulation: resolution of the meep grid. + erosion_dilation: erosion (negative) or dilation (positive) transform to + apply to the final projected design in "meep units"; The resolution + scaling is applied automatically. Returns: The projected and smoothed output. @@ -763,6 +770,9 @@ def smoothed_projection( >>> rho_filtered = conic_filter(rho, filter_radius, Lx, Ly, resolution) >>> rho_projected = smoothed_projection(rho_filtered, beta, eta_i, resolution) """ + if resolution_simulation is None: + resolution_simulation = resolution + # Note that currently, the underlying assumption is that the smoothing # kernel is a circle, which means dx = dy. dx = dy = 1 / resolution @@ -795,10 +805,13 @@ def smoothed_projection( # The distance for the center of the pixel to the nearest interface d = (eta - rho_filtered) / rho_filtered_grad_norm_eff + # Account for an erosion or dilation + d = npa.maximum(0, d) + # Only need smoothing if an interface lies within the voxel. Since d is # actually an "effective" d by this point, we need to ignore values that may # have been sanitized earlier on. - needs_smoothing = nonzero_norm & (npa.abs(d) < R_smoothing) + needs_smoothing = nonzero_norm & (npa.abs(d + erosion_dilation) < R_smoothing) # The fill factor is used to perform simple, first-order subpixel smoothing. # We use the (2D) analytic expression that comes when assuming the smoothing From 1bc3452b9423330057c4ce47ed192f684223397b Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Tue, 2 Sep 2025 21:22:46 -0700 Subject: [PATCH 2/3] Simplify the erosion and dilation perturbation --- python/adjoint/filters.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/adjoint/filters.py b/python/adjoint/filters.py index 10b70238d..f52e0f25d 100644 --- a/python/adjoint/filters.py +++ b/python/adjoint/filters.py @@ -806,12 +806,12 @@ def smoothed_projection( d = (eta - rho_filtered) / rho_filtered_grad_norm_eff # Account for an erosion or dilation - d = npa.maximum(0, d) + d = d + erosion_dilation # Only need smoothing if an interface lies within the voxel. Since d is # actually an "effective" d by this point, we need to ignore values that may # have been sanitized earlier on. - needs_smoothing = nonzero_norm & (npa.abs(d + erosion_dilation) < R_smoothing) + needs_smoothing = nonzero_norm & (npa.abs(d) < R_smoothing) # The fill factor is used to perform simple, first-order subpixel smoothing. # We use the (2D) analytic expression that comes when assuming the smoothing From 03614f137ec2bb88dcfe87226cecf4a3a8a7c124 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Thu, 4 Sep 2025 20:24:06 -0700 Subject: [PATCH 3/3] Modify the filtered field to include the erosion/dilation --- python/adjoint/filters.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/python/adjoint/filters.py b/python/adjoint/filters.py index f52e0f25d..ed5d5e032 100644 --- a/python/adjoint/filters.py +++ b/python/adjoint/filters.py @@ -778,8 +778,6 @@ def smoothed_projection( dx = dy = 1 / resolution R_smoothing = 0.55 * dx - rho_projected = tanh_projection(rho_filtered, beta=beta, eta=eta) - # Compute the spatial gradient (using finite differences) of the *filtered* # field, which will always be smooth and is the key to our approach. This # gradient essentially represents the normal direction pointing the the @@ -802,12 +800,14 @@ def smoothed_projection( ) rho_filtered_grad_norm_eff = npa.where(nonzero_norm, rho_filtered_grad_norm, 1) + # Account for an erosion or dilation + rho_filtered = rho_filtered - erosion_dilation * npa.where(nonzero_norm, rho_filtered_grad_norm, 0) + + rho_projected = tanh_projection(rho_filtered, beta=beta, eta=eta) + # The distance for the center of the pixel to the nearest interface d = (eta - rho_filtered) / rho_filtered_grad_norm_eff - # Account for an erosion or dilation - d = d + erosion_dilation - # Only need smoothing if an interface lies within the voxel. Since d is # actually an "effective" d by this point, we need to ignore values that may # have been sanitized earlier on.