How to find the slope of certain lines matching random points

Question

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:

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).

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?

Answer 1

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:

Note.
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)

Answer 2

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
        img.putpixel((i,j),(r,g,b))

thresh = 200
fn = lambda x : 255 if x > thresh else 0
r = img.convert('1')
r.save('/home/test/Desktop/sbVut.bmp')

"""
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[0].set_axis_off()

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[1].axis('image')

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_axis_off()
ax[2].set_title('Detected lines')

plt.tight_layout()
plt.show()

And the results

Answer 3

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.

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