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.

This answer is not a complete answer to the exact question but I leave it as it does explain some principles and works well in another related situation.

*

*This answer is based on the assumption that the grid is regular (equal distance between grid lines). In the solution this assumption is necessary when we apply the fit with a sine wave to detect the clusters of points.
The advantage of working with this assumption is that it works well for large and  sparse matrices. (but of course it only works of the grid is truly regular). An example is the image below that is made with by turning the problem below a bit more difficult by having the standard deviation of the noise at 20% of the distance between 30 by 20 grid lines and by having only 10% sampled points per grid point.



*The principle of the answer still stands. It is based on finding the clusters in a rotated image (one could use other methods than the sine wave applied here) and find the correct rotation by an optimisation.
The disadvantage is that the optimisation requires a lot of computation. Performing clustering at different angles. The advantage is that the method works without assumptions.
An example of such assumption. If the envelope of your grid points is forming a rectangle (this is not clear from your description as you say that points are missing), then you could quickly find an initial starting point by minimising the size of the smallest enveloping rectangle around the points, which is easy to find using the minimum and maximum coordinates after rotation)


*This method is similar to the method with the Haugh transform that Whuber describes. The image below shows how it is similar. We rotate the image and compute the x and y-coordinates. Those are plotted in the image as function of the angle.

The trick is to find the angle where the points are most strongly clustered around single centers (which depict the grid lines). With the Haugh transform this is done by looking for points with high contrast/density. The method below does this a bit less efficient by making a seperate computation for every seperate angle.
The image also shows how it is much easier to detect the optimal angle by looking for the smallest bandwidth of the coordinates.

In the case that a fitting is performed by computing the grid and a cost function based on distances, then a trick to speed up would be to rotate the points and find the distances to a rectangle with fixed orientation, instead of rotating the rectangle and finding the distances to the points.
In this way you do not need to find a nearest neighbour among $N \times M$ grid points, and only look at distance to the $N+M$ grid lines. The squared distance of the distance to the grid point is the sum of the squared distance to the grid lines.
So a part of the algorithm could be something like

*

*Rotate the grid points with some angle $\theta$



*Fit horizontal and vertical grid lines seperately to $x$ and $y$ coordinates of the points.
This just needs to take into account the distance and offset-position of the grid lines. The number of grid lines (size of the square) you do not to worry about. Consider at this moment the grid to be infinitely large.
I did a quick search but could not find anything. I imagine that there should be some standard algorithms/functions doing this, e.g. some Gaussian kernel that is convoluted with a regular grid. Otherwise if you do it yourself, then you might use an autocorrelation function to find an initial estimate for the distance between grid lines. You could also fit a Gaussian kernel with some number of points (for which there are many programs available), and use the fitted points to figure out the grid.
In the code example below we fit a sine wave such that the following function is maximized $$\sum_{k=1}^n \sin\left(\text{offset}+2 \pi \frac{x_k}{\lambda} \right)$$ the advantage of the sine wave is that it allows simple algebraic computations to find the optimal $\lambda$ and $\text{offset}$.

After reading the answer by Dan Piponi, I realize that this is effectively computing a Fourier transform.


*Compute some cost function using the distances of the points to the vertical and horizontal grid lines (possibly the sine wave cost function of the previous paragraph could be used, but I imagine that using something else, like squared distance, might result in less variance of the error).
In the example below we rescale and shift the points again such that we get grid lines at $x=0,1,2,\dots$ and we can use round off function to easily compute the distance.

Then you repeat the above with different values of $\theta$ to find the optimum for that parameter

In the above graph there is a clear minimum at -5 degrees, which is how the data was generated.
One remaining problem is that for different angles there are also minima and these are due to the tendency of adding more grid lines than the true model (overfitting). You need some way to limit the total amount of grid lines. For this you need prior information.
The fit above is an improvement from the code below. It is using the scaling to add a penalty for adding too many grid lines (this will scale the image more strongly and increase the distances). Instead of computing the distance as
### compute distance
distance[k] = sum((ttx-round(ttx))^2*modx$l_max^2) +
              sum((tty-round(tty))^2*mody$l_max^2)
}

the above image uses
### compute distance
distance[k] = sum((ttx-round(ttx))^2) +
              sum((tty-round(tty))^2)
}

And the distance is computed in the scaled frame (where grid line distance are standardized to be 1) instead of multiplying with l_max to scale back to the original.
So the distance is relative to the grid line distance.

Example r-code:
### settings
set.seed(1)
nx = 10
ny = 10
theta = 5*2*pi/360
dx = 1.5
dy = 1.5
sigma_x = 0.1
sigma_y = 0.1
p_sel = 0.7

### generate data
kx = rep(1:nx,each=ny) - (nx+1)/2+0.3
ky = rep(1:ny,times=nx) - (ny+1)/2
x = kx*dx + rnorm(nx*ny,0,sigma_x)
y = ky*dy + rnorm(nx*ny,0,sigma_y)
rx = cos(theta)*x-sin(theta)*y
ry = sin(theta)*x+cos(theta)*y
sel = which(rbinom(nx*ny,1,p_sel) == 1)
sx = rx[sel]
sy = ry[sel]

plot(sx,sy, main = "original data")


### example of rotated data
alpha = -6*2*pi/360
tx = cos(alpha)*sx-sin(alpha)*sy
ty = sin(alpha)*sx+cos(alpha)*sy
plot(tx,ty, main = "rotated data")


### function to find a regular grid for 1d data
### fitting by maximizing a sine wave
### we search for the offset of the wave and the wave length
fitpoints = function(s,lmd_range, plt = TRUE) {

   ### initial search for wavelength

   cost = sapply(lmd_range, function(z) sum(sin(s/z*2*pi))^2 + sum(cos(s/z*2*pi))^2 )
   l_max = lmd_range[which.max(cost)]
 
   if(plt) {
      plot(lmd_range, cost, type = "l", main = "find optimal lambda \n distance between lines")
   }

   ### given optimal wave length 
   ### search for offset
   ### 
   ### this makes use of the linear sum of cosine and sine waves
   ### sin(x+b) = cos(b)sin(x) + sin(b)cos(x)
   
   ss = sum(sin(s/l_max*2*pi))
   sc = sum(cos(s/l_max*2*pi))

   M = matrix(c(ss,sc,sc,-ss),2,byrow=TRUE)
   factor = solve(M) %*% c(sqrt(ss^2+sc^2),0)
   xs = seq(min(s),max(s),0.01)

   if(plt) {         
      plot(xs, factor[1] * sin(xs/l_max*2*pi) + factor[2] * cos(xs/l_max*2*pi), type="l", main = "sine wave with highest values for points")    
      points(s, factor[1] * sin(s/l_max*2*pi) + factor[2] * cos(s/l_max*2*pi), pch = 19)
    }

    angle = Arg(complex(real=factor[1],imaginary=factor[2]))
    offset = -angle/2/pi*l_max+0.25*l_max

    ### return result    
    return(list(l_max = l_max, offset = offset))
 }

modx = fitpoints(tx,seq(0.25, 5, 0.01))
mody = fitpoints(ty,seq(0.25, 5, 0.01), plt = 0)

ttx = (tx-modx$offset)/modx$l_max
tty = (ty-mody$offset)/mody$l_max

plot(ttx,tty, pch = 21, main = "rescaled and shifted data \n ready to get distance to lines by using roundoff function")
for (i in -30:30) {
   lines(c(i,i),c(-100,100),lty=2)
   lines(c(-100,100),c(i,i),lty=2)
}

### find optimal angle

angle_range = seq(-20,10,0.2)*2*pi/360
distance = angle_range*0 ### create an empty vector

for (k in 1:length(angle_range)) {
   alpha = angle_range[k]

   ### rotate
   tx = cos(alpha)*sx-sin(alpha)*sy
   ty = sin(alpha)*sx+cos(alpha)*sy
   ### optimize lines
   modx = fitpoints(tx,seq(0.5, 5, 0.01), plt = 0)
   mody = fitpoints(ty,seq(0.5, 5, 0.01), plt = 0)
   ### shift and rescale
   ttx = (tx-modx$offset)/modx$l_max
   tty = (ty-mody$offset)/mody$l_max
   ### compute distance
   distance[k] = sum((ttx-round(ttx))^2*modx$l_max^2) +
                 sum((tty-round(tty))^2*mody$l_max^2)
}

plot(angle_range*360/2/pi, distance, type = 'l' , main = "distance for grids fitted with different rotation")
points(angle_range*360/2/pi, distance, pch = 19)


Sidenote: The above does not compute a beginning or end of the rectangle. However, after finding $\theta$ and the grid lines, you can find out easily where the rectangle ends and starts (I would not place the end and start of the rectangle in the optimisation, these parameters do not change the cost function and just slow down the fitting).
