Triangulation methods and uncertainty

This is the basic midpoint method: it computes the point in space that minimizes the distance between the two observation rays. This is the simplest method, but also produces the most bias. Not recommended. Implemented in mrcal.triangulate_geometric() (in Python) and mrcal_triangulate_geometric() (in C).

Described in the "Triangulation Made Easy" paper above. The method is a close approximation to a reprojection error minimization (the "right" approach above) if we have pinhole lenses. Implemented in mrcal.triangulate_lindstrom() (in Python) and mrcal_triangulate_lindstrom() (in C).

Described in the "Closed-Form Optimal Two-View Triangulation Based on Angular Errors" paper above. Minimizes the L1 norm of the observation angle error. Implemented in mrcal.triangulate_leecivera_l1() (in Python) and mrcal_triangulate_leecivera_l1() (in C).

Described in the "Closed-Form Optimal Two-View Triangulation Based on Angular Errors" paper above. Minimizes the L-infinity norm of the observation angle error. Implemented in mrcal.triangulate_leecivera_linf() (in Python) and mrcal_triangulate_leecivera_linf() (in C).

Described in the "Triangulation: Why Optimize?" paper above: this is the "Mid2" method. Doesn't explicitly minimize anything, but rather is a heuristic that works well in practice. Implemented in mrcal.triangulate_leecivera_mid2() (in Python) and mrcal_triangulate_leecivera_mid2() (in C).

Described in the "Triangulation: Why Optimize?" paper above: this is the "wMid2" method. Doesn't explicitly minimize anything, but rather is a heuristic that works well in practice. Similar to leecivera_mid2, but contains a bit of extra logic to improve the behavior for points very close to the cameras (not satisfying \(r \gg b\)). Implemented in mrcal.triangulate_leecivera_wmid2() (in Python) and mrcal_triangulate_leecivera_wmid2() (in C).

We want to capture the effect of two different sources of error:

• Calibration-time noise. We propagate the noise in chessboard observations obtained during the chessboard dance. This is the noise that we propagate when evaluating projection uncertainty. This is specified in the --q-calibration-stdev argument to mrcal-triangulate or in the q_calibration_stdev argument to mrcal.triangulate(). This is usually known from the calibration, and we can request the calibrated value by passing a stdev of -1. See the relevant interface documentation (just-mentioned links) for details.
• Observation-time noise. Each triangulation processes observations \(\vec q\) of a feature in space. These are noisy, and we propagate the noise. As with the calibration-time noise, this noise is assumed to be normally distributed, independent in \(x\) and \(y\). This is specified in the --q-observation-stdev argument to mrcal-triangulate or in the q_observation_stdev argument to mrcal.triangulate(). A common source of these pixel observations is a pixel correlation operation where a patch in one image is matched against the second image. Corresponding pixel observations observed this way are correlated: the noise in \(\vec q_0\) not independent of the noise in \(\vec q_1\). I do not yet know how to estimate this correlation, but the tools are able to ingest and propagate such an estimate: using the --q-observation-stdev-correlation commandline option to mrcal-triangulate or the q_observation_stdev_correlation argument to mrcal.triangulate().

A big point to note here is that repeated observations of the same feature have independent observation-time noise. So these observation-time errors average out with multiple observations. This is not true of the calibration-time noise however. Using the same calibration to observe a feature multiple times will produce correlated triangulation results. So calibration-time noise is biased, and it is thus essential to make and use low-uncertainty calibrations to minimize this effect.

The test-triangulation-uncertainty.py test generates synthetic models and triangulation scenarios. It can be used to produce an illustrative diagram:

test/test-triangulation-uncertainty.py  \
  --do-sample                           \
  --cache write                         \
  --observed-point -2 0 10              \
  --fixed cam0                          \
  --Nsamples 200                        \
  --Ncameras 2                          \
  --q-observation-stdev-correlation 0.5 \
  --q-calibration-stdev 0.2             \
  --q-observation-stdev 0.2             \
  --make-documentation-plots ''

Here we have two cameras arranged in the usual left/right stereo configuration, looking at two points somewhere ahead. We generate calibration and observation noise, and display the results in the horizontal plane. The vertical dimension is insignificant here, so it is not shown, even though all the computations are performed in full 3D. For each of the two observed points we display:

• The empirical noise samples, and the 1-sigma ellipse they represent
• The predicted 1-sigma ellipse for the calibration-time noise
• The predicted 1-sigma ellipse for the observation-time noise
• The predicted 1-sigma ellipse for the joint noise

We can see that the observed and predicted covariances line up nicely. We can also see that the observation-time noise acts primarily in the forward/backward direction, while the calibration-time noise has a much larger lateral effect. This pattern varies greatly depending on the lenses and the calibration and the geometry. As we get further out, the uncertainty in the forward/backward direction dominates for both noise sources, as expected.

In the above plot, the uncertainties are displayed in the coordinate system of the left camera. But, as described on the projection uncertainty page, the origin and orientation of each camera's coordinate system is subject to calibration noise:

So what we usually want to do is to consider the covariance of the triangulation in the coordinates of the camera housing, not the camera coordinate system. We achieve this with "stabilization", computed exactly as described on the projection uncertainty page. We can recompute the triangulation uncertainty in the previous example (same geometry, lens, etc), but with stabilization enabled:

test/test-triangulation-uncertainty.py  \
  --do-sample                           \
  --cache write                         \
  --observed-point -2 0 10              \
  --fixed cam0                          \
  --Nsamples 200                        \
  --Ncameras 2                          \
  --q-observation-stdev-correlation 0.5 \
  --q-calibration-stdev 0.2             \
  --q-observation-stdev 0.2             \
  --stabilize                           \
  --make-documentation-plots ''

We can now clearly see that the forward/backward uncertainty was a real effect, but the lateral uncertainty was largely due to the moving camera coordinate system.

As mentioned above, the calibration-time noise produces correlations (and thus biases) in the triangulated measurements. Since the test-triangulation-uncertainty.py command triangulates two different points, we can directly observe these correlations. Let's look at the magnitude of each element of \(\Var {\vec p_{01}}\) where \(\vec p_{01}\) is a 6-dimensional vector that contains both the triangulated 3D points: \(\vec p_{01} \equiv \left[ \begin{array}{cc} \vec p_0 \\ \vec p_1 \end{array} \right]\). If we had only observation-time noise, \(\vec p_0\) and \(\vec p_1\) would be independent, and the off-diagonal terms in the covariance matrix would be 0. However, we also have calibration-time noise, so the errors are correlated:

As before, the exact pattern varies greatly depending on the lenses and the calibration and the geometry, but calibration-time noise always creates these correlations. To reduce these correlations and the biases they cause: lower the uncertainty of your calibrations by dancing better

When propagating noise, mrcal makes the very common assumption that everything is locally linear. This makes things simple, and is right most of the time. However, when running the triangulation routines with near-parallel rays, this assumptions can break down.

Let's run another simulation, but observing a more distant point, with more observation-time noise, no calibration-time noise, and gathering more samples:

test/test-triangulation-uncertainty.py  \
  --do-sample                           \
  --cache write                         \
  --observed-point -200 0 2000          \
  --fixed cam0                          \
  --Nsamples 2000                       \
  --Ncameras 2                          \
  --q-observation-stdev-correlation 0.5 \
  --q-observation-stdev 0.4             \
  --stabilize                           \
  --make-documentation-plots ''

The range to the observed point:

The two points in the synthetic world are at \((\pm 200, 0, 2000)m\) so the true range is ~ \(2010m\). We see that the calibration-time noise has little effect here. More importantly, we also see that the predicted distribution of the range to the point is gaussian (as we assume), but the empirical distribution is not gaussian: there's a much more significant tail on the long end. This makes sense. If the observation rays are near-parallel, small errors that make the rays more parallel push the range to infinity; while small errors that bring the rays together have a more modest, finite effect.

Similarly, when we look at the distance between our two points we get this distribution:

We see the same asymmetric non-gaussian distribution. Empirically I observe this distance-between-points distribution become more non-gaussian, faster than the range-to-point distribution.

At this time I do not know how much this matters or what to do about it, but these limitations are good to keep in mind.