Fitting an Orthogonal Grid to Noisy Points I have a list of coordinates that are meant to form an orthogonal grid that could be rotated. The grid is not necessarily uniform. The rotation is not typically greater than 10°. The coordinates are noisy and some are missing. I want to fit an orthogonal grid to the coordinates and return a list of the grid's vertices.

I have an algorithm that I wrote in Python that was able to find the grid in the image above and many like it, but was hoping to find a more robust algorithm in Python or R.
Current algorithm

*

*Find a reference gridline to compare with the remaining points by
choosing two neighboring points in the first or last row with a slope
closest to zero.

*Calculate the distances between this reference line and the remaining points.

*Segment points into groups w.r.t. the calculated distances based on group range. Each group should represent one gridline.

*Repeat steps 1 to 3 for the 90° rotated set of the same points. Combine results.

*Create a parallel slopes OLS model to determine linear equations for
the gridlines.

*Rotate back the gridlines for the previously rotated points to their original
orientation.

*Calculate the intersection points.

I used NumPy for the data, KDTree from SciPy for segmenting, and ols statsmodels.formula.api for the parallel slopes model.
Example Dataset
Note: This is not the same points from image because I generated these randomly at a later time.
pts = [(104, 131), (240, 136), (580, 183), ( 88, 234),
       (396, 277), (199, 431), (367, 451), (534, 464),
       ( 29, 554), (171, 627), (342, 628), (493, 638),
       ( 10, 739), (144, 747), (138, 927), (472, 966)]



An Update in the TLDR Category[but some commentors seemly wanted more info]
Purpose
I was hoping to keep the question more general because I find it to be a interesting statistics problem. Keeping it general keeps the assumptions to a minimum. For example, my first solution assumed the grid would be uniformed and aligned with the axes of the coordinate system in a Euclidean plane (i.e. not rotated). I was wrong, so my solution failed in many cases. A more detailed question could also introduce biases and lead to going off on tangents (e.g. data collection).
That said, this algorithm is meant to determine a grid of objects in a digital image (no gridlines, just objects). The data points come from the results of contour or edge detection using OpenCV. Examples include an image with a grid of headshots, a street grid, a lattice in a SEM image, etc. The points are noisy or missing because different shadings can have significant effects on contour detection. Another example is using face detection in locating headshots. The headshot photos (rectangles) might line up on a grid, but the faces are not aligned or not detected.
The images typically don't contain more than 200 points, rotated more than 10 degrees, or missing more than 20% of points. The origin of the coordinate system is the upper-left corner of the image. I'm still keeping the assumption of orthogonality, but that could be an issue if an image is skewed.
Python Code
I originally asked this question on stackoverflow two years ago (July 17, 2020). I didn't receive an answer until three weeks ago and that was for converting my initial code based on lists — that I never even used — to numpy. I then posted my solution answering my long forgotten question after I had some renewed interest.
My code runs in 0.018 seconds for the example points provided, and 0.30 seconds for a 12x16 grid of points. It runs at around a 99% success rate for a Gaussian distribution of errors of around 3% of the avg x and avg y-coordinates. In other words, 1% of the time there will be a few points too far out on the tails of the distribution.
Issues
As with any model, problems arise when the residual noise level becomes too high. Regarding the 1% fail rate of my code, removing outliers would probably help since a few more missing points would not likely have a significant impact. Noise can make finding the correct reference line an issue which comes down to finding the correct angle of rotation of the grid. A possible improvement could be to segment the points w.r.t. slopes relative to their k-NN's. Noise or an extraneous point can also be an issue when segmenting the points w.r.t. to distance from the reference point in that it could create an additional gridline.

As an adjustable approach, you could try

*

*Rotate all the points by a fixed angle (around the origin or centre of mass, whatever suits).

*Project the rotated points to the x-axis and y-axis. This is just the x-values and y-values of the rotated points.

*Use k-means or a method of clustering points on a 1d line. K-means lets you specify the number of clusters, which here will correspond to the number of parallel lines in the non-uniform grid.

*The centres of the clusters says where the grid lines should go.

*Rotate the grid lines backwards by the same angle to line them up with the original points.

With this approach you can explicitly set the number of grid lines in each direction. If searching for a solution, the parameters to adjust are reduced to a float, and two integers (the fixed angle, the number of vertical lines, and the number of horizontal lines).
Here is an example solution from some hasty code. I eyeballed the angle and number of lines each way.

import matplotlib.pyplot as plt
import numpy as np
from sklearn.cluster import KMeans


def fit_grid_to_points(
    points, angle: float, vertical_lines: int, horizontal_lines: int
):
    rotated_xs, rotated_ys = rotate_points(points, angle)

    x_centres = find_centres(rotated_xs, vertical_lines)
    y_centres = find_centres(rotated_ys, horizontal_lines)

    # Rotate the centres back. Match up with the original alignment.
    x_centre_points = np.hstack((x_centres, np.zeros_like(x_centres))).T
    y_centre_points = np.hstack((np.zeros_like(y_centres), y_centres)).T

    rotated_back_x_centres = rotate_points(x_centre_points, -angle)
    rotated_back_y_centres = rotate_points(y_centre_points, -angle)

    return rotated_back_x_centres, rotated_back_y_centres


def main():
    pts = np.array(
        [
            (104, 131),
            (240, 136),
            (580, 183),
            (88, 234),
            (396, 277),
            (199, 431),
            (367, 451),
            (534, 464),
            (29, 554),
            (171, 627),
            (342, 628),
            (493, 638),
            (10, 739),
            (144, 747),
            (138, 927),
            (472, 966),
        ]
    ).T

    angle = -0.14  # radians
    vertical_lines, horizontal_lines = 4, 6

    res_xs, res_ys = fit_grid_to_points(pts, angle, vertical_lines, horizontal_lines)
    display_found_grid(pts, res_xs, res_ys, angle, vertical_lines, horizontal_lines)


def find_centres(values, number_of_centres: int):
    kmeans = KMeans(n_clusters=number_of_centres)
    # Reshape because this needs a dummy dimension.
    kmeans.fit(values.reshape(-1, 1))
    return np.array(kmeans.cluster_centers_)


def rotate_points(points, angle):
    """Rotate the points around the origin by `angle` radians.
    Modified from
    https://stackoverflow.com/questions/34372480/rotate-point-about-another-point-in-degrees-python
    """
    xs, ys = points[0, :], points[1, :]
    rotated_xs = np.cos(angle) * xs - np.sin(angle) * ys
    rotated_ys = np.sin(angle) * xs + np.cos(angle) * ys
    return rotated_xs, rotated_ys


def display_found_grid(points, res_xs, res_ys, angle, vertical_lines, horizontal_lines):
    xs, ys = points[0, :], points[1, :]
    plt.scatter(xs, ys)

    for x_centre in np.array(res_xs).T:
        plt.axline(x_centre, slope=1 / angle, color="black", alpha=0.6)

    for y_centre in np.array(res_ys).T:
        plt.axline(y_centre, slope=-angle, color="black", alpha=0.6)

    plt.title(
        f"Fitting a non-uniform ({horizontal_lines}x{vertical_lines}) grid to points\n"
        + f"using a fixed angle of: {angle:0.3} (radians)"
    )
    plt.show()


if __name__ == "__main__":
    main()

The quality of the solution will depend on the clustering algorithm. Here I used kmeans as a quick fix. A more carefully chosen algorithm might be smarter about ensuring minimum distances between clusters which will guarantee the lines aren't pushed too close.
