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.

How about this approach:
Your grid can be described by seven parameters: Consider the lower left rectangle in your grid. You need five parameters to describe the four corners of this rectangle, and two parameters for the number of rectangles in both directions.
Next, you design a cost function that takes as input those seven parameters and your data (the points) and computes for each of your data points the squared distance to the nearest grid point. The sum of those squared distances would be returned as cost. Finally, you might want to add an extra term to penalize the grid vertices count: Imagine a grid with extremely many vertices, then there would be a vertex near any of your data points, reducing the cost to almost zero, but that is probably not what you want.
Finally, find the minimum of your cost function and the belonging seven parameters, e.g. with the R function optim, if you cannot find anything more appropriate.
