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.

Because the streets are on an irregular orthogonal grid, the Fourier Transform solutions--although clever--will likely fail.  The following approach exploits two ideas:

*

*Using the Hough transform to identify sets of points that tend to line up.


*Exploiting the assumed orthogonality of the grid to augment the data.
The link explains and illustrates the Hough transform and the Mathematica code I am using.  Here is the gist of it, as extracted from the second paragraph:

... the Hough transform depicts sets of lines. A line in the
plane can be parameterized by its slope, $x,$ and its distance, $y,$ from
a fixed origin. A point in this $x,y$ coordinate system thereby
designates a single line. Each point in the original plot determines a
pencil of lines passing through that point: this pencil appears as a
curve in the Hough transform. When features in the original plot fall
along a common line, or near enough to one, then the collections of
curves they produce in the Hough transform tend to have a common
intersection corresponding to that common line. By finding these
points of greatest intensity in the Hough transform, we can read off
good solutions to the original problem.

Please refer to the rest of that post for further explanation and code.  Here, I'll simply take you through the workflow with your example (the dots in the left image).
Import the image and crop it.
img = Import["https://i.stack.imgur.com/uAjLw.jpg"]
i = ColorNegate[Binarize[img]];
crop2 = ImageCrop[ImageTake[i, {25, 700}, {1, 1050}], 940];

Combine the image with its 90 degree rotation.
(Optionally, combine it with all three rotations at 90, 180, and 270 degrees.)
I also downsample it for speed and because coarsening the resolution is a crude kind of windowed smooth that helps identify points that line up.  Various forms of additional smoothing will appear further below without comment.  The amount of smoothing is a matter of trial and error but good results are quickly found.
crop3 = ImageResize[crop2, 80];
(* crop3 = ImageAdd[crop3, ImageRotate[crop3]]; *)
crop3=ImageAdd[crop3, ImageRotate[crop3, \[Pi]]] 


Notice how much more clearly the orthogonal geometry of the street intersections comes through in this combined image.
[Edit: I made a mistake here.  The commented-out line performs the 90 degree rotation; the final line rotates by 180 degrees.  I intended to comment out the last line only.  The result of this error is to produce an inferior solution--the one shown below--but it still illustrates the ideas, so I have left the figures as is.]
Compute the Hough transform.
hough2 = Radon[crop3, Method -> "Hough"]  // ImageAdjust


Smooth the Hough transform to help detect peaks.
blur = ImageAdjust[Blur[ImageAdjust[hough2, {1, 0}], 5]]


Identify the "morphological components" (isolated high points).
comp = MorphologicalComponents[blur, 0.85]  // Colorize


Find the angles associated with component centroids.
width = ImageDimensions[blur][[1]];
slopes =  Module[{x, y, z}, ComponentMeasurements[comp, "Centroid"] /. 
          Rule[x_, {y_, z_}] :>  Round[((y - 1/2)/(width - 1))  180., 0.1]
  ]

{101., 93.4, 7.1, 93.4, 101.}

The output varies a little depending on choices made previously, but the variation is consistent with the amount of noise apparent in the original image.
Rotate the original to make the grid appear isothetic (parallel to the coordinate axes).
j = ImageRotate[crop2, -slopes[[1]] / 180 * \[Pi]];
ImageResize[j, 150]


Smooth row and column sums to plot the marginal distributions of intensity.
"Streets" run left to right; "avenues" run straight up and down.
ListPlot[ListConvolve[{1,2,3,2,1}, (Plus @@ #)& /@ ImageData[j]], Joined->True, PlotRange-> {Full, Full},
PlotLabel->"Streets"]


You can post-process the locations of the peaks to identify the coordinates of the "street" centerlines.  The values along the horizontal axis are pixels indexing the cropped image I have been using.
Here is a plot for the transpose of the image data (the code is essentially the same):

Standard errors for these location estimates can also be obtained from the half-widths of the peaks in these "spectra."
At this juncture we have obtained the angle of the street grid and locations of all identifiable streets and avenues in the image.  Further processing will depend on the data structures you need and what you intend to do with this information.

I consider this to be a great result, because it has succeeded in identifying the numbers of streets and avenues and estimating the resulting $1+4+6=11$ parameters with appropriate precision from just $17$ rather uncertain point locations.
This method still works without the second idea (of combining the image with rotated versions of itself), when it can identify gridlines on a non-orthogonal street system, but it requires more data to achieve the same level of precision.
