# How to cluster parts of broken line made of points?

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?

• This appears to be a minor generalization of the problem discussed at stats.stackexchange.com/questions/33078. – whuber Nov 7 '18 at 15:08
• @whuber yeah maybe i should highlight the fact that i would really need to use a clustering algorithm such as the ones offered by scikit-learn or tensorflow – Federico Gentile Nov 7 '18 at 15:51
• One problem with your data is that in different lines constituting your "zigzag" density of points is different. Must it be so in your data? – ttnphns Nov 7 '18 at 16:21
• ...in other words, your data consists of clusters-to-be of very different densities. – ttnphns Nov 7 '18 at 16:23
• Please clarify what you mean by "supposed to be continuous:" does that imply you are given a sequence of points rather than an unordered set of them? How accurately are these points known--is there perhaps some random error in them? – whuber Nov 7 '18 at 16:42

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.

• Thanks a lot for the answer! it is very clear!! only one more request if possible: would you be able to point some python library able to achieve such transformation? So far i just found links to CV2 but it's more suitable for image processing while my problem is a bit different – Federico Gentile Nov 7 '18 at 19:41
• I don't think you need any image processing software--I didn't use any for the example. I'll post the code I did use. – whuber Nov 7 '18 at 19:43
• CASH is a clustering algorithm based on the Hough transform extending this to multiple dimensions. I would consider piecewise linear regression here though, not clustering! – Has QUIT--Anony-Mousse Nov 18 '18 at 8:37
• @Anony-Mousse I considered that when first reading this question. Since the "lines" are given parametrically, they do not necessarily represent a response to a regressor: for instance, there's nothing to prevent them from being vertical. A regression method would fit a single point to a vertical stack of points, whereas that's not what is being asked for here. That's why I ended up recommending the Hough transform technique. – whuber Nov 18 '18 at 17:50