Eugene W.
Eugene W.

Reputation: 446

Split the upper edge line and lower edge line for a given cloud of points?

I have a list of 2D points x,y. And I need to find a smooth curve for the upper and lower edges (red and blue curves, correspondingly). See the picture below: enter image description here

Here I've found a good example, where the outer edge of x,y points is detected. Using these I have work I have:

# https://stackoverflow.com/a/50714300/7200745
from scipy.spatial import Delaunay
import numpy as np

def alpha_shape(points, alpha, only_outer=True):
    """
    Compute the alpha shape (concave hull) of a set of points.
    :param points: np.array of shape (n,2) points.
    :param alpha: alpha value.
    :param only_outer: boolean value to specify if we keep only the outer border
    or also inner edges.
    :return: set of (i,j) pairs representing edges of the alpha-shape. (i,j) are
    the indices in the points array.
    """
    assert points.shape[0] > 3, "Need at least four points"

    def add_edge(edges, i, j):
        """
        Add an edge between the i-th and j-th points,
        if not in the list already
        """
        if (i, j) in edges or (j, i) in edges:
            # already added
            assert (j, i) in edges, "Can't go twice over same directed edge right?"
            if only_outer:
                # if both neighboring triangles are in shape, it's not a boundary edge
                edges.remove((j, i))
            return
        edges.add((i, j))

    tri = Delaunay(points)
    edges = set()
    # Loop over triangles:
    # ia, ib, ic = indices of corner points of the triangle
    for ia, ib, ic in tri.vertices:
        pa = points[ia]
        pb = points[ib]
        pc = points[ic]
        # Computing radius of triangle circumcircle
        # www.mathalino.com/reviewer/derivation-of-formulas/derivation-of-formula-for-radius-of-circumcircle
        a = np.sqrt((pa[0] - pb[0]) ** 2 + (pa[1] - pb[1]) ** 2)
        b = np.sqrt((pb[0] - pc[0]) ** 2 + (pb[1] - pc[1]) ** 2)
        c = np.sqrt((pc[0] - pa[0]) ** 2 + (pc[1] - pa[1]) ** 2)
        s = (a + b + c) / 2.0
        area = np.sqrt(s * (s - a) * (s - b) * (s - c))
        circum_r = a * b * c / (4.0 * area)
        if circum_r < alpha:
            add_edge(edges, ia, ib)
            add_edge(edges, ib, ic)
            add_edge(edges, ic, ia)
    return edges

def find_edges_with(i, edge_set):
    i_first = [j for (x,j) in edge_set if x==i]
    i_second = [j for (j,x) in edge_set if x==i]
    return i_first,i_second

def stitch_boundaries(edges):
    edge_set = edges.copy()
    boundary_lst = []
    while len(edge_set) > 0:
        boundary = []
        edge0 = edge_set.pop()
        boundary.append(edge0)
        last_edge = edge0
        while len(edge_set) > 0:
            i,j = last_edge
            j_first, j_second = find_edges_with(j, edge_set)
            if j_first:
                edge_set.remove((j, j_first[0]))
                edge_with_j = (j, j_first[0])
                boundary.append(edge_with_j)
                last_edge = edge_with_j
            elif j_second:
                edge_set.remove((j_second[0], j))
                edge_with_j = (j, j_second[0])  # flip edge rep
                boundary.append(edge_with_j)
                last_edge = edge_with_j

            if edge0[0] == last_edge[1]:
                break

        boundary_lst.append(boundary)
    return boundary_lst[0]

#generating of random points
N = 1000
r = 1 - 2*np.random.random((N,2))
r_norm = np.linalg.norm(r, axis=1)
points = r[r_norm <= 1] 
plt.figure(figsize=(10,10))
plt.scatter(points[:,0], points[:,1], color='k', s=1)

# Computing the alpha shape
edges = alpha_shape(points, alpha=1, only_outer=True)
#order edges
edges = stitch_boundaries(edges)
plt.axis('equal')
edge_points = np.zeros((len(edges),2))
k=0
for i, j in edges:
    edge_points[k,:] = points[[i, j], 0][0] , points[[i, j], 1][0]
    k += 1
plt.plot(edge_points[:,0],edge_points[:,1])

#theoretical/expected edges
# xx = np.linspace(-1,1, 100)
# yy_upper =  np.sqrt(1 - xx**2)
# yy_lower = -np.sqrt(1 - xx**2)
# plt.plot(xx, yy_upper, 'r:')
# plt.plot(xx, yy_lower, 'b:')

Right now the cloud of points is black. Blue line is obtained from the algorithm above. computed edge

UPDATE: the starting and final points can be chosen the most left point (OR by hand it is no problem) I am expecting the following result: enter image description here

Upvotes: 3

Views: 226

Answers (1)

Iddo Hanniel
Iddo Hanniel

Reputation: 2056

You can use the fact (see the scipy.spatial.Delaunay documentation) that "for 2-D, the triangle points are oriented counterclockwise". Therefore, the outer edge points constructed in the Alpha shape will always be oriented counterclockwise i.e., the inner side of the shape will be to their left (if this was not certified in the documentation we could have checked it ourselves and flipped if necessary, but here there is no need).

This means that the points along the polygon between the leftmost point and the rightmost point are the lower hull, and the points between the rightmost and the leftmost are the upper hull. The code below implements this idea.

min_x_ind = np.argmin(edge_points[:, 0])
max_x_ind = np.argmax(edge_points[:, 0])
if min_x_ind < max_x_ind:
    lower_hull = edge_points[min_x_ind:max_x_ind+1, :]
    upper_hull = np.concatenate([edge_points[max_x_ind:, :], edge_points[:min_x_ind+1, :]])
else:
    upper_hull = edge_points[max_x_ind:min_x_ind+1, :]
    lower_hull = np.concatenate([edge_points[min_x_ind:, :], edge_points[:max_x_ind+1, :]])

and the result can be visualized for example with:

plt.plot(upper_hull[:,0],upper_hull[:,1], "r", lw=3)
plt.plot(lower_hull[:,0],lower_hull[:,1], "b", lw=3)

enter image description here

Upvotes: 1

Related Questions