implement faster sphere-based subpixel averaging

[stevengj]

Right now, to support subpixel averaging, we have a function that computes the overlap fraction of the intersection of an object with a box; it's quite complicated and often involves a numerical quadrature step. Recently, we've realized that there is a much simpler way to accomplish an equivalent task. The two key ideas are:

• We only care about subpixel smoothing asymptotically at high resolutions, where the pixels are small compared to the object, and only to 1st-order accuracy in the resolution. In this limit, almost everywhere we can approximate the object within a pixel as a plane (just from the normal vector and the closest point of the object to the center). This makes it much easier to work with.
• Once you have a plane, you could compute its intersection with a box analytically — @wsshin actually implemented this in Julia, but it is complicated because there are so many cases (depending on which edges of the box intersect the plane). However, since subpixel smoothing corresponds to convolving the structure with a smoothing window, there is no particular need to convolve with a cube smoothing window that matches the pixels/voxels. Instead, you can can convolve with a spherical smoothing window, which makes the problem trivial: it is easy to compute the intersection of a sphere with a plane via the volume of a spherical cap. (This is, in fact, essentially what we are doing with the new smoothed-projection scheme for TopOpt in Meep: First-order subpixel smoothing for density-based TO meep#2741.)

Not only is this more elegant (making it easier to add new kinds of shape) and faster, but it is also more accurate because it avoids the numerical quadrature step that we currently have. This is important if we want to use finite differences for shape optimization: for example, suppose you want the derivative of a solution with respect to a perturbation of a sphere radius R, a basic step in an adjoint method would be to compute $\partial(\varepsilon^{-1})/\partial R$, which we might approximate by a finite difference of the discretized $\varepsilon^{-1})$ grid for a small change $R + \delta R$, but doing this accurately requires the subpixel smoothing to be computed accurately (ideally analytically) so that it is sensitive to very small $\delta R$ up to the limits of floating-point precision.