Lecture 21 Notes: Relative Orientation, Binocular Stereo, Structure from Motion, Quadrics, Camera Calibration, Reprojection

6.801/6.866: Machine Vision, Lecture 21 1 Lecture 21: Relative Orientation, Binocular Stereo, Structure from Motion, Quadrics, Camera Calibration, Reprojection Today, we will continue discussing relative orientation, one of our main problems in photogrammetry. We concluded our previous lecture with a discussion of critical surfaces - that is, surfaces which make solving the relative rotation problem difficult. In general, these critical surfaces are quadric surfaces, which take the form:  2  2  2 x y z ± ± ± =1 a b c This form is simple: (i) Centered at the origin, (ii) Axes line up. In general, the form of these quadric surfaces may be more complicated, with both linear and constant terms. The signs of these qudratic terms determine the shape, as we saw in the previous lecture: 1. + + + (ellipsoid) 2. + + − (hyperboloid of one sheet) - this is a class of critical surfaces 3. + − − (hyperboloid of two sheets) 4. − − − (imaginary ellipsoid) For the solutions to the polynomial R2 + · · · − R = 0 (the class of hyperboloids of one sheets), the solutions tell us why these critical surfaces create solution ambiguity: 1. R = 0: I.e. the solution is on the surface - the origin of the righthand system is on the surface. Two interpretations from this; the first is for binocular stereo, and the second is for Structure from Motion (SfM): (a) Binocular Stereo: The righthand system lies on the surface. (b) Structure from Motion: In Structure from Motion, we “move ont the surface” in the next time step. 2. R = −b: This solution is the origin of the lefthand system (we move the solution left along the baseline from the righthand system to the lefthand system). Two interpretations from this; the first is for binocular stereo, and the second is for Structure from Motion (SfM): (a) Binocular Stereo: “The surface has to go through both eyes”. (b) Structure from Motion: We must start and end back on the surface. 3. R = kb: Because this system (described by the polynomial in R, then we have no scaling, giving us a solution. This suggests that the entire baseline in fact lies on the surface, which in turn suggests: • This setting/surface configuration is rare. • This surface is ruled - we can draw lines inscribed in the surface. This suggests that our critical surface is a hyperboloid of one sheet. It turns out that we have two rulings for the hyperboloid at one sheet (i.e. at every point on this surface, we can draw two non-parallel lines through the surface that cross at a point). 1 Hyperboloids of one sheet, as we saw in the previous lecture, are not the only types of critical surfaces. Intersections of planes are another key type of critical surface, namely because they are much more common in the world than hyperboloids of one sheet. Analytically, the intersection of planes is given by the product of planes: (a1 X + b1 Y + c1 Z + d1 )(a2 X + b2 Y + c2 Z + d2 ) = 0 This analytic equation describes the intersection of two planes - gives us a quadric surface. In order for this intersection of planes to not have any constants, one of the planes must be the epipolar plane, whose image is a line. The second plane is arbitrary and can take any form. Practically, we may not run into this problem in this exact analytic form (i.e. a surface that is an exact instance of a hyperboloid of one sheet), but even as we venture near this condition, the noise amplification factor increases. As it turns out, using a higher/wider Field of View (FOV) helps to mitigate this problem, since as FOV increases, it becomes less and less likely that the surface we are imaging is locally similar to a critical surface. One way that we can increase the FOV is through the use of spider heads, which consist of 8 cameras tied together into a rigid structure. Though it requires calibration between the cameras, we can “stitch” together the images from the 8 cameras into a single “mosaic image”. Figure 1: A spider heads apparatus, which is one technique for increasing the FOV, thereby reducing the probability of surfaces resembling critical surfaces. Next, we will switch our discussion to another of our photogrammetry problems - interior orientation. 1.1 Interior Orientation Recall that interior orientation was another one of the 4 photogrammetry problems we study, and involves going from 3D (world coordinate) to 2D (image coordinates). Specifically, we focus on camera calibration. Have we already touched on this? Somewhat - with vanishing points. Recall that the goal with vanishing points was to find (x0 , y0 , f ) using calibration objects. This method: • Is not very accurate, nor general across different domains/applications • Does not take into account radial distortion, which we need to take into account for high-quality imaging such as aerial photography. What is radial distortion? We touch more on this in the following section. 1.1.1 Radial Distortion This type of distortion leads to a discrepancy between what should be the projected location of an object in the image, and what it actually projects to. This distortion is radially-dependent. This type of distortion becomes more apparent when we image lines/edges. It manifests from having a center of distortion. The image of an object does not appear exactly where it should, and the error/discrepancy depends on the radius of the point in the image relative to the center of distortion. This radial dependence is typically approximated as a polynomial in the radius r: r = ||r||2 δx = x(k1 r2 + k2 r4 + · · · ) δy = y(k1 r2 + k2 r4 + · · · ) 2 Figure 2: Radial distortion is in the direction of the radius emanating from the center of projection to the point in the image plane where a point projects to in the absence of radial distortion. Note that the error vector δr is parallel to the vector r - i.e. the distortion is in the direction of the radial vector. Many high-quality lenses have a defined radial distortion function that gives the distortion as a function of r, e.g. the figure below: Radial Distortion Function 1.5 e(r) 1 0.5 0 0 1 2 3 4 5 r These functions, as we stated above, are typically approximated as polynomial functions of r. How do we measure radial distortion? One way to do so is through the use of plumb lines. These lines are suspended from the ceiling via weights, which allows us to assume that they are straight and all parallel to one another in the 3D world. Then, we take an image of these lines, with the image centered on the center lines. We can estimate/assess the degree of radial distortion by identifying the curvature, if any, of the lines as we move further away from the center lines of the image. Figure 3: Plumb lines are a useful way of estimating radial distortion, because they allow us to estimate the degree of curvature between lines that should be straight and parallel. Plumb lines allow us to estimate two types of radial distortion: (i) Barrel Distortion and (ii) Pincushion Distortion, depicted below. The radius of curvature from these lines allows us to measure k1 (the first polynomial coefficient for r2 ): 3 Figure 4: Barrel and pincushion distortion, which we can absorb from plumb lines. A more subtle question: is it better to go from )a) undistorted → distorted coordinates, or from (b) distorted → undistorted coordinates? It is often desirable to take approach (b) (distorted → undistorted) because we can measure the distorted coordinates. Wse can use series inversion to relate the distorted and undistorted coodinates. This affects the final coordinate system that we do optimization in, i.e. we optimize error in the image plane. 1.1.2 Tangential Distortion and Other Distortion Factors There are other types of distortion that we can analyze and take into account: 1. Tangential Distortion: This is a problem we have solved, but to contrast this with radial distortion, rather than the distortion acting radially, the distortion acts tangentially (see figure below). With tangential distortion, our offsets take the form: δx = −y(1 r2 + 2 r4 + · · · ) δy = x(1 r2 + 2 r4 + · · · ) Figure 5: Tangential distortion acts tangentially relative to the radius - in this case, in a counterclockwise fashion (as we rotate θ). 2. Decentering: If the center of distortion is not the principal point (center of the image with perspective projection), then we get an offset that depends on position. This offset is typically small, but we still need to take it into account for high-quality imaging such as aerial photography. 3. Tilt of Image Plane: If the image plane is tilted (see the figure below), then magnification will not be uniform and the focus will not be perfect. In practice, to fix this, rather than changing the tilt of the image plane, you can instead insert a compensating element to offset this tilt. 4 Figure 6: Tilt of the image plane is another factor that can affect distortion of images, as well as variable magnification/scaling across the image. Oftentimes, it can be corrected via a compensating element that offsets this tilt, rather than realigning the image plane. Taking these effects into account allows for creating more sophisticated distortion models, but sometimes at the expense of overfitting and loss of generalizability (when we have too many degrees of freedom, and our estimates for them become so precariously configured that our model is not applicable for a slight change of domain). These distortion models can also be used for nonlinear optimization protocols such as Levenberg-Marquadt (LM) or GaussNewton (GN) (see lecture 20 for a review of these). 1.2 Tsai’s Calibration Method This calibration method relies on using a calibration object for which we have precise knowledge of the points on our 3D points of the object. In camera calibration, correspondences are defined between points in the image and points on the calibration object. Here, we encounter the same reason that vanishing point-based methods are difficult - it is hard to directly relate the camera to calibration objects. When we take exterior orientation into account, we generally perform much better and get higher-performing results. Exterior orientation (2D → 3D) seeks to find the calibration object in space given a camera image. Adding exterior orientation also increases the complexity of this problem: • Interior orientation has 3 Degrees of Freedom • Exterior orientation has 6 Degrees of Freedom (translation and rotation) Combined, therefore our problem now possesses 9 DOF. We can start solving this problem with perspective projection and interior orientation. 1.2.1 Interior Orientation Beginning with perspective projection, we have: xI − xO Xc = f Zc xI − x O Xc = f Zc Where: 1. (Xc , Yc , Zc )T are the camera coordinates (world coordinate units) 2. (xI , yI ) denote image position (row, column/grey level units) 3. (xO , yO , f ) denotie interior orientation/principal point (column/grey level and pixel units) Can we modify the measurements so we are able to take radial distortion into account? We need a good initial guess, because, numerical, iterative approaches that are used to solve this problem precisely exhibit 5 multiple mimima, and having a good initial guess/estimate to serve as our “prior” will help us to find the correct minima. From here, we can add exterior orientation. 1.2.2 Exterior Orientation As we are solving for translation and rotation, our equation for exterior orientation equation becomes (written as a matrix-vector problem): ⎡ ⎤ ⎡ ⎤⎡ ⎤ ⎡ ⎤ r11 r12 r13 Xs tx Xc ⎣ Yc ⎦ = ⎣r21 r22 r23 ⎦ ⎣ Ys ⎦ + ⎣ty ⎦ Zc r31 r32 r33 Zs tz Xc = RXs + t The figure below illustrates the desired transformation between these two reference frames: Figure 7: Exterior orientation seeks to find the transformation between a camera and a calibration object. Where the matrix R and the translation t are our unknowns. From this, we can combine these equations with our equations for interior orientation to find our solutions. 1.2.3 Combining Interior and Exterior Orientation Combining our equations from interior and exterior orientation: • Interior Orientation: Xc xI − xO = f Zc xI − xO Xc = f Zc • Exterior Orientation: ⎡ ⎤ ⎡ Xc r11 ⎣ Yc ⎦ = ⎣r21 Zc r31 r12 r22 r32 ⎤⎡ ⎤ ⎡ ⎤ r13 Xs tx r23 ⎦ ⎣ Ys ⎦ + ⎣ty ⎦ r33 Zs tz Xc = RXs + t Then we can combine these equations into: Xc r11 Xs + r12 Ys + r13 Zs + tx xI − xO = = r31 Xs + r32 Ys + r33 Zs + tz f Zc yI − yO Yc r21 Xs + r22 Ys + r23 Zs + ty = = r31 Xs + r32 Ys + r33 Zs + tz f Zc Written this way, we can map directly from image coordinates to calibration objects. How can we find f and tz , as well as mitigate radial distortion? 6 Let’s start by looking at this problem in polar coordinates. Radial distortion only changes radial lengths (r), and not angle (θ). Changing f or Zc changes only radius/magnification. Let’s use this factor to “forget” about the angle. We can also remove additional DOF by dividing our combined equations by one another: xI − xO = yI − y0 = r11 Xs +r12 Ys +r13 Zs +tx r31 Xs +r32 Ys +r33 Zs +tz r21 Xs +r22 Ys +r23 Zs +ty r31 Xs +r32 Ys +r33 Zs +tz r11 Xs + r12 Ys + r13 Zs + tx r21 Xs + r22 Ys + r23 Zs + ty We can see that combining these constraints removes our DOF tz , giving us 8 constraints rather than 9. Cross-multiplying and term consolidation gives: (XS yI0 )r11 + (YS yI0 )r12 + (ZS yI0 )r13 + yI0 tx − (XS xI0 )r21 − (YS xI0 )r22 − (ZS xI0 )r23 − xI0 ty = 0 A few notes on this consolidated homogeneous equation: • Since our unknowns are the rij values, we have a linear combination/equation of our unknowns! The coefficients they are multiplied by are either our observed image positions or the positions of the correspondences/keypoints on the calibration object. • Note that we subtract out an estimate of the center of projection. The estimation error of this center of projection affects adjacent points the most, as a slight miscalculation can greatly affect the calculations above - and as a result we throw out all correspondences that are -close (for some  > 0 to the center of the image. • Each correspondence between image coordinates and the points on the calibration object gives one of these equations: “This point in the image corresponds to another point on the calibration object”. As we have asked for the other frameworks and systems we have analyzed: How many correspondences do we need? Let’s consider which unknowns appear in the equations, and which do not: • DOF/unknowns in the equations: r11 , r12 , r13 , r21 , r22 , r23 , tx , ty (8 unknowns) • DOF/unknowns not in the equations: r31 , r32 , r33 , tz , f (5 unknowns) Can this inform of how we can solve this problem? Let’s consider the following points: • We are not yet enforcing orthonormality of the rotation matrices, i.e. RT R = I. • We can solve for the last row of the rotation matrix containing the unknowns r31 , r32 , r33 through finding the cross product of the first two rows of this matrix. • Having 8 unknowns suggests that we need at least 8 correspondences, but because this equation is homogeneous, there is scale factor ambiguity (i.e. doubling all the variables does not change the righthand side value of the equation). The solution to account for this scale factor ambiguity is to set/fix one of these variables to a constant value. Fixing an unknown means we now have 7 DOF (rather than 8), and therefore, we actually need 7 correspondences. Solving with 7 (or more) correspondences gives us a multiple of the solution. Solving with exactly 7 correspondences gives us 7 equations and 7 unknowns, which we can solve with, for instance, Gauss-Jordan elimination. More correspondences are desirable because they allow us to estimate both solutions and the error associated with these solutions. 0 0 0 0 0 0 • These solutions give us values for the 8 unknowns in the equation above: r11 , r12 , r13 , r21 , r22 , r23 , tx , ty = 1 (where ty is the fixed parameter to account for scale ambiguity). However, as aforementioned these values are scaled versions of the true solution, and we can compute this scale factor by noting that because the rows of our (soon to be) orthonormal rotation matrix must have unit length - therefore the scale factor is given by computing the length of the rows of this matrix: s= p 0 2 r11 1 0 2 + r0 2 + r12 13 0 2 r21 1 0 2 + r0 2 + r22 23 OR s= p A good sanity check: these two variables should evaluate to approximately the same value - if they do not, this is often one indicator of correspondence mismatch. 7 With this scale factor computed we can find the true values of our solution by multiplying the scaled solutions with this scale factor: 0 0 0 0 0 0 s(r11 , r12 , r13 , r21 , r22 , r23 , tx , ty = 1) → r11 , r12 , r13 , r21 , r22 , r23 , tx , ty Note that we also have not (yet) enforced orthonormality of the first two rows in the rotation matrix, i.e. our goal is to have (or minimize): 0 0 0 0 0 0 r11 r21 + r12 r22 + r13 r23 =0 To accomplish this, we will use the idea of “squaring up”. 1.2.4 “Squaring Up” Squaring up is one method for ensuring that vectors are numerically/approximately orthogonal to one another. This technique answers the question of “Given a set of orthogonal vectors, what is the nearest set of orthogonal vectors?” The figure below highlights that given vectors a and b, we can find the nearest set of rotated vectors a0 and b0 using the following, which states that adjustments to each vector in the pair are made in the direction of the other vector: a0 = a + kb b0 = b + ka Since we want this dot product to be zero, we can solve for the scale k by taking the dot product of a0 and b0 and setting it equal to zero: a0 · b0 = (a + kb)(b + ka) = 0 = a · b + (a · a + b · b)k + k 2 + a · b = 0 This equation produces a quadratic equation in k, but unfortunately, as we approach our solution point, this sets the second and third terms, which k depends on, to values near zero, creating numerical instability and consequently making this a poor approach for solving for k. A better approach is to use the approximation:   a·b k≈− 2 It is much easier to solve for this and iterate over it a couple of times. I.e. instead of using the approach above, we use: x= −b ± 2c √ b2 − 4ac Where we note that the “standard” quadratic solution is given by: √ −b ± b2 − 4ac x= 2a It is important to consider both of these approaches because of floating-point precision accuracy - in any case, as we approach the solution point, one of the two forms of these solutions will not perform well, and the other form will perform substantially better. With this, we are able to iteratively solve for k until our two vectors have a dot product approximately equal to zero. Once this dot product is approximately zero, we are done! We have enforced orthonormality of the first two rows of the rotation matrix, which allows us to find an orthonormal third row of the rotation matrix by taking a cross product of the first two rows. Therefore, we then have an orthonormal rotation matrix. 1.2.5 Planar Target Planar target is an example calibration object that possesses desirable properties, namely, as we will see below, that it reduces the number of correspondences we need from 7 to 5. Additionally, planar targets are generally easy to make, store, and allow for high camera calibration accuracy. They are typically mounted on the side of wheels so that they can be rotated. They are typically given analytically by the plane: Zs = 0, i.e. the XY-plane where Z is equal to 0. 8 Figure 8: A geometric view of a planar target, which is a calibration target whose shape is geometrically described by a plane. As a result of this, we no longer need to determine the rotation matrix coefficients r13 , r23 , and r33 (i.e. the third column of the rotation matrix, which determines how Zs affects Xc , Yc , and Zc in the rotated coordinate system. Mathematically: ⎡ ⎤ ⎡ ⎤⎡ ⎤ ⎡ ⎤ Xc r11 r12 r13 Xs tx ⎣ Yc ⎦ = ⎣r21 r22 r23 ⎦ ⎣ Ys ⎦ + ⎣ty ⎦ Zc r31 r32 r33 0 tz As a result of this, our equations become: Xc r11 Xs + r12 Ys + tx xI − xO = = r31 Xs + r32 Ys + tz f Zc yI − yO Yc r21 Xs + r22 Ys + ty = = r31 Xs + r32 Ys + tz f Zc And, as we saw before, dividing these equations by one another yields a simplified form relative to our last division for the general case above: xI − x O = yI − y0 = r11 Xs +r12 Ys +tx r31 Xs +r32 Ys +tz r21 Xs +r22 Ys +ty r31 Xs +r32 Ys +tz r11 Xs + r12 Ys + tx r21 Xs + r22 Ys + ty As we did before for the general case, cross-multiplying and rearranging gives: (XS yI0 )r11 + (YS yI0 )r12 + yI0 tx − (XS xI0 )r21 − (YS xI0 )r22 − x0I ty = 0 Now, rather than having 8 unknowns, we have 6: r11 , r12 , r21 , r22 , tx , ty . Because this is also a homogeneous equation, we again can account for scale factor ambiguity by setting one of DOF/parameters to a fixed value. This in turn reduces the DOF from 6 to 5, which means we only need 5 correspondences to solve for now (compared to 7 in the general case). As before, if we have more than 5 correspondences, this is still desirable as it will reduce estimation error and prevent overfitting. One potential issue: if the parameter we fix, e.g. ty , evaluates to 0 in our solution, this can produce large and unstable value estimates of other parameters (since we have to scale the parameters according to the fixed value). If the fixed parameter is close to 0, then we should fix a different parameter. 1.2.6 Aspect Ratio This is no longer a common issue in machine vision now, but when aspect ratios were set differently in analog with photo-diode arrays for stepping in the x and y directions, this required setting an additional parameter to analytically describe the scale of x relative to y, and this parameter cannot be solved without a planar target. 1.2.7 Solving for tz and f Now that we have solved for our other parameters, we can also solve for tz and f : 9 • These unknowns do not appear in other equations, but we can use our estimates from our other parameters to solve for tz and f . • We can again use the same equations that combine interior and exterior orientation: (r11 Xs + r12 Ys + r13 Zs + tx )f − x0I tz = (r31 Xs + r32 Ys + r33 Zs )xI0 (r21 Xs + r22 Ys + r23 Zs + ty )f − yI0 tz = (r31 Xs + r32 Ys + r33 Zs )yI0 Since we only have two unknowns (since we have found estimates for our other unknowns at this point), this now becomes a much simpler problem to solve. With just 1 correspondence, we get both of the 2 equations above, which is sufficient to solve for our 2 unknowns that remain. However, as we have seen with other cases, using more correspondences still remains desirable, and can be solved via least squares to increase robustness and prevent overfitting. • One problem with this approach: We need to have variation in depth. Recall that perspective projection has multiplication by f and division by Z; therefore, doubling f and Z results in no change, which means we can only determine tz and f as a ratio tfz , rather than separately. To remove this consequence of scale ambiguity, we need variation in depth, i.e. the calibration object, such as a plane, cannot be orthogonal to the optical axis. Figure 9: In order to obtain f and tz separately, and not just up to a ratio, it is crucial that the calibration object does not lie completely orthogonal to the optical axis - else, these parameters exhibit high noise sensitivity and we can determine them separately. It turns out this problem comes up in wheel alignment - if machine vision is used for wheel alignment, the calibration object will be placed at an angle of 45 or 60 degrees relative to the optical axis. 1.2.8 Wrapping it Up: Solving for Principal Point and Radial Distortion To complete our joint interior/exterior orientation camera calibration problem, we need to solve for principal point and radial distortion. As it turns out, there is unfortunately no closed-form solution for finding the principal point. Rather, we minimize image error using nonlinear optimization (e.g. Levenberg-Marquadt or Gauss-Newton). This error is the 2D error in the (x, y) image plane that describes the error between predicted and observed points in the image plane: Figure 10: Image error that we minimize to solve for our principal point and radial distortion. We can formulate the total error from a set of observations as the squared errors from each correspondence. Note that xI refers to the observed point, and xP refers to the predicted point after applying our parameter estimate transformations to a point on the calibration object corresponding to an observed point. Ideally, this error is zero after we have found an optimal set of parameters, i.e. for each i: xIi − xPi = 0 yIi − yPi = 0 10 Where do xPi and yPi come from? They come from applying: • Rotation matrix estimate (R) • Translation vector estimate (t) • Principal Point Estimate (x0 , y0 ) • Radial Distortion Polynomial (k1 , k2 , · · · ) Using these equations, we can take a sum of squared error terms and solve for our parameters using nonlinear optimization with Levenberg-Marquadt (LM). This can be implemented using LMdif in the MINPACK package (Fortran), as well as other packages in other languages. Recall that LM performs nonlinear optimization with L2 regularization. One small problem with this approach: LM only works for unconstrained optimization, which is not the case when we only use orthonormal rotation matrices. To overcome this, Tsai’s Calibration Method used Euler Angles, and we can use either Gibb’s vector or unit quaternions. Comparing/contrasting each of these 3 rotation representations to solve this problem: • Euler Angles: This representation is non-redundant, but “blows up” if we rotate through 180 degrees. • Gibb’s Vector: This representation is also non-redundant, but exhibits a singularity at θ = π. • Unit Quaternions: This representation is redundant but exhibits no singularities, and this redundancy can be mitigated by adding an additional equation to our system: o o q·q−1=0 Finally, although LM optimization finds local, and not global extrema, we have already developed an initial estimate of our solution through the application of the combined interior/exterior orientation equations above, which ensures that extrema that LM optimization finds are indeed not just locally optimal, but globally optimal. 1.2.9 Noise Sensitivity/Noise Gain of Approach As we have asked for other systems, it is important to consider how robust this approach is to noise and perturbations. A few notes/concluding thoughts: • As we saw above, this approach does not produce good estimates for our parameters tz and f if the calibration object is orthogonal to the optical axis. In addition to not being able to separate the estimates for these parameters, both of these parameters (because they are defined by a ratio) have a very high/infinite noise gain. • Since these methods are largely based off of numerical, rather than analytical methods, we can best evaluate noise gain using Monte Carlo methods, in which we use noise that is well-defined by its statistical properties, e.g. zero-mean Gaussian noise with finite variance denoted by σ 2 . We can analyze the estimated noise gain of the system by looking at the parameter distribution in the parameter space after applying noise to our observed values and making solution predictions. • Higher-order coefficients of radial distortion are highly sensitive to noise/exhibit high noise gain. 11