Fitting an Orthogonal Grid to Noisy Points

Question

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

1. 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.
2. Calculate the distances between this reference line and the remaining points.
3. Segment points into groups w.r.t. the calculated distances based on group range. Each group should represent one gridline.
4. Repeat steps 1 to 3 for the 90° rotated set of the same points. Combine results.
5. Create a parallel slopes OLS model to determine linear equations for the gridlines.
6. Rotate back the gridlines for the previously rotated points to their original orientation.
7. 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.

Answer 1

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

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

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

Answer 2

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.

Answer 3

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

Answer 4

We can turn the problem of finding the points in a noisy grid into the problem of finding the peaks of a smooth function, which can be much easier.

If your points are $p_i$, define the function $f$ of 2D vector $k$ by $f(k)=|\sum_{i=1}^n\exp(i k\cdot p_i)|$. For your example set of points you'll get the function whose "heat map" looks like

(You can think of this as the Fourier transform of the function defined by placing Dirac deltas at the points.)

There's a 3x3 grid of (bright white) local maxima in the middle with 0 at the centre that form a parallelogram. The vectors giving the midpoint of the right side and the top side give a pair of vectors, call them $g$ and $h$. In this case we have roughly $g=(0.039,0.006)$ and $h=(-0.005,0.039)$. These define the "reciprocal" lattice to the one you want. You can get the original lattice by taking the inverse of the matrix whose rows are $g$ and $h$ and scaling by $2\pi$.

Here's some Mathematica code that plots the original points along side the points in the new perfectly regular lattice

This doesn't give you the 'phase' of the points so you'll have to slide the lattice around until it fits.

Both the problem of finding bright peaks and sliding around a fixed lattice to get a fit are easier than the original problem. So this may make things easier for you.

On the negative side: this method is robust up to a certain amount of noise/deformation and then it quickly falls apart so you need to start with your worst examples and see if this approach can handle them. For example, I chose not to use this method for finding grid-points in a lens calibration method because even though the image looked grid-like, there was no actual regular grid that fit over the entire image due to the lens distortion.

If it's ambiguous which peaks to use you can, of course, write code that tries more than one pair of vector candidates.

Update: I think you don't need to do any shifting yourself. We can use the phase of the sum of the exponentials to extract the shift. Again, here is some Mathematica code:

Answer 5

As an adjustable approach, you could try

1. Rotate all the points by a fixed angle (around the origin or centre of mass, whatever suits).
2. Project the rotated points to the x-axis and y-axis. This is just the x-values and y-values of the rotated points.
3. 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.
4. The centres of the clusters says where the grid lines should go.
5. 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.