Subpixel Averaging

[bruxillensis]

Is anyone working on this? Would it be alright if I did? Any effect on other work/issues ongoing right now?

[l-berg]

Hi. In general, any work on subpixel averaging would be highly appreciated.

I am working on some minor adjacent details: Currently, FDTDX uses linear interpolation of the inverse permittivities to interpolate between materials (for continuous devices and mulit-material objects that use a voxel_mask). I am changing this to use permittivities instead of inverse permittivities, and also adding etching support for devices (i.e., support devices that do not fully define their grid_shape, but leave some of the existing materials intact).

These changes should hopefully come later this week, are only minor in scope and do not affect any interfaces.

[l-berg]

Not sure whether this is relevant, but I also created a smooth version of the Polygon. It uses shapely to compute the fraction of how much each pixel is covered by the polygon. Since shapely scales poorly, I use the matplotlib Path test for most pixels, and shapely only for border pixels.

import shapely
from matplotlib.path import Path
from scipy.ndimage import binary_dilation, binary_erosion
def polygon_to_mask_hybrid(
    boundary: tuple[float, float, float, float],
    resolution: float,
    polygon_vertices: np.ndarray,
    safe_margin: int = 2
) -> np.ndarray:
    """
    Generate a 2D fractional mask, highly optimized for large grids.
    Mixes fast binary point-in-polygon with exact Shapely boundary intersections.
    """
    assert polygon_vertices.ndim == 2
    assert polygon_vertices.shape[1] == 2
    min_x, min_y, max_x, max_y = boundary

    # 1. Create coordinate arrays
    x_coords = np.arange(min_x, max_x + 0.5 * resolution, resolution)
    y_coords = np.arange(min_y, max_y + 0.5 * resolution, resolution)
    X, Y = np.meshgrid(x_coords, y_coords, indexing="ij")

    # 2. Fast Binary Masking (Core Interior/Exterior)
    points = np.column_stack((X.ravel(), Y.ravel()))
    polygon_path = Path(polygon_vertices)
    inside_polygon = polygon_path.contains_points(points, radius=1e-3*resolution)
    binary_mask = inside_polygon.reshape(X.shape)

    # 3. Identify the Boundary Band using Scipy Morphology
    # Dilation pushes the boundary out, erosion pulls it in. 
    # The difference between them is our "uncertainty band" where edges lie.
    dilated = binary_dilation(binary_mask, iterations=safe_margin)
    eroded = binary_erosion(binary_mask, iterations=safe_margin)
    boundary_mask = dilated ^ eroded  # XOR gives us only the boundary shell

    # Initialize the output float array with our fast binary results
    fractional_mask = binary_mask.astype(float)

    # Get the 2D array indices of the boundary pixels
    rows, cols = np.where(boundary_mask)

    # 4. Run exact Shapely calculation ONLY on the boundary pixels
    if len(rows) > 0:
        # Extract physical coordinates for just the boundary pixels
        cx = X[rows, cols]
        cy = Y[rows, cols]
        
        half_res = resolution / 2.0
        
        # Create vectorized Shapely boxes for the boundary pixels
        boxes = shapely.box(
            cx - half_res, cy - half_res, 
            cx + half_res, cy + half_res
        )
        
        poly = shapely.Polygon(polygon_vertices)
        
        # Calculate exact intersections
        intersections = shapely.intersection(boxes, poly)
        areas = shapely.area(intersections)
        
        # Overwrite the boundary band with exact fractions (0.0 to 1.0)
        pixel_area = resolution ** 2
        fractional_mask[rows, cols] = areas / pixel_area

    return fractional_mask