Adjoint optimization area constraint

[Paul18811]

Hello everyone!

Currently, I am trying to optimize a structure using the adjoint solver. To make sure that the resulting structure is also manufacturable, I want to implement a constraint on the linewidth and area of the solid and void regions. The constraint on the linewidth is already accomplished as shown in the examples. However, I don't seem to be able to find a function for the area constraint. The constraint that I am talking about is described in the following article:

A. M. Hammond, A. Oskooi, S. G. Johnson, and S. E. Ralph, “Photonic topology optimization with semiconductorfoundry
design-rule constraints,” Opt. Express 29(15), 23916–23938 (2021).

Does anyone know where I can find the function that I am looking for?
Thank you in advance!

[smartalecH]

We haven't added a "plug-and-play" function, like we did with the linewidth/linespacing constraints (yet). But hopefully we can add a tutorial sometime soon.

In the meantime, it's rather straightforward to implement yourself. I just used the findcontours function from skimage. You can do something like this:

from skimage import measure
from scipy import ndimage

''' create a function that 'builds' the appropriate constraint function
and corresponding gradient (differentiable using autograd)'''

def area_constraint(x,min_area):
    # x ........... design variables after mapping
    # min_area .... minimum area constraint (meep units)

    dA = 1/resolution * 1/resolution # differential area
    N = x.size # total number of elements

    # pad to take care of boundary conditions
    pd = int(resolution*min_area*4)
    y_p = np.pad(x,(pd,pd),mode='edge')
    y_p = np.pad(y_p,(10,10),mode='constant',constant_values=0)
    
    # pull the contours
    temp = np.array(y_p)
    contours = measure.find_contours(temp, 0.5)
    
    # build the indicator function
    I = np.zeros(x.shape)
    for c in contours: # loop through contours
        r_mask = np.zeros_like(y_p, dtype='bool')
        r_mask[np.round(c[:, 0]).astype('int'), np.round(c[:, 1]).astype('int')] = 1
        r_mask = ndimage.binary_fill_holes(r_mask)
        area = np.sum(y_p[r_mask].flatten())*dA
        # only worry about a contour if it is smaller than the threshold and sufficiently gray
        if (np.mean(y_p[r_mask].flatten()) > 0.6) and (area < min_area):
            I[mpa.centered(r_mask,x.shape)] = 1
    
    Ja = lambda x_a: npa.dot(x_a.flatten(),I.flatten())/N
    a = Ja(x)
    g = grad(Ja)(x)
    
    return a, g

I haven't ran this... so you might need to tweak some things. But it gives you the basic idea.

(Feel free to submit a PR with a new tutorial if this ends of working for you).

[Paul18811]

Thank you for your reply! Then I will try to implement it myself and will let you know when it works.

[Paul18811]

I implemented the minimum area constraint based on your example. It has some limitations (described below) when there is a large hole inside an island or vice versa. In the end, I didn't use this in my adjoint simulations. It turned out that the robust formulation with a length scale restriction was good enough to prevent small islands/holes in my optimizations.
I don't know whether the code below is good enough for a PR.

from skimage import measure
from scipy import ndimage
from scipy.signal import convolve2d

def area_constraint(x,min_area)
    # x ........... design variables after mapping
    # min_area .... minimum area constraint (meep units)

    dA = 1 / resolution * 1 / resolution  # differential area
    N = x.size  # total number of elements

    y_p = np.pad(x, (10, 10), mode='constant', constant_values=0) #adding extra cells to be able to find contours that not end at the boundary

    # pull the contours
    temp = np.array(y_p)
    contours = measure.find_contours(temp, 0.5)

    #plot the contours
    fig, (ax0,ax1) = plt.subplots(1,2,figsize=(8,8))
    ax0.imshow(temp,plt.cm.gray)
    ax1.imshow(temp,plt.cm.gray)
    for n, contour in enumerate(contours):
        ax1.plot(contour[:,1],contour[:,0], linewidth=2)
    ax1.axis('image')
    ax1.set_xticks([])
    ax1.set_yticks([])
    plt.savefig('contours.png')

    #calculating the amount of solid and void in the contours to distinguish between islands and holes
    area_solid05 = []
    area_void05 = []
    saved_mask05 = []
    I = np.zeros(x.shape)
    for e in contours:  # loop through contours
        r_mask3 = np.zeros_like(y_p, dtype='bool')
        r_mask3[np.round(e[:, 0]).astype('int'), np.round(e[:, 1]).astype('int')] = 1
        r_mask3 = ndimage.binary_fill_holes(r_mask3)
        saved_mask05.append(r_mask3)
        area_solid05.append(np.sum(y_p[r_mask3].flatten()) * dA)
        area_void05.append(np.abs(np.sum(y_p[r_mask3].flatten()) - y_p[r_mask3].size) * dA)

    #adding kernels that will add hole/void space, if the hole/void is too small and needs to be enlarged
    kernel_island = np.asarray([[False, True, False],
                [True, True, True],
                [False, True, False]])
    kernel_hole = np.asarray([[True, False, True],
                [False, False, False],
                [True, False, True]])

    #four categories, so two differentiations. Is it a hole or island? Is it smaller than half the minimum area or is larger than
    #half the minimum area? (and smaller than the general minimum area)
    for n in range(len(area_solid05)):
        if (area_solid05[n] >= area_void05[n]) and (area_solid05[n] <= 0.5 * min_area): #then we have an island that is too small
            y_p[saved_mask05[n]] = 0
        if (area_solid05[n] >= area_void05[n]) and (area_solid05[n] < min_area) and (area_solid05[n] > 0.5 * min_area):#then we have an island that is too small
            y_p[saved_mask05[n]] = 1
            area = y_p[saved_mask05[n]].size * dA
            while area < min_area: #loop that keeps adding material to the island until it is large enough
                result = convolve2d(saved_mask05[n].astype(int), kernel_island.astype(int), mode='same').astype(bool)
                y_p[result] = 1
                area = y_p[result].size * dA
        if (area_solid05[n] < area_void05[n]) and (area_void05[n] <= 0.5 * min_area):#then we have a hole that is too small
            y_p[saved_mask05[n]] = 1
        if (area_solid05[n] < area_void05[n]) and (area_void05[n] < min_area) and (area_void05[n] > 0.5 * min_area): #then we have a hole that is too small
            y_p[saved_mask05[n]] = 0
            area = y_p[saved_mask05[n]].size * dA
            while area < min_area: #loop that keeps adding void space to the hole until it is large enough
                result = convolve2d(saved_mask05[n].astype(int), kernel_hole.astype(int), mode='same').astype(bool)
                y_p[result] = 1
                area = y_p[result].size * dA

    #Note that the above part sees a countour as a hole if the void area is larger than solid area inside the contour.
    #If there would be a large hole inside a contour of a solid, the area of the hole can be larger than that of the solid.
    #The method above would then falsely assume that the solid contour is a hole.
    #It is not expected that this aspect will drastically limit the usage of this function.
    #To be sure, one could inspect the plotted contours to see if such a situation has happened.


    x = y_p[10:-10,10:-10] #removing the extra cells added at the beginning

    #checking the result:
    #opt.update_design([x.flatten()])
    #plt.figure(2)
    #opt.plot2D(True, plot_monitors_flag=False)
    #plt.savefig('patched.png')

    return x

Here you can see two images before and after. They are in binary, but the above code also works in gray scale.