Estimating absolute position from approximate relative distances

Viewed 209

My problem is the following:

I have a list of relative XY-distances between points scattered in a plane that have been measured with some amount of error. I have defined one point as being the origin and I would like to estimate the absolute position of all the other points in that frame of reference. Each the point have at least two paths that links them to the origin and due to the measurement error the absolute position I would get from each path is different. I would like to know :

  • does this problem have a name ?
  • do you know of an algorithm to make a good guess at the absolute positions ?

The rest of the text is only details that you may skip over if the problem is already clear for you.


The context of this problem is that I want to build a composite image from a lot of smaller images that have been taken in sequence but at irregular positions. I've computed a first set of distances by computing the cross-correlation between each image in the sequence and thus made a first rough estimate of the absolute positions. There is a significant amout of drift in these absolute positions and images that should be perfectly overlapping are not. Using this I then computed the cross-correlation between each overlapping image so I now have a network of distances between the images and I hope to be able to use that to refine the absolute positions and get a better final image.

Here is the first composite image I made: enter image description here

The base images are sligtly transparent and put on a white background. Several images stacked at the same place = darker image, no data = white. The colorfull snake is the sequence I used to compute the first absolute positioning, each segment illustrate the distance measured between the image it is on and the next one, and it starts in the bottom left corner. In the bottom left corner we can see two images that cover the same corner feature but that are placed with a 200 pixel error in both x and y directions.

I tried to use automatic stitching tools but they all failed because of the general lack of contrasted distinctive features on many images.

2 Answers

I do not know of a name for this problem, but https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html can find the solution that minimizes the sum of the squares of the errors.

You already have fixed the first point at the origin. You can insist that the second point is on the x-axis. That then leads to a list of points like this: (0, 0), (0, y1), (x2, y2), (x3, y3), ..., (xn, yn). Now you simply write a function that takes [y1, x2, y2, x3, yx, ..., xn, yn] and returns the sum of the squares of the errors. Optimize it. The BFGS method is probably a good choice for that.

To simplify the description, I will assume that the points are represented by complex numbers.

Your problem can be formulated like this:

We have n+1 points, the one arbitrarily set to (0,0) and n *unknown" points, and you have m noisy observations, m being (much) greater than n.

An observation correspond to a (slightly) erroneous estimation of the difference between to points i and j:

Y[k] = X[i] - X[j]

Then, the set of relations can be represented as:

Y = A X + N

Where X is the vector of unknown points, vector Y correspond to the observations, A is the matrix defining the observations and N represent the errors of the observations.

A first approach is the LS (Least Square) one.

By using the Moore-Penrose pseudo-inverse, we can get the LS estimation by:

X1 = (A^H A)^(-1) A^H Y

Where A^H is the Hermitian transpose of A. As here A is real, it corresponds to the transpose of A. This estimation minimizes

min norm of (AX - AX1)

Another approach is to assume the N is an additive Gaussian noise of variance s2. The we can get the MMSE/Wiener estimate by:

X2 = (A^H A + s2 I)^(-1) A^H Y

Where I is the identity matrix of size n. This estimates minimizes:

min E (norm (X2 - X)) = min MSE (Mean Square error)

where E(.) represents the statistical mean (expectation). This is a stochastic method (we minimize an expectation), on the contrary of the LS method.

Sometimes, we cannot assume that the noise is Gaussian, or we don't know the variance. Moreover, it may appear that the LS estimate is not totally stable, depending of the singular values of matrix A. Then, a third method is sometimes used:

X2 = (A^H A + lambda I)^(-1) A^H Y

Where lambda is a (small) relaxation coefficient to avoid the possible instability of the LS method.

Note: if matrix inversion is too complex, and if the observation errors are rather small, it is certainly possible to use iterative methods to improve the estimations.

Related