# Compute grid orientation, spacings, skew and origin from a set of 2D points

I have a collection of 2D points that are originated from a regular grid (rectangular grid most of the time, but possibly skewed as well). Each point can be located on the grid from its Row and Column numbers or i, j indices in a matrix-like fashion.
However, it is possible that a majority of these points have shifted from their original grid location to a final location but there is no guarantee this final location sits exactly on a grid node.
Shifted points should get their row, column numbers updated to reflect their final location (or rather the closest grid node) but there is no guarantee for this. It is possible they have retained their initial indices.
Finally, the overall shape of the points set cannot be assumed to be rectangular-like as rows may very well be of unequal length.

Here's an example:

I know the points x, y coordinates as well as their Row, Columns numbers and I am trying to retrieve the following information:
- the grid aspect: di, dj (row height and column width)
- the grid orientation (rows may not be aligned with the x-axis)
- the grid origin (the x, y coordinates of the cell with i=0, j=0)
- the grid angle (if the grid is not rectangular)

di, dj are fairly easy to compute as I know the row and column numbers, I can compute the distance between consecutive points on each row and the median of all intervals will generally give me very accurate results.
My first attempt to compute the orientation was similar: by computing the orientation of segments between consecutive points and smoothing out changes of orientation I got some very accurate results as well there.
Computing the origin was cumbersome however but I recently discovered Procrustes analysis (see links below) and by applying a procrustes superimposition to the x, y coordinates and i, j indices, I'm not only getting the grid orientation but also a full set of parameters (Rotation, reflexion, scale and translation) to directly transform i, j indices to x, y coordinates so it is very straight forward to compute the grid origin but whilst I'm still getting very good estimates for the orientation, the origin precision is not very satisfactory despite my attempt at detecting points that have shifted to far out from a grid node and removing them from the computation.

So I guess my question is: is there a more appropriate way to achieve this ?

My problem seems related to the question measuring the regularity of a grid and possibly Detecting grid orientation and properties as well but I don't understand how to apply the structure factor hinted at in the later to my problem.
Other similar questions are these ones:
- Detect grid in 2D plane
- Lattice fitting to points
- Lattice fitting/regression
- Fitting data to a square lattice (discrete points by multiple parameters)
Unfortunately, non of these questions ended up being answered by other than their respective OP...

References:
+ Procrustes analysis
+ Orthogonal Procrustes problem