Multiple views

1Rays that should meet

One image fixes the ray on which a point of space lies, but not where on that ray. The note on epipolar geometry used a second image to find where the point can appear in it. This note goes two steps further. Given the pixels of a point in two or more images, it finds the point itself. And given its pixels in two images, it predicts where the point appears in a third.

A point M seen in image i as the pixel m_i lies on the ray from the optical center C_i through m_i. It lies on the ray of every image that sees it, so it is where those rays meet. Finding a point from the rays of its images is called triangulation, and finding it this way is 3D reconstruction. Here that word means locating a point. It is not what computed tomography reconstructs, which is what an object is made of at each of its points.

Real pixels are measured, and the cameras are estimated, so the rays almost never meet exactly, not even two of them. The question becomes which point of space best explains the pixels.

The example

The example on this page is made for it, with round numbers where possible. An X-ray source, the optical center, can stand anywhere on a horizontal circle of radius 1000 mm around a vertical axis Z, and each image it takes is 2000 × 1600 pixels. Image 1 is taken from one point of that circle, image 2 from a quarter turn around it, and image 3 from 175° around, almost across the object from image 1. The object is a flat bar 45 mm long, pointed at both ends, which stands 12° from vertical with its lower tip level with the source. Its two tips are clicked in each of the three images, and, as clicks by hand are, each click is a few pixels off. Since the bar is known, what the images give can be checked against where its tips really are. Lengths are in millimeters and image coordinates in pixels.

2Two equations per view

For one image, write its camera P with entries p_jk, the pixel as m = [x y 1]^T and the point as M = [X Y Z 1]^T. The projection λm = PM is three equations:

The unknowns are X, Y and Z, and λ too, a new one for each image. The third equation gives λ outright: it is the depth of M in front of the camera, in millimeters in this example, since the first three entries of the third row of P form a unit vector. Putting it into the first two removes it:

That is why each image gives two equations: it brings three equations and one new unknown, λ, and one equation is spent on it. The two left are linear in X, Y and Z, because the pixel (x, y) is a known number.

Each of them has a picture. The first holds for every point whose image has this x, whatever its y: the points that project onto one column of the image. They form a plane through the optical center. The second is the plane through the center that projects onto the row y. The ray of m is where the two planes meet. So one image gives a line, not a point, and a second image is needed to fix a point on it.

3Least squares from n views

Each of n images gives its two rows, with its own camera and pixel. Stacked, they form one system Q[X Y Z]^T = r, where Q is 2n × 3 and r has 2n entries, the two rows of image i in places 2i − 1 and 2i. As in the note on camera calibration, least squares gives the point that makes the sum of the squared entries of Q[X Y Z]^T − r smallest:

The hat marks an estimate. Call that point M. It exists when Q has rank 3, that is, when the 2n planes do not all share one line.

Two views and the epipolar constraint

Two images give four equations for three unknowns, one more than needed, so the four have a common solution only for some pairs of pixels. Write them as [Q | −r][X Y Z 1]^T = 0: a 4 × 4 system with right side 0, which has a solution other than 0 only when det [Q | −r] = 0. Each row of [Q | −r] is xp₃ − p₁ or yp₃ − p₂ for one of the two pixels, with p_j the rows of its camera. So the determinant is linear in m₁ and linear in m₂, which is the form of m₂^TFm₁, and it vanishes for every pair of corresponding points. It is the epipolar constraint, times a factor that depends only on the cameras. For images 1 and 2, with F = [BC₁]_×BA⁺ as in the note on epipolar geometry, the factor is −6.40·10⁷, whatever the pixels. It is the determinant det G of that note, with λ₁ and λ₂ already eliminated.

So two rays meet exactly when their pixels satisfy the epipolar constraint. The clicks of the example do not. In image 2, they lie 10.83 and 6.86 pixels from the epipolar lines of their partners in image 1, so their rays miss each other. Least squares splits the difference: with images 1 and 2 only, the upper tip's M projects 5.55 pixels from its click in image 1 and 5.41 from the one in image 2.

Views that face each other

The rays of two images must cross at a clear angle. The sources of images 1 and 3 are 175° apart, on opposite sides of the bar, so the two rays through a point of it are almost one line: they meet at about 172°. Along that line the two images barely see a difference, and a few pixels of error move M far along it. With only these two images, the lower tip comes out 10.99 mm from where it is and the upper tip 3.76 mm, and the bar 48.21 mm long. Images 1 and 2 give 45.34 mm, and images 2 and 3 give 44.43 mm, with both tips within about a millimeter of where they are. In the equations, the columns of Q are nearly dependent: its singular values are 11319, 11292 and 786.1, against 11325, 8014 and 7987 for images 1 and 2.

More views

With three images, Q is 6 × 3. The errors of the clicks partly cancel, and no single pair of images decides the result. Images 1, 2 and 3 together give 45.17 mm, with the tips 0.51 and 0.56 mm from where they are.

4What the fit minimizes

Least squares makes QM − r small, writing M for [X Y Z]^T, but what is that residual? Take the first row of image i. Its left side minus its right side is (p₃ · M)x_i − p₁ · M, writing p_j · M for the product of row j of P_i with [X Y Z 1]^T. The first product is the depth λ_i of M, and the second is λ_ix_i, where m_i = (x_i, y_i) is where M projects into image i. So each row of the residual is λ_i(x_i − x_i) or λ_i(y_i − y_i), a depth in millimeters times a pixel error, and least squares minimizes

That is the error in pixels weighted by the squared depth, in mm·px: an image whose source is farther from the point weighs more. The quantity with a plain meaning is the reprojection error, the distance in pixels from each clicked m_i to the projection m_i of the point found, with λ_im_i = P_iM. The best point in that sense minimizes

No formula gives that point, because x_i and y_i divide by the depth, which depends on M. It is found by iterating, with a gradient method for example, starting from M. When all the depths are alike, the two answers nearly agree. Here they run from 985 to 1020 mm: minimizing J moves the tips by 0.026 and 0.024 mm, and changes the length by 0.028 mm.

5A third view from two

Now three images, with cameras A, B and C. Given corresponding points m₁ and m₂ in the first two, where does their point of space appear in the third? Epipolar geometry already answers it. m₃ lies on the epipolar line of m₁ in image 3, and on that of m₂:

where F₁₃ is the fundamental matrix of images 1 and 3, from A and C, and F₂₃ that of images 2 and 3, from B and C. So m₃ is where the two lines cross, and the point where two lines cross is their cross product, as in the note on homogeneous coordinates:

It is also a test. m₁, m₂ and m₃ correspond exactly when m₃ is that point: the condition is necessary and sufficient.

When the two lines are one

The crossing is not defined when ℓ₁₃ and ℓ₂₃ are the same line. ℓ₁₃ is where image 3 cuts the epipolar plane of M, C₁ and C₃, and ℓ₂₃ is where it cuts the plane of M, C₂ and C₃. Both planes contain M and C₃, and they are the same plane when it holds C₁ and C₂ as well. That happens in two ways:

• The three optical centers are on one line, as when one camera takes the three images while moving in a straight line. This is common, and then it happens for every point.
• The centers are not on one line, but M lies on the plane through the three of them.

Near either case the lines cross at a shallow angle, and a small error in m₁ or m₂ moves their crossing far. And even when the lines cross well, the method takes two steps, two lines and then their crossing, rather than one formula.

In the example

Every source here is on one circle, so the plane through any three of them is the plane of the circle, Z = 0. It passes through every optical center, so each camera sees it edge-on, as a line, and it is the same line in every image: the row v = 800, level with each camera's axis.

The bar's lower tip lies on that plane. Predicted in image 3 from its clicks in images 1 and 2, the two epipolar lines cross at only 0.46°, and their crossing, (298.48, 797.28), lands 842.56 pixels from where the tip was clicked in image 3, (1141, 805). Right on that row, one pixel of error in m₂ moves the crossing by more than a hundred pixels, and from the tip's true pixels the two lines are both the row v = 800, with no single crossing at all. The plane is also why images 1 and 3 alone put this tip 10.99 mm off: for a point on it, the rows of both images give the same plane, Z = 0, and only their columns, nearly parallel, place it along the rays. The upper tip is 44 mm above the plane. There the lines cross at 37.93°, and their crossing lands 25.30 pixels from its click.

6The trifocal tensor

A way around the epipolar lines starts from the three projections themselves. Pick a 4 × 4 matrix H whose first three rows are A, and change coordinates to M̃ = HM. Since A = [I | 0]H,

and λ₁m₁ = ÃM̃, λ₂m₂ = B̃M̃, λ₃m₃ = C̃M̃: the same projections, in coordinates where the first camera is as simple as it can be. H must be invertible. Adding the row [0 0 0 1] under A works whenever C₁ is a finite point, because then the left 3 × 3 block of A is invertible.

The ray of m₁

λ₁m₁ = [I | 0]M̃ says that the first three components of M̃ are λ₁m₁. With the fourth equal to 1,

As λ₁ runs over all depths, M̃ runs along the ray of m₁. Finding the point means finding λ₁.

The depth from one coordinate

Project this M̃ into image 2 with B̃. Write b̃_j for the first three entries of row j of B̃, as a row, and b̃_j4 for its last one. The first coordinate of the pixel is x₂ = (λ₁b̃₁m₁ + b̃₁₄) / (λ₁b̃₃m₁ + b̃₃₄), and solved for λ₁:

This is the depth at which the ray of m₁ projects onto the column x₂ of image 2: where it crosses the plane through C₂ that images onto that column, as in section 2. The same step in image 3, with the rows c̃_k of C̃ and the column x₃, gives λ₁ = (c̃₁₄ − c̃₃₄x₃) / (x₃c̃₃m₁ − c̃₁m₁).

Four trilinearities

If m₁, m₂ and m₃ correspond, the two depths are the same λ₁. Setting them equal, multiplying out and collecting the terms in 1, x₂, x₃ and x₂x₃ gives an equation that is linear in m₁, in x₂ and in x₃:

Each T^jk here is a column of three numbers, T^jk = [T₁^jk T₂^jk T₃^jk]^T, with

for i, j, k = 1, 2, 3. Using y₂ instead of x₂, or y₃ instead of x₃, gives three more equations of the same kind. The four are the trilinearities:

The 27 numbers T_i^jk form the trifocal tensor of the three images. Like F, they depend only on the cameras, not on the points, so they are computed once for all the points. And like F, they are homogeneous: another fourth row for H multiplies all 27 by one number, which changes none of the four equations. The four equations are independent, and they hold only for corresponding points: when all four hold, m₁, m₂ and m₃ correspond.

The same equations from a determinant

As with F in the note on epipolar geometry, the same equations come from the projections written as one system, Gv = 0. Here G is 9 × 7, with the rows of Ã, B̃ and C̃ in its first four columns and m₁, m₂ and m₃ in its last three, and v = [M̃; −λ₁; −λ₂; −λ₃]. Corresponding points give a solution v ≠ 0, so every 7 × 7 matrix made of 7 of the 9 rows of G has determinant 0. Leaving out two rows of one camera keeps only the other two cameras, and gives their epipolar constraint, such as m₃^TF₂₃m₂ = 0. Leaving out one row from each of two cameras keeps all three, and those determinants give the trilinearities.

7Transfer with the tensor

Knowing m₁ and m₂, the first two trilinearities are linear in x₃ and in y₃. Solved for them, they give the point in image 3 directly:

Why it works: x₂ fixes λ₁, which is where the ray of m₁ crosses the plane through C₂ and the column x₂, and that point of space, projected by C, is m₃. It is a reconstruction from one image and a half, followed by a projection into the third. No epipolar lines are crossed, so the plane of the three centers does no harm. For the bar's lower tip, the transfer gives (1152.43, 803.88), 11.49 pixels from the click in image 3, where the crossing of the epipolar lines was 842.56 pixels off. What is left is the clicks' own error: in images 1 and 2, they are 6.86 pixels from each other's epipolar line.

Its own weak spot

If the epipolar line of m₁ in image 2 is the column x₂ itself, the plane through C₂ and that column holds the whole ray of m₁, and x₂ cannot fix λ₁. The last two trilinearities, which use y₂, then do the job. In the example, the epipolar lines are close to rows, since the epipoles lie on the row v = 800, mostly far to the sides of the images, so the form with x₂ is the one to use. The form with y₂ fails on that row, just as the crossing of epipolar lines does. For images whose sources face each other, such as 1 and 3, the epipoles fall inside the images, at (1349.29, 800) in image 1. Taken as the first two images, the form with x₂ then fails for pixels m₁ close to the column through that epipole, whose epipolar lines are columns too.

Is it a third match?

m₁, m₂ and m₃ correspond if m₁ and m₂ satisfy the epipolar constraint, and the three satisfy the first two trilinearities, or the last two. The first two use x₂ but never y₂, and the epipolar constraint between m₁ and m₂ is what checks y₂. With measured points, both tests get a tolerance:

with F₁₂m₁ = [a b c]^T, as in the note on epipolar geometry. The epipolar test comes first because the transfer never looks at y₂. Paired the wrong way, the upper tip's click in image 1 with the lower tip's in image 2, the clicks are 350.93 pixels off each other's epipolar line. But the bar stands upright in image 2, so its two clicks there have almost the same x₂, and the transfer lands at (1068.56, 447.49), 12.86 pixels from the upper tip's click in image 3, as if the three matched.

8Four views and more

With a fourth image, the same tests chain. m₁, m₂, m₃ and m₄ correspond if m₁, m₂, m₃ correspond and m₂, m₃, m₄ correspond: the shared m₂ and m₃ fix one point of space, so all four are its images. That makes three checks: the epipolar constraint for m₁ and m₂, two trilinearities for images 1, 2 and 3, and two for images 2, 3 and 4, or in practice the distance from m₄ to the point transferred from m₂ and m₃. Four views have tensors of their own, 81 numbers with 16 equations, but the chain makes them unnecessary. Chains like this follow a detected object through a sequence of X-ray images, where each new image must agree with the two before it.

9Summary

The key ideas next to their NumPy equivalents. In code, the cameras are mat_a, mat_b and mat_c, or p for any one of them, the pixels are homogeneous vectors such as m1 = np.array([x1, y1, 1]), and mat_q and mat_r are Q and r.