How to cluster parts of broken line made of points?

Question

I am studying clustering techniques and i am pretty new at this topic. Here is my problem: I created a 5 lines which are made of points. This lines are supposed to be continuous and they look like this:

# Create data
n0, x0, y0 = 17, np.linspace(0,  11, n0), np.linspace(10, 3,   n0)
n1, x1, y1 = 35, np.linspace(11, 27, n1), np.linspace(3,  3,   n1)
n2, x2, y2 = 15, np.linspace(27, 35, n2), np.linspace(3,  15,  n2)
n3, x3, y3 = 4,  np.linspace(-1, 0,  n3), np.linspace(10, 10,  n3)
n4, x4, y4 = 10, np.linspace(35, 46, n4), np.linspace(15, 15.5,n4) 

# Plot data
plt.figure()
plt.plot(x0, y0, 'o', color='grey', markersize=1)
plt.plot(x1, y1, 'o', color='grey', markersize=1)
plt.plot(x2, y2, 'o', color='grey', markersize=1)
plt.plot(x3, y3, 'o', color='grey', markersize=1)
plt.plot(x4, y4, 'o', color='grey', markersize=1)
plt.show()

My goal is to use a clustering technique in order to be able to cluster each line that i have created in order to recognize the 5 different lines presented in the plot.

To do so i have opted for the GaussianMixture clustering algorithm which i thought it could be suitable for this sort of data distribution (a line could be seen as a very skewed distribution maybe). Here is what i wrote:

# Prepare data for clustering
X = np.concatenate((x0,x1,x2,x3,x4))
Y = np.concatenate((y0,y1,y2,y3,y4))
data = np.column_stack((X,Y))

# Cluster
from sklearn.mixture import GaussianMixture
gmm = GaussianMixture(n_components=5, covariance_type='full')
X_ = gmm.fit(data)
y_ = X_.predict(data)
print(set(y_))

And the output that i got is something like this (it changes every time i run the code though...):

# Plot clusters
plt.figure()
plt.plot(data[y_==0][:,0], data[y_==0][:,1], 'o', color='red')
plt.plot(data[y_==1][:,0], data[y_==1][:,1], 'o', color='blue')
plt.plot(data[y_==2][:,0], data[y_==2][:,1], 'o', color='lime')
plt.plot(data[y_==3][:,0], data[y_==3][:,1], 'o', color='orange')
plt.plot(data[y_==4][:,0], data[y_==4][:,1], 'o', color='cyan')
plt.show()

As you can see the colors should represent the 5 clusters (aka lines) that i have originally created but apparently the output is not what i want.

Could you please provide a better way to approach this problem? In case what i am doing is partially correct, what am i mistaken?

Answer 1

The data appear to be a sequence of $(x,y)$ coordinates and the question concerns how to cluster them into linear segments.

Solution

A simple method uses a brute-force (but efficient) version of the Hough Transform.

The basic idea is that every line in the plane can be parameterized by its distance to a fixed origin and the direction perpendicular to it. Each ordered pair of distinct points, $(x_i,y_i)$ and $(x_{i+1},y_{i+1}),$ determines a line and therefore may be encoded as a (distance, angle) ordered pair. Sequences of points lying on the same line can then be easily recognized as stationary sequences of (distance, angle) pairs. Simply cluster these.

For numerical stability, it is wise to use an origin adapted to the data--perhaps their barycenter.

There are many ways to proceed from this point. For instance, a kernel density estimate of the (distance, angle) pairs is tantamount to the Hough transform itself. After removing outliers, any reasonable clustering method for 2D Euclidean points should readily identify all the segments.

Remark

If the data do not have a definite sequence, compute the Hough transform for all possible pairs $(x_i,y_i)$ and $(x_j,y_j).$ This is more computationally demanding and will produce a much messier transform, but nevertheless any lines formed by more than a few of the points will correspond to the highest-density regions of the transform.

---

Example

The question presents a sequence of 81 points (not all distinct). Here are portions of the derived information described above:

            x         y id  Distance      Angle
1   0.0000000 10.000000  1  7.999634 -2.1375255
2   0.6875000  9.562500  1  7.999634 -2.1375255
3   1.3750000  9.125000  1  7.999634 -2.1375255
...
16 10.3125000  3.437500  1  7.999634 -2.1375255
17 11.0000000  3.000000  1       NaN  0.0000000
18 11.0000000  3.000000  2  3.703704 -1.5707963
19 11.4705882  3.000000  2  3.703704 -1.5707963
...
51 26.5294118  3.000000  2  3.703704 -1.5707963
52 27.0000000  3.000000  2       NaN  0.0000000
53 27.0000000  3.000000  3  7.812028 -0.5880026
54 27.5714286  3.857143  3  7.812028 -0.5880026
...
66 34.4285714 14.142857  3  7.812028 -0.5880026
67 35.0000000 15.000000  3  6.164931  1.7088024
68 -1.0000000 10.000000  4 -3.296296 -1.5707963
69 -0.6666667 10.000000  4 -3.296296 -1.5707963
70 -0.3333333 10.000000  4 -3.296296 -1.5707963
71  0.0000000 10.000000  4 -6.102943 -1.4288993
72 35.0000000 15.000000  5 -7.610268 -1.5253730
...
80 44.7777778 15.444444  5 -7.610268 -1.5253730

The regions of stationary (distance, angle) values are clear and the breaks between them are obvious.

Here is a plot of the (distance, angle) pairs, colored according to the original segments created in the question:

I applied a hierarchical clustering solution (using the hclust function in R with the default complete linkage) and used that to classify the original data (except for the last point). Here is the classification, shown with colors:

Code

I know this question is oriented towards Python, but I used R for the examples so I'm confident this code works :-). To do the clustering by line segment, just apply any suitable clustering algorithm to the Hough transforms of the data increments:

f <- function(x, origin=c(0,0)) {
  x <- t(t(x) - origin)                       # Center data at the origin
  dx <- rbind(apply(x, 2, diff), c(0.0, 0.0)) # Compute differences
  v <- cbind(-dx[, 2], dx[, 1])               # Compute the normal directions
  n <- sqrt(v[, 1]^2 + v[, 2]^2)              # Prepare to normalize the vectors `v`
  cbind(Distance=rowSums(x * (v/n)), Angle=atan2(v[, 2], v[, 1]))
}

The function f accepts an $n\times 2$ array of $(x,y)$ data and an optional origin. It returns the (distance, angle) coordinates for each increment from $\mathbf{x}_i = (x_i,y_i)$ to $\mathbf{x}_{i+1}=(x_{i+1},y_{i+1})$ (angles in radians). The formulas, in succession, are

$$\eqalign{ \mathbf{v}_i &= (\mathbf{x}_{i+1}- \mathbf{x}_{i})^\perp\\ \text{distance}_i &=\mathbf{x}_{i} \cdot \mathbf{v}_i / |\mathbf{v}_i|\\ \text{angle}_i &= \arg(\mathbf{v}_i). }$$

The $^\perp$ operator rotates a vector through 90 degrees via $(x,y)^\perp = (-y,x)$ and the $\arg$ operator is the argument (angle) of the complex number $x+i\,y,$ usually computed using some version of an arctangent function.