Signed-distance functions for MaterialGrid Levelsets

[smartalecH]

We've recently added subpixel smoothing to MaterialGrids. One nice side effect of this new feature is the ability to optimize when β=∞ (i.e. shape optimization).

Often, shape optimization is used when we already have a good initial guess and just want to improve the performance of that device by perturbing the boundaries. In order to supply an initial guess as a MaterialGrid, however, the underlying design variables must be slowly varying, which is never the case when we are starting with a binary structure...

Luckily, we can create a signed-distance function (e.g. using the fast-marching method) from our binary initial guess. The scikit-fmm python package does this rather well. For example, suppose we would like to supply a simple circle as our initial guess:

design_shape = mp.Vector3(1,1,0)
design_region_resolution = 1000
Nx = int(design_region_resolution*design_shape.x)
Ny = int(design_region_resolution*design_shape.y)
x = np.linspace(-0.5*cell_size.x,0.5*cell_size.x,Nx)
y = np.linspace(-0.5*cell_size.y,0.5*cell_size.y,Ny)
xv, yv = np.meshgrid(x,y)
design_params = np.sqrt(np.square(xv) + np.square(yv)) < rad

If we want an equivalent representation which is amenable to subpixel smoothing, we can do the following:

def make_sdf(data):
    '''
    Assume the input, data, is the desired output shape
    (i.e. 1d, 2d, or 3d) and that it's values are between
    0 and 1.
    '''
    # create signed distance function
    sd = skfmm.distance(data- 0.5 , dx = 1)

    # interpolate zero-levelset onto 0.5-levelset
    x = [np.min(sd.flatten()), 0, np.max(sd.flatten())]
    y = [0, 0.5001, 1]
    f = interpolate.interp1d(x, y, kind='linear')

    return f(sd)

design_params_sdf = make_sdf(design_params)

which results in a nice, smoothly varying function:

Now, if we project the above signed-distance function..

plt.imshow(mpa.tanh_projection(design_params,np.inf,0.5))

we can even verify that the difference between the projected levelset and the initial starting condition is indeed 0:

plt.imshow(mpa.tanh_projection(design_params,np.inf,0.5)-design_params)

cc @oskooi this might make a nice tutorial someday (assuming we combine it with some actual shape optimization).

[yugeniom]

Hi there! I am trying to implement something like idea 3 of the proposed tutorial 2011.
As far as I understand, including the signed-distance-function remapping not just during the grid initialization, but as part of the objective function, requires it to be autograd-differentiable for the backpropagation to work. I assume skfmm package isn't part of autograd, so it won't be differentiable out of the box; do you know some differentiable alternative of the same fast-marching algorithm? Does maybe autograd include it somehow?
thanks!

PS
Is the need for such a signed-distance function related to the ability to reconstruct the orginal binary profile with the maximum accuracy maybe? Do you think some other autograd-provided smooth remapping - say the mpa conic filter you use in some meep adjoint examples - could do the job somehow as well?

[smartalecH]

@oskooi had the same question as you, so it's worth going through things carefully.

There are many reasons to use a level set, but for the sake of this discussion, let's focus on just two (which are somewhat disjoint):

Case 1: (this one is not as relevant to your problem)

Imagine you've got some shape for your design. Perhaps this design came from some physical intuition you have, or even from another optimization approach. But now you just want to do some "touch-up" optimization on this design -- maybe by perturbing the boundaries until your FOM is achieved (this is sometimes referred to as shape optimization...). In this case, you don't care about the underlying parameterization. In fact, the more general the parameterization, the better, because you don't want to limit the expressivity of the final design (other than by using fabrication constraints, perhaps).

To do this within meep's framework, we need to map our original design to a format that meep can manipulate (e.g. using material grids). A common way to do this is using signed distance functions (SDF), which create a level-set approximation to some binary representation. SDFs have lots of cool properties, like unit slope, etc. But the important thing is we can use them to take any arbitrary starting guess, and do even further arbitrary shape optimization on it by perturbing its boundaries.

Again, we don't care about the parameterization in this case. We just want to start with a particular design and ensure flexibility throughout the optimization.

Because of this, we actually don't need the SDF generating function to be differentiable. We only use it to generate our initial density variables, and we never use it again!

Case 2: This one is more relevant to your example.

Suppose you do care about the underlying parameterization. This happens when the underlying structure has distinct properties that we want to manipulate. For example, a grating coupler can be defined by duty cycle, etch depth, chirping profile, etc (although free-form grating couplers are really cool).

In this case, we need to come up with a level set that is meaningfully parameterized (we can't just generate an SDF in this case because that isn't tied to any parameters!).

Let's imagine we have a 1D grating coupler we want to optimize. For simplicity sake, we can assume we are only optimizing the grating period (Λ). We can describe our grating structure using the following function,

$$ f(\Lambda) = \frac{1}{2}\mathrm{sin}(\frac{2\pi}{\Lambda}) + \frac{1}{2} $$

where $f$ varies from 0 to 1 (the range of values accepted by our material grid). Luckily, this is a very smooth function, so we can directly pass the output of this function to our MaterialGrid. If we set $\beta=\infty$, then the subpixel smoothing routine will project $f$ for values $\ge 0.5$ (resulting in a binary grating structure).

Now, since the material grid densities are smoothly varying, meep can compute the gradient $\frac{d \mathrm{FOM}}{d f}$, and then we can manually backpropagate the rest of the gradient using autograd's chain rule (like we usually do with mapping functions).

Of course, you typically want to do much more sophisticated parameter optimizations in 2D and 3D. So the challenge becomes crafting a levelset function that is smoothly varying on the output, but meaninfully depends on the parameters you care about.

For example, to extend the 1D grating example to also include duty cycle for multiple layers, you could take a look at Andy Michael's paper.