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

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

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.

• A problematic issue is that there can be multiple solutions. When you rotate the rectangle some degrees it still fits but half the grid points are not linked to a data point. So you need some prior information that allows to select which option is best. How many data points are typically missing? Can it be more than half? Can you have multiple data points per gridpoint? Jul 8 at 8:15
• With the random samples that I am generating, I randomly eliminate around 30% of the points. I'm not expecting it to be that high with the actual data where a typical dataset might be around 150 points. There should only be one point per gridpoint, but if the points remain in the group of the same gridline, it shouldn't matter with my current algorithm. Jul 8 at 8:30
• Could you elaborate on what you mean by "not necessarily uniform"? Would that mean the spacings in the x direction are the same, the spacings in the y direction are the same, but the two values could differ, as in a city street grid? Or would it mean the spacings in one or both coordinate directions can vary from one block to the next? If it's the latter, your problem might be practically insoluble except when the variance of the noise is tiny. Do you only have point coordinates, or would you also have street name (or their analogs) for each point? Are the errors independent or not?
– whuber
Jul 8 at 16:38
• Overall, so much potentially useful information is lacking here that we are a little hamstrung. If you would disclose the nature of these data--what they represent, how they are obtained, and how you process them to produce coordinates--we likely could be much more helpful.
– whuber
Jul 8 at 16:39
• Are the gridlines supposed to be equal distance? Jul 8 at 18:07

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]]; *)


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.

• Impressive. This is very creative. I'll have to take a closer look. Do you think a Hough transform could be used without using an image? Otherwise, I would do everything using data points until you need the Hough transform. This would eliminate reading, adding 180° image and downsizing. You could make your own image. It might look like this in Python: img=np.zeros([1000,600,3], dtype=np.uint8); img.fill(1); for pt in pts: cv2.circle(img, pt, 3, (255, 255, 255), -1). Jul 9 at 1:26
• @Jakub, I imagine that this transform can also be used without using an image, and apply the transform to you set of points. Jul 9 at 9:02
• This method is actually similar to the methods that rotate the image and look at the projection of the points onto the x and y axis. Those projections onto the x and y axis could also be turned into an image when we plot the coordinates as function of the rotation, and then look for the places with a high density. The clever part is that by using an image, it is easier to find the right angle. (When I was thinking about the problem I thought of looking at each vertical line in the plot of the Hough plot and try to detect clusters, e.g by k-means or a fitting a density curve with a kernel) Jul 9 at 9:16
• @Sextus Please follow the link to my previous answer (in the second line of this post) for an abstract. The advantages of using an image are extensive: the ability to use efficient methods based on the Radon Transform, to blur and smooth things, etc. You do not need a high-resolution image, provided you bin the counts rather than just the presence of points. A cell size noticeably, but not hugely, smaller than the minimum grid spacing ought to work. You do want to use square cells to enable efficient, accurate rotations.
– whuber
Jul 9 at 13:26
• @SextusEmpiricus Fair enough: for the reader's convenience, I have now embedded the second paragraph of that post at the beginning of this one.
– whuber
Jul 9 at 14:08

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.

• The number of rectangles creates a non-continuous cost function, which is difficult to fit with optim, and needs some special attention. One way could be to use with distance between grid lines (and skip the parameters for the size of the rectangle) but then there might still be a problem of local Optima. Jul 8 at 5:54
• Sounds good, but I don't know how I would determine the initial values of the seven parameters. Any recommendations? I'm guessing the parameters you're thinking of are x, y, a, b, theta, m, n where (x, y) are the coordinates of a point; (a, b) are the lengths of the sides of the rectangle; theta is the angle of the grid; and (m, n) are the number of rectangles in both directions. What if it's not a uniform grid? For the best fit, it's unlikely to be a uniform grid. Jul 8 at 6:00
• @SextusEmpiricus Yes, or one tries some continuous relaxation. But, practically, it would probably be easiest to just try some fixed choices for the grid. Jul 8 at 6:04
• @Jakub For initialization, choose those roughly from your data. This depends on what parameters you choose. If the grid is not uniform, then, depending on the kind of nonuniformity, you need more parameters to describe it. Jul 8 at 6:09
• @frank I don't know how to choose the initial values of the parameters from my data. I could provide a sample dataset, if that would help you explain it to me. There could be different distances between all the grid lines is the kind of nonuniformity that I was thinking. Jul 8 at 6:21

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

• Any row or column could have only one point and this includes the first and last row. My initial starting point is the point with the min or max y-coordinate depending which has a slope closest to zero with one of it's nearest neighbors. I'm not sure what you meant in the 2nd paragraph about rotating a rectangle. I'm currently not rotating any rectangles or finding distances to rectangles. I find distances to a line which I expect to be the first or last row. I will look at more of this tomorrow. Thank you. Jul 8 at 8:21
• @Jakub your algorithm sort of goes through the points line by line and finding clusters of points. The alternative is to find the parameters (describing the rectangular grid) by optimizing a cost function that uses the distance between data points and grid points. The trick for this optimization is to seperate the distance function in a horizontal and vertical part. This is a bit similar to what you do by fitting horizontal and vertical lines seperately, only it is more convenient to do this fitting when you rotate the points. Jul 8 at 9:45
• What do you mean by "marginal projection of the points"? Jul 8 at 15:22
• @SextusEmpiricus I'm open to an optimization vs the OLS that I'm using, but speed is important. Maybe opt can be faster, but usually not. I'm not so sure that the rows and columns are independent (separate) in my algorithm. My OLS model does them at the same time, not separately. I'm not sure if this is being pedantic, but I'm not clustering. I'm sorting the pts by distance and segmenting them by range. This is pt-by-pt, because the lines (groups) are not known. I thought this would be faster than clustering because it's one-dimensional and it's unlikely that I will see more than 200 points. Jul 8 at 21:58
• @Jakub, I had the impression that you were scanning somehow through the points based on distances with an initial reference line (the segmentation that you use is not clear to me and this is how I imagined it). Also the OLS model is not clear to me. I can see how you could fit parallel lines to different groups, but it is not clear to me how you can ensure the lines to be perpendicular. (and I guess now that you are not looking for a regular grid?)... Jul 8 at 22:31

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:

• Did you visually locate the local maxima for g and h? I would have to automate finding them and I'm not sure how to do this. Can this be modified to fit a non-uniform grid? Jul 17 at 7:38
• I wasn't sure what you mean by non-uniform grid. The step along one axis doesn't have to have the same length as the step along the other. But all the steps along one axis need to be the same. So the grid must be a subset of points of the form mu+nv for fixed u and v (integer m and n). So might not work for you. Finding the local maximum entails a brute force search through pixels using whatever a priori information you already have (eg. no more than 10 degrees from an axis, the two maxima need to be around 90 degrees apart). Jul 19 at 0:06
• Being able to handle unequal step-sizes is important for my objective. The example I gave did have equal step-sizes for each axis, so I can understand the confusion. However, the noise alone made for unequal divisions. The result shows the two lines on the right are good fits and the two lines on the left are not. The first gridline is to the left of the points and the second gridline is to the right of the points. I plotted it with lines which made it more clear. A parallel slopes OLS can handle that, but it has the weakness of having to group the points where your solution does not. Jul 19 at 3:11

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

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

• I tried to use k means as well but it becomes problematic with larger amounts of points and also you would need to know the number of clusters in advance. Jul 8 at 18:06
• Interesting. In what way does it become problematic with larger amounts of points? I agree that figuring out the number of clusters can be tricky, but I'm sure some of the automatic methods would give a good initial guess. Something like increasing the number of clusters one by one until two cluster centres are too close might be enough for a first pass. There is a chance the maximum number of gridlines could be known apriori from when the data was gathered. Jul 8 at 18:13
• when I used a 10 by 10 grid and fitted ten clusters, then I did not get 10 evenly spaced clusters. The algorithm placed 3 clusters around a single center and missed two centers. This might also be related to optimizing of the k-means clustering, I did not dive into it deeper, Jul 8 at 18:17
• Hmm, thats a disappointing result. I tried using MeanShift instead, and out of the box with no tuning it picks out a decent 5x4 grid, with the last gridline being inbetween two sets of points. Adjusting the parameter bandwidth for it varied between too many lines and not enough. A rough trial of bandwidth = values.size * 2 seems to do the job, finding the same 6x4 grid that is pictured above automatically. I'm sure there is a very suitable clustering algorithm that the buck can be passed to. Jul 8 at 18:46
• Results look great! Thanks! I timed your code and got 0.037 seconds for this example and 0.068 seconds for a 10x15 grid example; so I don't see significant scaling issues. My code does the same examples in 0.018 and 0.029 seconds, respectively. For my purpose, it is not likely the grid will be any larger than 12x16. My algorithm does not require the values of the number of rows and columns or the angle of the grid. I use a simple loop to segment because its 1-dimensional and I thought a clustering algorithm might be too costly. I would have to compare success rates to see if this works better. Jul 8 at 20:31