Update SSP with erosions/dilations

[smartalecH]

Here we update the SSP algorithm to account for the following features/fixes:

• (TODO) Generalize the routine to work in 1D and 3D (not just 2D)
• (in-progress) Add resolution_simulation parameter to account for differences in meep and design grid resolutions
• (in-progress) Augment the algorithm to allow for precise erosions and dilations.

(cc @oskooi, @stevengj )

[smartalecH]

@stevengj I have the following test for our new erosion and dilation feature:

Lx = 2; Ly = 2
resolution = 50
eta_i = 0.5; eta_e = 0.75
lengthscale = 0.1
filter_radius = mpa.get_conic_radius_from_eta_e(lengthscale, eta_e)
Nx = onp.round(Lx * resolution) + 1
Ny = onp.round(Ly * resolution) + 1
rho = onp.random.rand(Nx, Ny)
beta = npa.inf
rho_filtered = mpa.conic_filter(rho, filter_radius, Lx, Ly, resolution)
rho_projected = mpa.smoothed_projection(rho_filtered, beta, eta_i, resolution)
erosion_dilation = 0.5
rho_projected_eroded = mpa.smoothed_projection(rho_filtered, beta, eta_i, resolution, erosion_dilation=-erosion_dilation)
rho_projected_dilated = mpa.smoothed_projection(rho_filtered, beta, eta_i, resolution, erosion_dilation=+erosion_dilation)

plt.figure(figsize=(12,4))
plt.subplot(1,3,1)
plt.imshow(rho_projected_eroded)
plt.subplot(1,3,2)
plt.imshow(rho_projected)
plt.subplot(1,3,3)
plt.imshow(rho_projected_dilated)

It's obviously not quite there yet.

I'm not sure how to tackle the following challenge: we only use our level-set approximation to identify the boundaries themselves. Anything not on a boundary, we fall back to our brute-force tanh_projection output (hence the where statement).

It seems that in order to do erosions/dilations, we really have to use the level-set across the whole domain, right?

[stevengj]

Your code doesn't seem right?

Let's consider the $\beta = \infty$ case first, where the tanh projection is just the step function and SSP boils down to $F(d)$. We should just be doing $F(d + \text{erosion\_dilation})$, no?

[stevengj]

I think this equation in the SSP paper needs to be modified too:

because it is relative to an interface at the wrong (old) position.

However, if you just use $F(d)$ directly for $\beta = \infty$, as I suggested above, you don't need to worry about this for now. That's why I think it will be easier to debug that case first.

My guess is that the corrected general formula should be something like:

$\hat{\rho}_{\pm} = P_{\beta,\eta}\left(\tilde{\rho} \pm (\hat{R} F(\mp d) \mp \Delta d) \, \Vert \nabla \tilde{\rho} \Vert \right)$

where $\Delta d$ is what you are calling erosion_dilation, i.e. you have replaced $d$ everywhere with $d_{\text{new}} = d_{\text{old}} + \Delta d$. … I have to check the continuity conditions on $\hat{\tilde{\rho}}$, though.

[stevengj]

Let me make an alternative suggestion.

1. Move the computation of rho_filtered_grad_norm_eff before the computation of rho_projected.
2. Add the following line before the computation of rho_projected: rho_filtered = rho_filtered - erosion_dilation * npa.where(nonzero_norm, rho_filtered_grad_norm, 0)
3. Remove all of your other changes to the code (e.g. now d will automatically include erosion_dilation)

[smartalecH]

@stevengj ok so that indeed produces an erosion or dilation:

We should do some experiments to show that the resulting transformation really reflects what the user prescribed, and doesn't depend e.g. on the "steepness" of the filtered field like $\eta$ does. The above transformation prescribed a 1% transformation (i.e. 1 pixel)

[stevengj]

As a test, you could compare it to non-differentiable morphological open/close operations (maybe even setting R_smoothing = 0 to get rid of the subpixel smoothing).

To avoid being too sensitive to discretization effects during testing, you'll probably want to make erosion_dilation ($\Delta d$) larger than a pixel. Note that we still want it to be smaller than the filtering radius, i.e. you want $\Delta x \ll \Delta d \ll \tilde{R}$.

Ideally the difference with morphological operations should go to zero as the resolution increases ($\Delta x \to 0$) for fixed $\Delta d$ and $\tilde{R}$.