I have a large set of random points, $(x_i,y_i)$, which, for unknown reasons, seem to align along certain lines:


I want to calculate the slope of those lines (also their other parameters, but I'm asking about the slope here)

Separating the lines algorithmically, seems to be too complex. I guess that I would need a classification algorithm, and those never give stable results.

That's why I am only talking about the slope.

What I did:

I made a synthetic sample generating random lines:

synthetic sample

The average point $\color{blue}{(c_x,c_y)}$ is drawn in blue, and on this chart, the lines have $-45$° angle.

I have the intuition that the sum of all points distances to a line passing by the blue point, should be maximal/minimal when the line matches the slope I'm looking for (or the perpendicular slope).

Distance to normal line

I calculated the distance summed of all points to a line passing for the center, with angle $\alpha$, for different angles (from $0$ to $360$), and I got this chart:


I drawn the point lines at $-45$°, and the distance seems to be maximal just around the normal angle $+45$, so it suggest that I could use a maximization solver to find the slope normal to the lines, by maximizing the distance of all points to a line passing by the center (blue point).

It seems to work for different angles I tried, but I'm not sure if that's a correct procedure.

Also, the distance is not minimal when the slope $\tan(\alpha)$ is parallel to the lines.

EDIT: I just noticed that if the points are clustered in two distant groups, independent of the orientation of the lines, the maximum distance would be the one separating the groups. Maybe Fourier transforms would be able to detect the orientation of the lines?

  • $\begingroup$ Note that in the diagram where the points form lines with a common (negative) slope, you have logarithmic $x$ and $y$ axis scales. So, the points $(x_i, y_i)$ themselves do not form lines; $(\log x_i, \log y_i)$ do. The way I'd approach this problem is collect points $(\log x_i, \log y_i)$ into sets, joining each point to a set if it has a neighbor close enough. given a suitable distance limit, you'll get several sets, each for a specific line, plus a few random sets in the scattered region. $\endgroup$
    – Glärbo
    Feb 23, 2021 at 16:18
  • $\begingroup$ @Glärbo By $(x_i,y_i)$ I meant the log, as drawn. You are describing a clustering algorithm, but there is no general criteria to decide the "suitable distance". The approach I describe, may also fail because of large scale structure. I wonder if a Fourier transform calculated at all angles would detect the $\alpha$ in a more robust way. $\endgroup$
    – Raxi Ral
    Feb 23, 2021 at 16:53
  • 1
    $\begingroup$ The Fourier transform idea sounds promising. $\endgroup$ Feb 23, 2021 at 17:15
  • 1
    $\begingroup$ I suspect you'll find that the rotated 2D Fourier transform needs similar analysis (as the clustering approach): you'll end up looking at $\alpha$ that maximizes the difference in high frequency components between the two axes – but what is that high frequency limit or range? The line intervals are not constant, nor are the pairwise nearest-neighbor distances between points in the same line. This is why I would use a discrete clustering approach, obtaining a fit for each cluster. $\endgroup$
    – Glärbo
    Feb 23, 2021 at 19:07
  • 1
    $\begingroup$ You use the term "seems". How much distinguishable are the lines? Did you try a 2D Fourier Transform? $\endgroup$
    – Moti
    Feb 24, 2021 at 5:32

3 Answers 3


Using the Hough transform (see here) over a bitmap extracted from the original figure and with some help from a python repository. The process begins with the furnished image sbVut.png

import numpy as np
from skimage.transform import hough_line, hough_line_peaks
import matplotlib.pyplot as plt
from matplotlib import cm
from PIL import Image

To transform the sbVut.png image into sbVut.bmp

img = Image.open('/home/test/Desktop/sbVut.png')
(w,h) = img.size    
for i in range(w):
    for j in range(h):
        (r,g,b) = img.getpixel((i,j))
        r = 255-r
        g = 255-g  
        b = 255-b

thresh = 200
fn = lambda x : 255 if x > thresh else 0
r = img.convert('1')

Begins the Hough transform

image = np.array(Image.open('/home/test/Desktop/sbVut.bmp'))

# Classic straight-line Hough transform
(h, theta, d) = hough_line(image)

# Generating figure 1
(fig, axes) = plt.subplots(1, 3, figsize=(15, 6),
                  subplot_kw={'adjustable': 'box'})
ax = axes.ravel()

ax[0].imshow(image, cmap=cm.gray)
ax[0].set_title('Input image')

ax[1].imshow(np.log(1 + h),
         extent=[np.rad2deg(theta[-1]), np.rad2deg(theta[0]), d[-1], d[0]], cmap=cm.gray, aspect=1/1.5)
ax[1].set_title('Hough transform')
ax[1].set_xlabel('Angles (degrees)')
ax[1].set_ylabel('Distance (pixels)')

ax[2].imshow(image, cmap=cm.gray)
for _, angle, dist in zip(*hough_line_peaks(h, theta, d)):
    y0 = (dist - 0 * np.cos(angle)) / np.sin(angle)
    y1 = (dist - image.shape[1] * np.cos(angle)) / np.sin(angle)
    ax[2].plot((0, image.shape[1]), (y0, y1), '-r')
ax[2].set_xlim((0, image.shape[1]))
ax[2].set_ylim((image.shape[0], 0))
ax[2].set_title('Detected lines')


And the results

enter image description here


As the points are dense, the segments will be isolated by a connected-components algorithm (or contouring in OpenCV). No particular care to take.

For every blob, you can fit a line, but I guess that just joining the extreme points is enough.

enter image description here

Below, the angles formed with the vertical (using the minimum area rectangle) for the 50 largest blobs.

enter image description here

  • $\begingroup$ In fact, separating the lines algorithmically seems very easy. $\endgroup$
    – user65203
    Apr 2, 2021 at 20:03

The method employed by us:

  1. Convert the picture issued by the OP to Windows Bitmap format
  2. Load the BMP file into a computer program, developed for the purpose
  3. Map the color content onto real (double precision) numbers
  4. Define a zero level of this mapping; in our case white = - 0.2
  5. Define empirically limitations for the area that looks interesting
  6. Use a contouring module for making isolines at level zero
  7. Use all sort of limitations for distinguishing relevant data from the rest
  8. Calculate second order moments (variances) for the relevant isolines
  9. Take a look at the theory in the webpage Two-dimensional Moments
  10. Especially take notice of the formula $\;\tan(2 \theta) = 2 \sigma_{xy}/(\sigma_{xx} - \sigma_{yy})$
  11. Calculate mean value $\pm$ spread of relevant angles $\alpha=-\theta$ and Output
The numerical end result obtained is:
Alpha = 38 +/- 2 degrees
Picture produced:

enter image description here

A more robust method may be to consider instead the angles of the minor axes of the ellipses of inertia $\;\sigma_{yy}(x-\mu_x)^2-2\sigma_{xy}(x-\mu_x)(y-\mu_y)+\sigma_{xx}(y-\mu_y)^2=\sigma_{xx}\sigma_{yy}-\sigma_{xy}^2\;$ with the x-axis.

Free accompanying source-only software (Delphi Pascal) has been made available at the site: MSE publications / references 2021
Disclaimer. Anything free comes without guarantee :-(and without referee)


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .