Lattice fitting/regression

Question

I have a bunch of measured points in an image (as shown in the picture below), and I want to find the best square-sized lattice to fit those measured points, such that the sum of distance of each measured point from its closest lattice point reaches a minimum. Note that the measured points are not necessariy "complete" with regard to the ideal lattices points, meaning that some of the ideal lattice points may not have corresponding measured points (actually I also want to use the ideal lattice to "recover" the missing measured points).

It seems the square-sized lattice can be determined by the following parameters:

• pitch (horizontal pitch == vertical pitch, since it's a square-sized lattice)
• angle
• an origin point (seems it's not unique, we can just pick any ideal lattice point as the origin)

But I stuck at formalizing the distance between a measure point and its closest lattice point...

A picture with measure points:

So far I find two related posts:

1. Lattice fitting to points This one is almost the same as my problem, and the "solution" given by the poster himself should be viable for me too. But I want to find a "better" one.

2. The other one Matrix least squares with rows instead of columns seems have a nicer solution, but the problem is a little bit different.

Answer 1

As I did not find a general solution to my problem, I tried to find a solution for my specific problem at hand. Compared to general questions I asked, what's specific in my problem are that:

• Although some points may be missing from the measured data, there is no "outlier" points in my measured points, meaning that all measured points has a corresponding lattice point (plus a little bit "noise").
• The rotation is limited in a small range, e.g. (-5,5) degree
• The number of columns and rows is also limited within a known range, e.g. there are [15,16] columns; and there are more columns than rows.

The idea of my current solution is roughly:

• Project all points to a cartesian coordiantes rotated by a angle $\theta$
• Clustering (segmenting) the X coordinates of all measure points on the rotated coordinates. the clustering algorithm also gives a measure of "withiness" of each cluster (i.e., how close the points are within each cluster).
• For the optimal $\theta$, the sum of the "withiness" of the points in each cluster reaches a minimum

I chosen the R package Ckmeans.1d.dp for clustering the X coordinates, and use Newton-iteration for find the optimal $\theta$.

So far it seems it works for me.

For one of my data set, the $\theta$ vs. total withinss plot shows that it's a "good" function:

The code snippet in R:

    ```{r}
epsilon <- 0.000001  # stop whenwhen f' <= epsilon
h <- 0.001           # step, in degree

# theta: the angle in degree
# x: is the original vector
# k: is the number of clusters 
f <- function(theta, x, k) {
  radian <- theta / 180 * pi
  x1 <- x * complex(1, cos(radian), sin(radian))
  Ckmeans.1d.dp(Re(x1), k=c(k,k))$tot.withinss
}

f.pitch <- function(theta, x, k) {
  radian <- theta / 180 * pi
  x1 <- x * complex(1, cos(radian), sin(radian))
  res <- Ckmeans.1d.dp(Re(x1), k=c(k,k))
  (res$centers[k] - res$center[1])/(k-1)
}

# f'
f1 <- function(theta, step, x, k) {
  f0 <- f(theta, x, k)
  f.h <- f(theta+step, x, k)
  (f.h - f0) / step
}

# f"
f2 <- function(theta, step, x, k) {
  f0 <- f(theta, x, k)
  f.h <- f(theta+step, x, k)
  f.2h <- f(theta+2*step, x, k)

  f1.0 <- (f.h - f0) / step
  f1.h <- (f.2h - f.h) / step

  (f1.h - f1.0) / step
}

theta0 <- 0
while(T) {
  tmp1 <- f1(theta0, h, A, k)
  tmp2 <- f2(theta0, h, A, k)
  if (abs(tmp1) < epsilon) {
    print(c(theta0, tmp1, tmp2))
    break
  }
  theta0 <- theta0 - tmp1/tmp2
  print(c(theta0, tmp1, tmp2))
}

# now we have the angle, calcuate the pitch
print(paste("theta=", theta0, "degree; pitch=", f.pitch(theta0, A, k)))
```

The output of the code:

[1]   -0.2044698 1743.6150540 8527.4959049
[1]   -0.20446626   -0.02982388 8527.71416073
[1] -2.044663e-01  5.844214e-10  8.527714e+03
[1] "theta= -0.204466264230956 degree; pitch= 116.161669012395"