Part I: Camera Calibration using 3D calibration object (50 points)
We want to calibrate the camera of a robot vehicle. We will use a linear method as described in the lectures.
We place a large cubic frame of size 4 meters on the road several meters in front of the vehicle.
of the eight corners of the cubic frame are dened with respect to a world coordinate system with its axes
parallel to the cube edges and with its origin at the center of the cube.
The world coordinates of the cube
2 2 2,
-2 2 2,
-2 2 -2,
2 2 -2,
2 -2 2,
-2 -2 2,
-2 -2 -2,
2 -2 -2
We detect the corresponding cube corners at the following pixel positions in the camera image:
422 323; 178 323; 118 483; 482 483; 438 73; 162 73; 78 117; 522 117
1. Draw the image points, using small circles for each image point.
2. Write a Matlab function that takes as argument the homogeneous coordinates of one cube corner and
the homogeneous coordinates of its image, and returns 2 rows of the matrix P (slide 30 of the Camera
Calibration pdf document). This matrix P will be used to compute the 12 elements of the projection
matrix M such that λpi = MPi
3. Use this Matlab function to generate 2 rows of the matrix P for each cube corner and its image and
obtain a matrix with 16 rows and 12 columns. Print matrix P.
4. Now we need to solve the system Pm = 0. Find the singular value decomposition of matrix P using
matlab svd function. The last column vector of V obtained by svd(P) should be the 12 elements in row
order of the projection matrix that transformed the cube corner coordinates into their images. Print
the matrix M.
5. Now we need to recover the translation vector which is a null vector of M. Find the singular value
decomposition of matrix M = UΣV
. The 4 elements of the last column of V are the homogeneous
coordinates of the position of the camera center of projection in the frame of reference of the cube
(as in slide 36). Print the corresponding 3 Euclidean coordinates of the camera center in the frame of
reference of the cube.
6. Consider the 3×3 matrix M’ composed of the rst 3 columns of matrix M. Rescale the elements of this
matrix so that its element m33 becomes equal to 1. Print matrix M’ . Now let the rotation matrices
be as dened in slide 38 where the axes e1, e2, e3 are the x, y, z axes respectively. Therefore M’ can
be written as M0 = KRT
7. We will perform the RQ factorization of M’ in several steps. First, nd a rotation matrix Rx that sets
the term at position (3,2) to zero when Rx is multiplied to M’. The cosine and sine used in this matrix
are of the form
cos(θx) = m33/sqrt(m2
33 + m2
sin(θx) = −m32/sqrt(m2
33 + m2
Note that the term at position (3, 2) would also be set to zero if the signs of cos(θx) and sin(θx) were
reversed, but this would lead to nding a negative focal length for the camera. So we should choose the
signs that leads to a positive focal length. Compute the angle θx of this rotation in degrees. Compute
matrix N = M’ ∗ Rx. Print Rx, θx and N.
8. The element n31 of N is small enough so that there is no need for a rotation Ry. However, element n21
is large and a rotation matrix Rz is needed to set it to zero. Compute the rotation matrix Rz using
cosine and sine of the form
cos(θz) = n22/sqrt(n
21 + n
sin(θz) = −n21/sqrt(n
21 + n
Compute the rotation angle θz in degrees. This angle is actually very small.
9. Since we factorized out Rz we can directly compute the calibration matrix K, how? Compute K and
rescale so that its element K33 is set to 1. Print K. What are the focal lengths of the camera in pixels?
What are the pixel coordinates of the image center of the camera?
Part II: Camera Calibration using 2D calibration object (70 points)
In this part we are going to implement camera calibration from multiple images of 2D planes. Additionally
we will learn how to augment images with virtual objects. We will follow the method described in the
book chapter on camera calibration by Zhengyou Zhang which was proposed in his paper Flexible Camera
Calibration by Viewing a Plane from Unknown Orientations – Zhang, ICCV99. Start by carefully reading
Section 2.4 of that Chapter.
Calibration Grid and images
We know the following about the grid. The grid is 9 squares in width and 7 in height. Each square is 30mm
x 30mm. If we select the bottom left corner of the grid to be the origin of the world coordinate system, and
the grid to be the plane corresponding to Z=0 , then we know the 3D coordinates of each corner in that
For calibration we will use four images: images2.png, images9.png, images12.png, images20.png
Corner Extraction and Homography computation (10 points)
First we want to extract the four corners of the calibration grid form each image. We will use the grid
corners to estimate the homographies relating two images. We can manually get the four grid corners from
each image. One way to let a user manually select points in matlab is using ginput function.
Once the 4 corners are extracted, compute the homography H that relates the grid 3d coordinates to
the corners. Use the function homography2d which is provided. Repeat this for the four images provided.
Report the computed H for each image [deliverable].
Computing the Intrinsic and Extrinsic parameters (30 points)
Now given the four homographies, follow the instructions in section 2.4.4 to compute the intrinsic parameters
and extrinsic parameters. We need to linearize the two constraints in Eq 2.19 and 2.20 into two equations in
a homogeneous system as in Eq 2.25.
Then solve for b and estimate the intrinsic parameters as described in
page 21. Print the computed matrix B and the intrinsic parameters. Then compute and print R, t for each
Verify that your rotation matrix is in fact a rotation matrix, print RT R, is it an identity as it should be?
We can enforce R to be a rotation matrix by SVD decomposition of R and setting the singular values to
ones, i.e., set the rotation matrix to UV T where R = UΣV
Print the new R and RT R after enforcing the rotation matrix constraints. [deliverable]
Improving accuracy (30 points)
Since we used four manually entered points to compute the homography. A small error in one of the points
will directly eect the computed homography. To x this we are going to estimate the homography from all
• First given the computed homographies from Section 2, compute the approximate location of each
grid corner in the image. (Hint : This can be done since we know the 3d locations of the grid corners
and the approximate homography. Call these points p_approx. Create a gure with the image and
approximate grid locations. Call this Figure 1 : Projected grid corners [deliverable]
• Second, using the provided Harris function detect Harris corners in the image and display them. Use
the following parameter values for the Harris detection : sigma = 2, thresh = 500, radius = 2.
[cim, r, c, rsubp, csubp] = harris(rgb2gray(im), sigma, thresh, radius, disp);
Here r is the y-coordinate of the Harris corner, c is the x-coordinate of the Harris corner, rsubp is
the y coordinate with subpixel accuracy, csubp is the x coordinate with subpixel accuracy. Use rsubp,
csubp. Create a gure with image and overlayed Harris corners. Call this Figure 2 : Harris corners .
• Third, compute the closest Harris corner to each approximate grid corner. (You may nd it useful to
use the dist2.m function provided). Let these closest Harris corners be p_correct. Create a gure with
the image and p_correct overlayed. Call this Figure 3 : grid points . [deliverable]
• Finally, compute a new homography from p_correct, print H [deliverable]
• Repeat this for the other three images. Then use the homographies to estimate K and R, t for each
image. Report your K, R’s, and t’s [deliverable]. Save your results, you will need to use them in
• Using the new computed H, compute the errors between points in p_correct and points you get by
projecting grid corners to the image (Hint there is no need to use R, t for projecting) . Call this
err_reprojection. Report your result. [deliverable]
• Now repeat the process using 4 images. Compare your results to your previous results and those of
part 2 [deliverable].
• Can you suggest a way this can be done automatically (i.e without rst letting the user manually select
the 4 corners) ?
Part III: Augmented Reality 101 (50 points)
Augmenting an Image (30 points)
Now we would like to use our computed homographies from part II to map a clip art image onto the grid
such that it seems to be part of the grid. The image should be synthesized such that the clip art bottom left
corners is the same as the grids (0,0) corner.
When tting the clip art you should rescale it to t the grid
while keeping the clip art aspect ratio. Using your computed homography nd a way to map your image
on the grid such that you image will have the same projective distortion as the grid. If the clip art have
any white pixels you should make these pixels appears transparent when overlayed over the grid.
image of the four images in Part II, create a gure with the original image and your virtual clip art overlayed
over the grid.
Your image should be one of the images provided in the clipart directory. To nd which clip
art you are supposed to map, take the last 4 digits of your ruid id add them up and use the clip art le
corresponding to the rst digit in the sum, i.e., if the rst 4 digits in your RUID are 5243, use the clip art
Augmenting an Object (20 points)
Now we would like to augment our images with 3D objects. For our purposes we are going to use a cube as
a virtual object. We will only render the cube as a wire frame and we would like its base to be locate on
the 3×3 grid of squares in the bottom left corner of our grid.
The cube should be standing up from the grid.
First print the 3D coordinates of the cube. Then, nd a way to use your computed H, K, R, t to synthesize
new images with the virtual cube inserted.
Extra credit (20 points)
Do one of the following :
1. Instead of augmenting a cube, augment a general mesh from a 3D le of your choice.
2. Can you nd a way to estimate the intrinsic and extrinsic parameters from only two images of the
grid. What assumptions on the intrinsic parameters are needed to achieve this. (Hint the answer can
be found in Sec 2.4)