Get list of coordinates contained in Polygon

[GeneralErrorOK]

I'm a very happy user of Shapely within our video analysis tool. This tool offers the user a possibility to analyze certain parts of a video stream. Our users can draw polygonal shapes with the video screen to define their regions of interest. We use Shapely to store these regions and to perform some basic math on them. Works like a charm!

But, we've stumbled upon an interesting issue. We use Nvidia's Optical Flow technology to look for motion. For every video frame, we get a 3-dimensional array with X and Y movement vertices. These are simpel floats that define how much this block of 4x4 pixels has moved compared to the previous videoframe. Very elegant, very simple.

But, we are only interested in the blocks (of 4x4 pixels) that fall within the RIO('s) the user has specified. So for every video frame we need to select these exact elements from the 3D array before we do any more calculations on them (thresholds, etc.).

For a 1920x1080 video the array has a shape of 480:270:2. We need the indexes of the first two dimensions that (times 4) match the coordinates of the Polygon the user defined (because the coordinates of the Polygon match video resolution, not OF array resolution).

I've solved this simply by (once) looping through all OF vertices, and constructing a new Polygon out of the indexes and then using the contains() method to verify whether or not this exact block of optical flow vertices falls within the user's ROI. Like so:

vector_coords = []
video_width, video_height, _, _ = self.configuration.video_specs
vector_array_width = video_width // crop_factor
vector_array_height = video_height // crop_factor
for vx in range(vector_array_width):
    for vy in range(vector_array_height):
        px = vx * crop_factor
        py = vy * crop_factor
        shell = (
            (px, py),
            (px+crop_factor, py),
            (px+crop_factor, py+crop_factor),
            (px, py+crop_factor),
            (px, py)
        )
        vector_bounds = Polygon(shell)
        if detection_area.polygon.contains(vector_bounds):
            vector_coords.append((vx, vy))

where vector_coords now contains a list of x/y coordinates we can later use to "cut" all movement information from every 3d array. We can cut these in mere milliseconds by simply intersecting the lists so that works fine.

But... the algorithm we use above is VERY slow. Understandably, because we create a Polygon, check if this new polygon is contained within an existing one and we do this 480 x 270 = 129600 times! On our machine it takes ~ 22 seconds which is unacceptable.

Question is: is there a faster way to get a list of all (integer) coordinate pairs of all points within a Polygon? Without having to create a new "candidate" Polygon and checking whether it fits inside? It works but its very, very inefficient.

[brendan-ward]

The contains predicate can be kind of slow depending on the inputs.

A few suggestions, not sure if they'll help in this case (sorry, no time to write code example):

• try to use vectorized operations where you can instead of for loops. E.g., create all your vector_bounds polygons at once then query those
• use a spatial index of the array of vector_bounds with an appropriate predicate (intersects and contains_properly should both be pretty fast); if you really need to use the contains predicate you can use it against the results that come from the tree first to cut down on the number of things you are testing with it
• use box to create your vector_bounds polygons

If you can transform your inputs that you use to create individual vector_bounds polygons into arrays of xmin, ymin, etc, then calling box(xmins, ymins, ...) is a vectorized method to create all those vector_bounds polygons in one vectorized method.

Also see contains_xy which is a point-in-polygon method that may or may not be faster if you query it 4 times (e.g., once with array of xmins, once with array of ymins, etc). I haven't done a test comparing the performance of this vs querying boxes. But even in this case, using the spatial index first will make things faster.

Note: contains_properly is able to use prepared geometries under the hood (both in STRtree query or used directly), so it is considerably faster than contains, but has slight differences regarding the relationship to the boundary of the polygon.

[GeneralErrorOK]

Thank you for the excellent guidance, brendan-ward! I did not fully understand everything you said - the spatial index STRTree part is hard to follow for me as non-math major ;), but I did manage to shave off 90% of the calculation time by avoiding the nested for loops and using list comprehensions and vectorized operations to calculate everything.

What I do now is:

• create lists of all xmin, xmax, ymin, ymax coordinates for the vector boxes (default is 4x4 pixels)
• create Polygons with the box() predicate (this is FAST!)
• check containment of all vector boxes in the detection area
• filter out all boxes that are not contained with numpy.extract()
• retrieve indexes (which is pixel coordinate // crop_factor) of the vector boxes from the created/filtered polygons

We've gone from ~22 seconds to ~2 seconds which is WAY more acceptable.

video_width, video_height, _, _ = self.configuration.video_specs

xmins = [px for px in range(0, video_width, crop_factor)] * (video_height // crop_factor)
ymins = [py for py in range(0, video_height, crop_factor)] * (video_width // crop_factor)
xmaxs = [px for px in range(0+crop_factor, video_width+crop_factor, crop_factor)] * (video_height // crop_factor)
ymaxs = [py for py in range(0+crop_factor, video_height+crop_factor, crop_factor)] * (video_width // crop_factor)

vector_boxes = box(xmins, ymins, xmaxs, ymaxs)
prepare(detection_area.polygon)
results = detection_area.polygon.contains_properly(vector_boxes)

vector_boxes_in_detection_area = np.extract(results, vector_boxes)

vector_coords = [(int(vector_box.bounds[0] // crop_factor), int(vector_box.bounds[1] // crop_factor)) for vector_box in vector_boxes_in_detection_area]

But I'm still wondering if I can optimize this even more. It seems like I'm doing too many operations, still.

[brendan-ward]

You should be able to use the spatial index query method, after constructing your boxes (it will also prepare your input geoms automatically):

from shapely import STRtree

tree = STRtree(vector_boxes)
indices = tree.query(detection_area.polygon, predicate="contains_properly")
vector_boxes_in_detection_area = vector_boxes[indices]

bounds = shapely.bounds(vector_boxes_in_detection_area)  # returns array: [[xmin1, ymin1, xmax1, xmax1], [xmin2, ymin2, xmax2, ymax2]]

vector_coors = bounds[:,:2] // crop_factor  # this should be equivalent to what you were doing before

The spatial index should make this faster because it does a fast check for overlap between the bounding box of the detection_area.polygon and each of your boxes, and then only submits those with overlapping bounding boxes to the contains_properly predicate.

[GeneralErrorOK]

Wow! Thanks! We went from 22 seconds to 2 seconds to ~300ms with your solution. That's incredible. Thanks a lot for your help!