How to create a mask based on points coordinates?

Assuming you are only dealing with rectangles and evenly-spaced points, a purely mathematical approach would be the following:

1. Compute the first grid vector v1 as points[1] - points[0]
2. Iterate over the points, always computing the difference between the current point and the last one. Compare that difference to the grid vector from step 1, if it is different, you hit a new row.
3. Once you have the start point of the second row, compute the second grid vector v2 as the difference between this start point and the very first point of the point list.
4. When you have both grid vectors, you can check if a point is on the edge by computing the four points point + v1, point - v1, point + v2, point - v2. If any of these four points is not in your list, the point is on the edge.

Note that this will only work as long as all numbers can be represented exactly, e.g. when you have integer coordinates.

The following approach can also handle float coordinates:

1. Iterate over the points and compute the difference vector.
2. Compute the dot product of the current difference vector and the previous one.
   • If the dot product is positive, we walked forward within the row
   • If the dot product is negative, the current point is at the start of a new row
3. Keep track of the last row start. All points that follow it are in the last row and therefore part of the border.

The dot product v1 @ v2 is equal to the lengths of the two vectors and the cosine of the angle between them multiplied together. The lengths are never negative, so when the dot product is negative, the angle has to be larger than 90 degrees, meaning the vectors show in different directions.

In Python code, it could look like this (untested):

last_vec = None
first_row = True
last_row_start = None
result = [points[0]]
for p1, p2 in itertools.pairwise(points):
    current_vec = p2 - p1
    if last_vec is None:
        result.append(p2)
        last_vec = current_vec
        continue
    if current_vec @ last_vec < 0:
        if not first_row:
            result.append(p1)
        result.append(p2)
        last_row_start = False
        first_row = False
    if first_row:
        result.append(p2)
# include last row
result.extend(points[point.index(last_row_start)+1:])

This might give surprising results: consider a parallelogram of points like this: [[0, 0], [1, 0], [2, 0], [3, 1], [4, 1], [5, 1], [6, 2], [7, 2], [8, 2]]. The second approach I gave would see this as a 1x9 rectangle.

Maybe another way ... in two major steps

1. find corners indices

To find the corners, I first thought of using argmax but the results differ depending on the orientation of the rectangle. So measuring the center then the distance from the center seemed safer to me.

centroid = coords.mean(axis=0)
dist = np.zeros(coords.shape[0])
for i in range(coords.shape[0]):
    dist[i]= np.linalg.norm(coords[i] - centroid)
corners =np.argwhere(np.round(dist,7) == np.amax(np.round(dist,7)))

in the example corners are therefore :

• left low corner index = 0
• left upper corner index = 42
• right low corner index = 5
• right high corner index = 47

2. compute cross product to find colinear points: Since we've got corners indices, it's now possible to compute cross product and so find all colinear points.

I saw two options here :

a. compute all cross products = 0

b. stop computation when the first null cross product is found because it's correspond to the delta we are searching

for i in range(left_low_corner+1, coords.shape[0]):
    vec_1 = coords[left_high_corner]-coords[left_low_corner]
    vec_2 = coords[i]-coords[left_low_corner]
    cross_p = np.cross(vec_1,vec_2)
    if abs(cross_p[2]) < 0.001:
        VD = i
        break

When VD is known, it's now easy to retrieve left and right sides indices

side_le_index = np.arange(left_low_corner,left_high_corner+delta_verti,delta_verti)
side_ri_index = np.arange(right_low_corner,right_high_corner+delta_verti,delta_verti)

This also now possible to get the HD step value and therefore indices of vertices of lower and upper sides:

nb_column = int(coords.shape[0]/len(side_le_index))
step = int((right_low_corner+1)/nb_column)
side_lo_index = np.arange(left_low_corner,right_low_corner+step, step) 
side_up_index = np.arange(left_high_corner,right_high_corner+step, step) 

it may seem a little crooked but it works whatever the orientation of the rectangle and whatever the horizontal and vertical deltas.

If anyone has any tips for improving the code or a more efficient way to approach the problem please let me know!