Python: Boundary points of a non-convex grid

Question

I have the following situation as shown in the Figure below:

I want to find out the grid points that surround the red points. The red points are trajectories of moving agents. So in many situations we have a bunch of points, therefore the solution should be as fast as possible.

The grid is plotted as points.

First step, I managed to reduce the number of grid points as shown below (plotted as x):

This is my code:

step = .5
gridX, gridY = np.meshgrid(np.arange(xmin-step, xmax+step, step), np.arange(ymin-step, ymax+step, step))
mask = False * np.empty_like(gridX, dtype=bool)
threshold = 0.5
for (x,y) in  zip(df_traj['X'], df_traj['Y']):
    pX = x * np.ones_like(gridX)
    pY = y * np.ones_like(gridY)
    distX = (pX - gridX)**2
    distY = (pY - gridY)**2
    dist = np.sqrt(distX + distY)
    condition = (dist < threshold)
    mask = mask | condition

gX = gridX*mask
gY = gridY*mask

Second step, and this is where I need a little help:

How can I efficiently filter out the inner points of the grid and keep only the "x-points" outside the "red area"?

EDIT

In this special case I have 92450 red points.

So, your input is given in the variable `df_traj`? And approximately how many points do you have? "A lot"? — user202729, Mar 05 '18 at 13:12
If I understand your question correctly I believe you are looking to find the [non-convex hull](https://stackoverflow.com/questions/9830218/determine-non-convex-hull-of-collection-of-line-segments) of your points — Cory Kramer, Mar 05 '18 at 13:13
@gbtimmon The OP said explicitly that the boundary is not convex. — user202729, Mar 05 '18 at 13:15
If I understood the problem correctly removing all points that have exactly four neigbors should do the trick, or did I miss anything? — MB-F, Mar 05 '18 at 13:20
The point (0.5, 0) has 4 neighbors, but some of them are "within" the red area. But this point should bot be removed.. — Tengis, Mar 05 '18 at 13:24
Really sounds like you are looking for a concave hull algorithm. I have implemented one for Python a while back: [hull](https://github.com/jsmolka/hull) (based on a paper in the repo). — jsmolka, Mar 05 '18 at 13:25
@jsmolka looks good. I was really hoping to find some simple "mask" as was the case with the "x" points. Thanks I'll have a look at your solution — Tengis, Mar 05 '18 at 13:28
@Tengis I agree with your intuition that there should be a simple solution considering that you already have reduced it to a discrete grid. It all comes down to distinguishing between points like (0.5, 0) and say (-2, -0.5). Seems to me it has to do with the curvature of the enclosing shape, but I'm not quite there yet :) — MB-F, Mar 05 '18 at 13:36

score 2 · Answer 1 · answered Mar 05 '18 at 14:02

I think if you just walk around the edge, since its a evenly spaced grid, it should work. No need for far more complicated non-convex-hull to handle pnts that can be anywhere. This isn't adapted to your code and I cheat with my data structures to make the code easy so youll have to handle that but it think as psuedocode it should work.

pnts = <<lists of points>>
edge_pnts = []
fpnt = pnt_with_min_x_then_min_y
cpnt = fpnt
npnt = None 

while npnt != fpnt:

    if   (cpnt[0] + 1, cpnt[1]    ) in pnts: npnt = (cpnt[0] + 1, cpnt[1]    )
    elif (cpnt[0] + 1, cpnt[1] + 1) in pnts: npnt = (cpnt[0] + 1, cpnt[1] + 1)
    elif (cpnt[0],     cpnt[1] + 1) in pnts: npnt = (cpnt[0]    , cpnt[1] + 1)
    elif (cpnt[0] - 1, cpnt[1] + 1) in pnts: npnt = (cpnt[0] - 1, cpnt[1] + 1)
    elif (cpnt[0] - 1, cpnt[1]    ) in pnts: npnt = (cpnt[0] - 1, cpnt[1]    )
    elif (cpnt[0] - 1, cpnt[1] - 1) in pnts: npnt = (cpnt[0] - 1, cpnt[1] - 1)
    elif (cpnt[0]    , cpnt[1] - 1) in pnts: npnt = (cpnt[0]    , cpnt[1] - 1)
    elif (cpnt[0] + 1, cpnt[1] - 1) in pnts: npnt = (cpnt[0] + 1, cpnt[1] - 1)
    else: raise ValueError("Oh no!")

    edge_pnts.append(npnt)
    cpnt = npnt

I think you're assuming convexity here. you should take into account the direction from which you came, and always try going e.g. CCW from it. For example, your algorithm will fail in the point (0.5,-0.5) in the OP's example — shayelk, Mar 05 '18 at 14:30

shayelk · Answer 2 · 2018-03-06T06:51:50.777

You just need to pick a point you know is on the hull (let's take the leftmost point among the topmost points), and assume you "got to it" from above (as we know there is no points above it). now while the next point is not in your list:

Try going CCW from the direction you came from.

The code looks like that:

matrix = [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
         [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
         [0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]

# Find the leftmost topmost point
first_point = None
for i in range(len(matrix)):
    if first_point:
        break
    for j in range(len(matrix[0])):
        if matrix[i][j]:
            first_point = [i, j]
            break

next_point = first_point
prev_direction = 'up'
next_direction_dict = {'up': 'left', 'left': 'down', 'down': 'right', 'right': 'up'}
opposite_direction = {'up': 'down', 'left': 'right', 'down': 'up', 'right': 'left'}
hull_points = []


def go_direction(point, direction):
    # Find the point to a given direction of a given point
    i = point[0]
    j = point[1]
    if direction == 'right':
        j += 1
    elif direction == 'up':
        i -= 1
    elif direction == 'left':
        j -= 1
    elif direction == 'down':
        i += 1
    else:
        raise ValueError
    return [i, j]


def find_next_point(matrix, point, prev_direction):
    next_direction = next_direction_dict[prev_direction]
    next_point = go_direction(point, next_direction)
    prev_direction = next_direction
    while not matrix[next_point[0]][next_point[1]]:
        next_direction = next_direction_dict[prev_direction]
        next_point = go_direction(point, next_direction)
        prev_direction = next_direction
    from_direction = opposite_direction[prev_direction]
    return next_point, from_direction

next_point, prev_direction = find_next_point(matrix, next_point, prev_direction)
while next_point != first_point:
    if next_point not in hull_points:
        hull_points.append(next_point)
    next_point, prev_direction = find_next_point(matrix, next_point, prev_direction)

Edit: Now also handles single point 'tentacles' by iterating until returning to the first point

armatita · Answer 3 · 2018-03-05T16:45:12.450

1

For non convex polygons, like your example, convex hull is not a solution. My recommendation is that, given you already have a discrete grid, that you simply attribute the value False to a bool grid cell when a sample occurs inside. Something like this:

import numpy as np
import matplotlib.pyplot as plt

# Generic data production
X, Y = np.random.normal(0, 1, 100000), np.random.normal(0, 1, 100000)
ind = np.where((X > 0) & (Y > 0))
X[ind] = 0
Y[ind] = 0    

# Generic grid definition
step  = 0.5
xmin, xmax = X.min(), X.max()
ymin, ymax = Y.min(), Y.max()
firstx = xmin-step/2
firsty = ymin-step/2
lastx  = xmax+2*step/2
lasty  = ymax+2*step/2
gridX, gridY = np.meshgrid(np.arange(firstx, lastx, step), np.arange(firsty, lasty, step))

# This is the actual code that computes inside or outside
bool_grid = np.ones(gridX.shape, dtype="bool")
bool_grid[np.int_(0.5+(Y-firsty)/step), np.int_(0.5+(X-firstx)/step)] = False

# Plot code
plt.scatter(gridX.flatten(), gridY.flatten(), marker="+", color="black", alpha=0.3)
plt.scatter(gridX[bool_grid].flatten(), gridY[bool_grid].flatten(), marker="+", s=90, color="green")
plt.scatter(X, Y, s=10, color="red")
plt.show()

, which results in the following (green crosses are True values):

NOTE: This is very fast but it has some limitations. If your data is not compact you'll have True values inside the shape (so holes are possible). It's possible to process the image to remove the holes however (a flood fill or a moving window based algorithm for example). Another possibility is to play with the resolution of the grid.

edited Mar 05 '18 at 16:45

answered Mar 05 '18 at 15:09

armatita

12,825
8
48
49

What if instead of analyzing the shape -- you analyzed the space surrounding the space and found all connected points that are outside. Then the holes problem goes away and you can find the edges quite quickly. You might need to add some padding to ensure there is a continuous outside, but thats trivial. – gbtimmon Mar 05 '18 at 16:33
Wait i dont understand how you are finding the edge at all here ? – gbtimmon Mar 05 '18 at 16:34
@gbtimmon That's because you are thinking in shapes (polygons) while the algorithm is working for points. Think of the cells of the grid like boxes. For each sample calculate its position on the grid and turn that "box" to false. It's an operation compatible with broadcasting. That is the reason why it's so fast. The shape will be a natural consequence. As for the holes, my note was simply a warning about something that might not be completely evident from the code. But as far as I know the "possible" holes can actually be something positive for the OP. – armatita Mar 05 '18 at 16:42
@gbtimmon From what I understood the OP is looking for the equivalent of the green cross markers in my example plot (so `gridX[bool_grid]` and `gridY[bool_grid]`). At least that's what I interpret from his question: *"How can I efficiently filter out the inner points of the grid and keep only the "x-points" outside the "red area"?"* – armatita Mar 05 '18 at 16:52
Thats to simple, he already has the inner points. So then its just all points - inner points. He wanted to filter inner points of his polygon to get a set of border points I thought. – gbtimmon Mar 05 '18 at 17:57

score 0 · Answer 4 · answered Mar 05 '18 at 16:40

Heres another though that came to my mind --

A flood fill of the space:

pnts = <<lists of points>>
seen = set()
edges = []
stack = (0,0)

while stack:
    ele = stack.pop()
    if ele in pnts:
        edges.append(ele)
    else:
        seen.add(ele)
        if (ele[0] + 1, ele[1]) not in seen: 
            stack.append(ele[0] + 1, ele[1])
        if (ele[0] - 1, ele[1]) not in seen: 
            stack.append(ele[0] - 1, ele[1])
        if (ele[0], ele[1] + 1) not in seen: 
            stack.append(ele[0], ele[1] + 1)
        if (ele[0], ele[1] - 1) not in seen: 
            stack.append(ele[0], ele[1] - 1)

Then you need to sort the points which shouldn't be too hard.

Python: Boundary points of a non-convex grid

EDIT

4 Answers4