Mapping as an optimization problem

Mathematica example implementation

I am generally interested in maps of cities and places. During these winter holidays, I finally got to try out something I was wondering about for years: Assume you want to determine the location of certain points, for example the trees and buildings on your property. Since all these points lie on a plane, it should in principle be possible to find all relative positions from just measuring distances. But can I actually do this without professional software and measurement equipment?

The equipment I had was a 50m measurement tape, and a 15m laser range finder. Both have their shortcomings. The measurement tape needs 2 people to operate. It is somewhat elastic, on a 10m distance, you can stretch it maybe one cm. When it is lying on the ground straight, but without being pulled forcibly, I found it to agree with the laser range finder to within 1cm. The laser, on the other hand, has a hard time with the shrubs and trees on the property. Even places that look perfectly clear turned out to be tricky because of high grass, tree branches, or uneven surfaces that were hard to pinpoint with the laser.

Having obtained a list of measurements, the next task is to compute the location of the points. Since we only measured distances between some adjacent points, our measurement can not tell us where on the globe these points are located. We need to arbitrarily fix the position of one of the points as a reference, I randomly labeled the points and chose point 1 to be at a fixed position. This still leaves the possibility to rotate the points around point 1. Coincidentally, our property has one fence that is visibly straight. This is very handy since it is a natural choice for one of the axes of our coordinate system. We can impose the condition that the fence is the $y$ -axis, that is, it has coordinate $x=0$ . Now several of the points are either in, or very close to, the fence. The $x$ -coordinates of all these points are known by explicitly measuring their distance to the fence, this reduces the number of coordinates that need to be computed. All remaining points can then in principle be fixed by measuring distances to multiple know points.

Having obtained a list of distances, how can we concretely compute the positions of the points? It took me a few attempts to find a method that is reasonably stable and reliable. The measured distances are a list of entries of the form $\{i,j,d_{ij} \}$ , where $i$ and $j$ are the indices of the two points, and $d_{ij}$ is the measured distance between them. Then, I define a function

$f = \left( (x_1-x_2)^2 + (y_1- y_2)^2 - d_{1,2}^2 \right)^2 + \left( (x_1- x_3)^2 + \ldots \right)^2.$

The goal is to choose the coordinates $x_i$ and $y_i$ of all points $i$ such that the function $f$ is minimized. Morally, the function $f$ measures how much the distances differ from the measured distances, hence ideally it should become zero, but due to measurement errors, it will never be exactly zero. The powers in $f$ are somewhat arbitrary, one can also choose them differently.

I have made a Mathematica notebook that implements this procedure and test it with random data. Obviously, the functionality of minimizing a polynomial function in many variables is also available in hundreds of other programs and libraries. Here is an example of 20 randomly generated points (black), and the starting points for the numerical minimization (red), which are also random:

The minimization procedure amounts to moving the red points towards the black ones. I fixed the coordinates of point 1, and the $x$ -coordinates of points 2,3,4,5, this leaves 34 coordinates free. Hence, we have a numerical minimization problem in 34 dimensions. In this particular run, I included a random subset of 124 out of the possible 180 measurements between the points, and 1% random error for each of them. Due to this error, one can not expect to find the true positions. The outcome of the minimization looks like this:

Most points are relatively close to the true position, but there are a few outliers. We can, for all 124 input measured distances, compute how much the computed distance between 2 points differs from the measured one. The relative differences are in the range of a few per cent, which is not much bigger than the measurement error.

A closer look shows that the central points are positioned relatively well, whereas the outer ones show a larger deviation. This is expected: The outer points are “measured from only one direction”, which is less accurate than having independent measurements into all directions. This can be nicely visualized in the following way: The function $f$ has unit m $^4$ , therefore $\sqrt[4]{f}$ can be measured in m or in cm. We now set a threshold, say $c=50$ cm, and we investigate point $i$ . Namely, we keep all other points fixed, and move point $i$ around. This will increase the value of $f$ (because we had previously determined the locations of all points such that $f$ is minimized). We mark in grey the area where $\sqrt[4]{f}$ increases by at most $c$ . To understand why this is sensible, assume that $i$ is incident to 10 measured distances. Then, moving $i$ by 1cm will change all of these distances: Those in the direction of motion will change by plus or minus 1cm, those perpendicular will not change much. Let’s say 4 were approximately in the direction of motion, then we have increased $\sqrt[4]{f}$ by 4cm. In the same situation, if there are 20 measurements, $\sqrt[4]{f}$ is already increased by 8cm. That is, the region where $\sqrt[4]{f}$ is less than a fixed $c$ will be smaller if the point is determined by more measurements. Additionally, this analysis is sensitive to the direction of the measurements: When all 10 measurements come approximately from the same direction, as is the case for our outlier points, then moving these points perpendicularly does barely increase $\sqrt[4]{f}$ . Hence, the uncertainty might be large in one direction, and small in another. This is in fact what we find:

In practice, there will always be errors in the measurements, or points that have been forgotten, etc. To remedy this, it is useful to include the measurements into the image. I draw a fine green line for a measurement that is consistent with the computed distance to within 3 cm, because more I can not hope for in realistic conditions, where often the location of a tree or fence is not even clear to within 5 cm: Exactly which corner of a pole do we want to measure?

When the computed difference is larger than the measured one by more than 3cm, I draw a red line, when it is shorter by more than 3cm, a blue line. I make these lines bold if the error is large. This way, one can print the image, and go and double-check precisely those measurements which are drawn bold. This plot also makes it particularly obvious why the two points to the left have large uncertainty in the vertical direction: There simply are no measurements to constrain them in that direction.

For the real-world measurements I did, there are around 135 points, and some 650 measurements between them. I was only interested in the points on our own property, but it turned out that accuracy can be greatly increased by including a few “outlier” points, such as street lamps or a tree on the neighbor’s backyard, in order to have measurements into all directions for the points I really care about. The colored plot of measurement deviations has been really helpful, often it revealed mistakes at previous measurements, e.g. misread a number, or measured two different corners of the same object at different times so that their location was not consistent.

Clearly, this method is slow compared to e.g. a lidar or aerial image. However, it is relatively robust for places with many obstacles and convoluted geometry. Moreover, it is nice that this global optimization takes into account all information to find the best possible solution. I have since also measured two neighboring properties, and this has implicitly improved the accuracy of the previous measurement because the border points and fence poles are now determined with more measurements.