cvFindExtrinsicCameraParams2() to track the rigid transform that represents the pose of a chessboard calibration object relative to a camera. You will also use the theory of the pinhole camera model to draw simple 3D graphics which depict the chessboard and camera from an arbitrary external viewpoint.- Get a baseline program running which displays live images from your camera, as you did in HW2. You may use one of the example programs that were given with HW2. Again, aim for a resolution of
and about 10FPS. Get your printout of the chessboard calibration object and make arrangements so that you can move either it or the camera while keeping the chessboard flat and fully visible in the image. Your program should have a feature that pauses/unpauses the live video by hitting the spacebar, as the example programs cvdemoandcvundistortfrom HW2 do. - Undistort the live image using the calibration data you have already measured for your camera. If you’ll be coding in C++, you may use the
cvundistortexample program; otherwise, use that as a reference to develop your own equivalent. Use these undistorted images for the remainder of the problem. - Define a right-handed 3D coordinate system
that moves with the chessboard and which represents object frame. The origin of is the lower left interior corner, the 
axis points to the right across the columns, and the 
axis points up across the rows. 
thus points out of the chessboard. Distance in is measured in units of the chessboard pitch, so each of the 
interior chessboard corners is defined by coordinates 
. For example, the coordinates of the lower left corner in are 
and those of the upper right corner are 
. - Define a second right-handed 3D coordinate system
that moves with the camera and which represents camera frame. The origin of is the center of projection of the camera in the pinhole model, points towards the right side of the image, points down in the image, and points straight forward along the optical axis. Distance in is also measured in units of the chessboard pitch. - For every captured frame, use the OpenCV functions
cvFindChessboardCorners()andcvFindExtrinsicCameraParams2()to calculate a
3D rigid transformation matrix 
that takes points in to points in . Only calculate for frames where cvFindChessboardCorners()returns nonzero, which indicates there is a high likelihood that the whole chessboard was detected correctly in the image. For frames where you calculate
, draw a live reconstruction of the chessboard in magenta (
) by calling cvLine()to draw
lines in a grid covering the interior chessboard corners. (Draw nothing for frames where you did not calculate .)There are a few ways to interpret this. In particular, you could just use the pixel results of
cvFindChessboardCorners()to get the image pixel locations of the endpoints of the lines, or you could calculate the 3D coordinates of those points in frame, transform to , and then backproject into the image using the camera matrix. Because of what we are going to do next, please do the latter (*). Specifically, for each line (i) calculate the 3D endpoints 
in , (ii) transform from to by calculating 
and similar for 
, (iii) backproject into the image by calculating 
and similar for 
, where 
is the 
camera intrinsic matrix
,and then perform the perspective division
to get the pixel 
where 
appears, and similar for 
.(*) The extent to which these two different approaches should give the same result depends precisely on how well the camera calibration models your actual camera, and on how accurate
cvFindChessboardCorners()was.On top of the drawing you made above, draw three colored line segments, all in 3D in the frame
Your drawings both of the chessboard and of these axes should be updated for every new video frame. If you have implemented things correctly, it should look like they are virtually “stuck to” the chessboard as it moves around in 3D.: one from to 
in red, another from to 
in green, and a third from to 
in blue. These represent the axis triad of itself. Again, do this only for frames where you calculated .Now introduce one more
rigid transform 
into the projection chain. We will use this as a navigation transform to render a representation of both the chessboard and the camera from a simulated external viewpoint. This viewpoint will be modeled as if we had a second camera, with the same calibrated pinhole model as your actual camera, looking at the whole scene (chessboard and camera) from “outside”. Let 
be the camera frame of this simulated camera. Then is defined as the transform that takes points from to , so the composite transform 
takes points from to .To make things interesting, you will handle keypresses that can change
, effectively allowing the simulated camera to “fly around”. We know that since is a 3D rigid body transform (i.e. a transform in 
) that it has six DoF, and that we can divide these up as three translations and three rotations. We will take the strategy of using six keys to control positive/negative translations along all three axes of , and four more keys to control positive/negative rotations about the and axes of (we don’t allow the virtual camera to rotate around the 
axis of , this is a common design in part because this kind of “rolling” camera motion can be visually confusing).Set
to the identity transform initially. Implement keyboard handlers (i.e. using cvWaitKey(), or if you are using one of theCvBaseimplementations from HW2, just overridehandleKeyExt()) for theiandkkeys that pitch the virtual camera by a rotation of
or 
about the axis of , respectively, and similarly make j,l(that’s the letter j and the letter l keys) yaw by
rotations about the axis of . Make e,dshift the virtual camera
along the axis of , s,fshift by along the axis of , and a,zshift by along the axis of . Finally, make rreset to the identity transform. See below for an explanation of why we use these keys.Add drawing code for (i) an RGB axes triad for frame
(this is in addition to the one you already implemented for frame ), and (ii) a very simple representation of the actual camera as a unit-altitude rectangular pyramid in yellow (
), also in , representing the horizontal and vertical FoV. I.e. the apex of this pyramid is 
in and the four base points are given by the four combinations
where
are the dimensions of your images in pixels (e.g. ).Whenever
is not the identity matrix (i.e. whenever the virtual viewpoint does not correspond to the actual camera view), hide the actual captured images and make all your drawings on a white background. Conversely, whenever is identity, show the captured images instead of a white background, but hide your axis triad and pyramid in frame . As before, for all frames where you computed the transform from to (i.e. whenever the chessboard was successfully detected, irrespective of ), draw the axis triad in frame and the reconstructed chessboard grid.Be sure that you can still change
while capture is paused.
cvCalcOpticalFlowPyrLK() to track feature points, and you will use cvKalmanPredict() and cvKalmanCorrect() to incorporate a simple model of point motion.Starting with a basic program that displays images from your camera (aim for
, 10FPS as usual), such as cvdemofrom HW2, add code that (i) adds and removes feature points where the user clicks the mouse on the image, displaying them as green dots, (ii) usescvCalcOpticalFlowPyrLK()to track the points from frame to frame, and (iii) removes existing feature points either when they are clicked (at their current position, wherever it may be) or when their tracking fails or degrades below a threshold that you determine is reasonable.Clearly, the provided
Your program should have a feature that pauses/unpauses the live video by hitting the spacebar, as the example programs likelkdemo.cwill be useful. Note that we will not need the ability to automatically initialize the feature points (cvGoodFeaturesToTrack()) or the “night mode” oflkdemo.cvdemofrom HW2 do. Make sure that adding and removing track points by clicking with the mouse works even when capture is paused.Define a
measurement covariance matrix
cvCalcOpticalFlowPyrLK(). Though it is possible to estimate
from the content of the image (see the end of the section on the Wikipedia entry for Harris corners), for our purposes we will take a simpler approach. We will set 
and 
and 
to constants. These are, respectively, the horizontal and vertical standard deviations of the error in the tracking results from cvCalcOpticalFlowPyrLK(); they are measured in pixels. Initialize
, and make the e,d,s,fkeys increment/decrement (s-,f+) and (e+,d-) by 0.1. (See below for an explanation of why we use these keys.) Clamp the minimum value of each to 0.1 and display the new settings for both on the console (i.e. with printf()orcoutorSystem.out.println()) whenever either is changed.- Add code that draws error ellipses in green for all of the points returned by
cvCalcOpticalFlowPyrLK(), which we will now call the measurements. To do this you will need to define a confidence limit
, which means that the actual location of each feature should be inside the corresponding ellipse with probability 
. Clearly, larger values of must correspond to larger ellipses. Handle 
as a special case: simply don’t draw any ellipses. Initialize 
and make the a,zkeys increment/decrement by 
, respectively. Clamp to the range 
, and display its new value on the console whenever it is changed. (Be sure to also re-draw all the error ellipses whenever changes, even when capture is paused.) Starting with the third frame received after a feature is added, it is possible to predict (or at least to make a reasonable guess about) the new location of the feature
given its two previous locations 
and 
(here we are using 
to index the sequence of video frames). There are various ways to do this depending on the assumptions we make about object motion. We will simply model the point to move with locally (locally in the temporal domain, that is) constant velocity in the image. We will also assume that the amount of time that passes between adjacent captured frames is constant. Thus, we initialize an estimate of the velocity of the feature from its first two observed locations as 
and we initialize an estimate of its position as 
. The state of the feature is now the 4D vector 
). (Note that we only explicitly calculate the state by these equations to initialize it; later state updates will be completely determined by the Kalman filter.)Use
cvKalmanPredict()to calculate
for every feature point that has been defined for at least 3 frames (because our model does not include any cross-correlation between features, you may use a separate CvKalmanstructure for each feature). Model the uncertainty in this prediction with a diagonal covariance matrix
of the form
Initialize
For every feature that has been present for at least 3 frames, draw the predicted location
(10 times larger than the initial setting of and because we are significantly less certain of our simple model than of the measurements from cvCalcOpticalFlowPyrLK()), and make thei,k,j,lkeys increment/decrement
(j-,l+) and 
(i+,k-) by 0.1. (See below for an explanation of why we use these keys.) Clamp the minimum value of each to 0.1 and display the new settings for both on the console whenever either is changed. and its error ellipse, which is given by the upper-left submatrix of and the current value of the confidence limit (use the same value of for all parts of this problem). Draw these in red.For every feature that has been present for at least 3 frames, use
cvKalmanCorrect()to combine the predicted state
and the observation 
(the output of cvCalcOpticalFlowPyrLK()for that feature) into an estimated state
and estimate covariance matrix 
. Display the estimated location 
and its error ellipse (given by the upper left submatrix of and the current value of ) in blue, and use the predicted values to set the input prevFeaturesfor the next call tocvCalcOpticalFlowPyrLK().
for a 2D Gaussian distribution
)
of the random variable modeled by the distribution will be found in the error ellipse 
(this notation simply means the ellipse 
is some function of 
, 
and 
coordinate axes), and 
, but if we did, that would effectively correspond to an ellipse of infinite size.
). The details are given 
:
since 
then this expands out to
and 
of the ellipse: for the horizontal case, plug in 
. The scale factors, as we have defined them, are thus the “stretchings” of the radius of the original circle to make the ellipse radii. What happens if one of the scale factors, say 
, is zero? Then the corresponding radius is zero, but 
and 
for an arbitrary CCW rotation 
. We actually multiply by 
to have the effect of rotating back to standard position, where the stretching is done:
since 
(naming the two columns of 
then 
is a unit vector pointing along the major axis, which has radius 
is a unit vector pointing along the minor axis, which has radius 
(and if 
we have a circle, which is a special case of an ellipse).
:
:
—it has to be decomposable into the form 
, where 
(symmetric) and 
for any nonzero 
, 
(i.e. square) matrix 
to some other vector 
, by simple left multiplication: 
. In general, the “output” vector 
will not point in the same direction as 
for some scalar 
It is not too hard to convert this statement into a prescription for finding such scaling values 
, which are called the eigenvalues of 
is the 
is not invertible, which is the same as saying its 
, and the roots of the polynomial will give the various eigenvalues 
are a basis for the space of solutions to the specific equation 
. This is a homogeneous system (right hand side 
) so it will, in general, have not just a single solution but an entire linear space of solutions.
is an eigenvector then so is 
for any 
; so we may normalize all the eigenvectors (this step is not strictly necessary, but assume we do so). If all the eigenvectors are linearly independent then it can be shown there are precisely 
for 
. If we re-index the corresponding eigenvalues as 
(allowing for multiple copies of eigenvalues with multiplicity greater than 1), then we can factor 
is an 
and 
.
. Assume we’ve already verified that both are non-negative. Now we can easily find 
and 
. NOTE: as we move on to covariance matrices below, will will find ourselves starting with an eigendecomposition of the inverse of 
calculated from 
. Thus 
.
exists and that 
. The probability distribution associated with a covariance matrix that is only positive semi-definite is not strictly a Gaussian, however, it is still possible to draw an error ellipsoid for it. Our method for drawing ellipsoids allows this case; the ellipsoid will have one or more zero radii in the direction(s) corresponding to its zero eigenvalues.)
(sets of 
, and 
increases.
of the exponent for which any desired fraction 
of the probability distribution is given by the 
is some as-yet unspecified function of 
). 
.
is the eigendecomposition of 
is easy to compute since the inverse of a diagonal matrix is trivially formed by taking the reciprocal of each entry on the diagonal.) This also shows that the inverse of a symmetric positive definite matrix is itself symmetric and positive definite. (If the matrix is only positive semi-definite, then replace the inverse with the pseudoinverse 
, which differs only in that 
is formed by taking the reciprocal of each entry on the diagonal of 
taken to be zero.)
and 
(here 
, so the ellipse radii are directly proportional to their square roots; note that a zero covariance matrix eigenvalue will correspond to a zero radius, and an infinite covariance matrix eigenvalue will correspond to an infinite radius), and the eigenvectors of 
(this can be derived from the development in Appendix A of [1]). But for other dimensions it is not even possible in general to write down 
just has 
and an 
covariance matrix