polyfillaa

sofia_redux.toolkit.image.fill.polyfillaa(px, py, xrange=None, yrange=None, start_indices=None, area=False)[source]

Finds all pixels at least partially inside a specified polygon

Find the number of vertices for each polygon and then loop through groups of polygons with an equal sides. Then for each group of similar sided polygons:

  1. Create a shared pixel grid to be used by all polygons in the group. The size of the grid is determined by the maximum range of polygon widths and heights in the group.

  2. For each polygon edge, calculate where it crosses the left and bottom edges of the each grid cell.

  3. Determine whether the vertices are inside the grid points. The crossings determined in step 2 can be used to determine whether the lower-left grid points are contained within a polygon. In this case, a grid point is said to be in the polygon if there are an odd number of intersection points with the polygon along the y-coordinate of the grid point.

  4. Remember that we have only calculated polygon crossings on the left and lower edges of each cell and whether the lower-left corner of a cell is enclosed in the polygon. To avoid duplicate calculations, take note of the following facts:

    1. If a lower-left grid point (gp) is enclosed in the polygon then the lower-right gp of the cell to the left, the upper-right gp of the cell to the lower-left, and the upper-left gp of the cell below must also be enclosed.

    2. If the polygon crosses the left edge of the cell, it must also cross the right edge of the cell to the left.

    3. If the polygon crosses the bottom edge of the cell, it must also cross the top edge of the cell below it.

  5. Given all of these points, it is clear that the maximum number of clipped polygon vertices will be equal to 2 times the number of polygon vertices of the input group of polygons since each edge can cross a maximum of 2 cell sides, or be clipped to where one vertex remains inside the cell and the other is located on the edge of the cell. If both vertices are inside the cell then that edge remains unchanged. Therefore the maximum number of clipped polygon vertices occurs when all polygon edges intersect 2 cell edges (imagine two superimposed squares and then rotating one by 45 degrees). We create a 3D output array containing the new vertices for each polygon and for each pixel and fill it in the following order:

    1. vertices that are inside

    2. grid points that are inside

    3. clipped vertices

  6. If we do not need to calculate area, then we can stop here since we have all the new vertices. However, if we do need to calculate area, then these points need to be arranged in clockwise or anti-clockwise order. This is done by ordering the points based on angle with respect to the center-of-mass of the points. As a side note, this takes up a significant amount of processing time (the sorting, not the angle). I attempted many alternate solutions to this by not changing the original order of the input vertices and clipping in-place. One of the most promising methods was to encode where a vertex was in relation to the cell and then clip based on that, ordering by keeping points in the following manner and looping around the edges of a cell in the order bottom -> left -> top -> right:

    outside -> inside = keep the inside and clipped vertices inside -> outside = clipped vertex only inside -> inside = keep the second inside point outside -> outside = keep both clipped vertices

    The problem occurs with whether grid points (corners) should be included or not, and where they are located in the final order of points. This can be achieved by encoding a vertex location relative to the cell as bits where 1 indicates outside and 0 indicates inside in the order bottom-left-top-right. So for example, 0000 indicates a point is inside the cell and 0100 indicates a vertex is to the left of the cell. codes containing two 1’s indicate corners, and a set of rules can then be established determining whether a gp is inside or outside based on the code combination of vertices for one edge.

    However, in order to achieve vectorization, this requires storing at least 16 times the amount of initial data in the cache (number of polygons * area of polygon * 16) which can be huge and clumsy. We can get around this with loops, but this is not efficient with Python. Therefore, it is quicker and safer to just use a sort on the angle and be done with it. If you think there’s a better solution then please feel free to implement it (and tell me too for my own curiosity).

  7. Calculate area using the shoelace formula:

    A = 0.5 * sum( x_i * y_(j+1)) - sum(x_(i+1) * y(j))|

  8. Finally, organize the results by lumping together polygons based on the number of cells enclosed within. This allows us to grab which cells belong to which polygons quickly.

Parameters:
pxarray_like of (int or float)

Contains the x-coordinates of the polygons. May be provided as a flat array if start_indices is provided. Alternitvely, a 2-level nested list may be provided with each sublist containing vertices for a single polygon.

pyarray_like of (int or float)

Contains the y-coordinates of the polygons. May be provided as a flat array if start_indices is provided. Alternitvely, a 2-level nested list may be provided with each sublist containing vertices for a single polygon.

xrangearray_like, optional

(2,) [lower x-limit, upper x-limit] array defining the x range on which the polygon is superimposed (number of columns). size of the pixel grid on which the polygon is superimposed Supplying it will clip all x results to within xrange[0]->xrange[1].

yrangearray_like, optional

(2,) [lower y-limit, upper y-limit] array defining the y range on which the polygon is superimposed (number of rows). size of the pixel grid on which the polygon is superimposed Supplying it will clip all y results to within yrange[0]->yrange[1].

start_indicesarray_like of int, optional

Multiple polygon shapes may be specified with the polygons parameter by specifying indices in the first dimension of polygons which should be considered a first vertex of a polygon. For example, the nth polygon consists of vertices at polygons[n:(n+1)]. Note that start_indices will be sorted and that the last index (start_indices[-1]) gives the last point belonging to start_indices[-2]. The last polygons vertex is not automatically appended. i.e. the last polygon is polygons[start_indices[-2]:start_indices[-1]], not polygons[start_indices[-1]:].

areabool, optional

if True, return an additional dictionary containing the area

Returns:
polygon, [areas](2-tuple of (dict or array)) or (dict or array)
A dictionary containing the output:

polygon index (int) -> tuples of (y, x) grid coordinates contained within each polygon

areas (optional if area is set to True)
polygon index (int) -> list of pixel areas in the same order as

given by output indices.

Notes

I vectorized the crap out of this thing and removed redundancies so it could run faster than the IDL version. The IDL version was loop driven which is generally a no no, but it this case it was very very well done and also able to use C compiled clipping code. Python looping using the same method could not even slightly keep up. For example, 50,000 polygons on a 256, 256 grid, where each polygon covered an area of 3x3 grid pixels took 20 seconds using the Python equivalent of the IDL code with all the speed saving tricks / data types and objects available. In comparison, IDL took ~3 seconds while this method took ~1 second. If you want to attempt speeding it up further then the main choking point is the sorting of angles when calculating area. If you want to use the old method for calculating output polygons/area, replace the main body of code after normalization to a regular pixel grid to loop through all polygons, then all pixels covered by the polygon and calculate the area using sofia_redux.toolkit.polygon.polyclip.