CS4277-CS5477-Assignment 1 Solved

Introduction
In this assignment, you will get to implement a panorama stitching algorithm. As discussed in the lectures, images taken using a purely-rotating camera are related by a homography, which can be estimated by at least four point correspondences. You will first implement a homography estimation algorithm, then add a RANSAC scheme to reject outliers. You will also implement a image warping function and use it to stitch together multiple images into a single coherent panorama.


Optional references: * HZ Book, Sections 4.1, 4.7, 4.8

Part 1: Homography Estimation from Point Correspondences
In this part, you will implement the Normalized Direct Linear Transformation (DLT) algorithm for homography estimation. This is a linear algorithm for determining the homography transformation H given at least four 2D to 2D point correspondences xi ↔ xi0 (in homogeneous coordinates). The goal is to compute the 3 × 3 homography matrix H such that Hxi = xi0 for each i.

1(a) Transformation using provided homography matrix
To ensure you understand the homography transformation, first implement the function that transform 2D points using a provided homography matrix. Implement the following function(s): transform_homography()

 


You should get the following if the function is implemented correctly:

Points after (translation)

[[1. 1.5 2. 2. 2. 1.5 1. 1. ]

[0.5 0.5 0.5 1. 1.5 1.5 1.5 1. ]]

Points after (rigid)

                  [[0.5                        0.933013 1.366025 1.116025 0.866025 0.433013 0.                          0.25           ]

                  [0.1                 0.35              0.6                       1.033013 1.466025 1.216025 0.966025 0.533013]]

Points after (affine)

[[0. 0.5 1. 1.25 1.5 1. 0.5 0.25]

                     [0. 0.                           0. 0.75 1.5 1.5 1.5 0.75]]

Points after (homography)

                  [[0.                    0.5                1.                   0.88              0.8                   0.466667 0.133333 0.08              ]

                 [0.                    0.                   0.                   0.6                1.                   1.                   1.                   0.6              ]]

1(b) Compute Homography Matrix
Now, implement the function to compute the homography matrix given pairs of correspondences using the normalized DLT algorithm.

Implement the following function(s): compute_homography()
•     Prohibited Functions:           cv2.findHomography(), cv2.getPerspectiveTransform(), np.linalg.solve(), np.linalg.lstsq()

•     You may use the following functions: np.linalg.svd()

Remember to use normalization, which will increase the robustness of the estimation.

 [3]: # Test case with exactly 4 points src1 = np.array([[0.0, 0.0], [331.0, 0.0], [331.0, 474.0], [0.0, 474.0]]) dst1 = np.array([[182.0, 94.0], [482.0, 48.0], [598.0, 466.0], [146.0, 533.0]]) homo1 = compute_homography(src1, dst1) # Test case with more than 4 points

src2 = np.array([[434.375, 93.625], [429.625, 190.625],

[533.625, 301.875], [452.375, 460.625],

[558.375, 188.625], [342.444, 362.596],

[345.625, 41.875], [341.625, 146.125]])

dst2 = np.array([[204.780, 92.367], [201.875, 190.625],

[297.125, 296.875], [224.446, 456.556],

[318.407, 192.155], [107.625, 371.375],

[109.875, 26.624], [106.625, 138.125]])

homo2 = compute_homography(src2, dst2)

 print('homo1:\n', homo1, '\n') print('homo2:\n', homo2)

You should obtain approximately the following (or some scale of it) if your function is implemented correctly:

homo1:

[[ 9.012217e-01 -1.792522e-01 1.820000e+02] [-1.394830e-01 5.490343e-01 9.400000e+01]

[-1.062811e-05 -7.075535e-04 1.000000e+00]]

homo2:

[[ 1.738142e+00 9.494425e-03 -4.466937e+02] [ 2.745643e-01 1.514698e+00 -1.213769e+02]

[ 1.160079e-03 2.069049e-05 1.000000e+00]]

 

Part 2: Image Warping using Homography
In this part, you will implement the image warping function. Given two images Isrc, Idst and the homography matrix H which relates the coordinates in the two images:

xdst = Hxsrc

, where xsrc and xdst are pixel coordinates of Isrc and Idst respectively, the objective is to warp Isrc onto Idst.

The obvious way of doing this is forward-interpolation, where for every pixel in Isrc, you compute where it should be in Idst and copy the pixel value over. However this will result in holes in the image (why?).

The better way is to do it in the reverse direction (i.e. backward interpolation), where we compute the corresponding coordinate in Isrc for every pixel in Idst. First create a HdstWdst × 2 matrix containing the coordinates of all the pixels in Idst, i.e,

                                                                                                                Wdst − 1        Wdst − 1T

,

                                                                                 0    1     2    . . .       Hdst − 2        Hdst − 1

where Wdst and Hdst are the width and height of Idst. Then use your transform_homography() function to compute the corresponding coordinates in Isrc. Lastly, for every pixel in the destination image, fill its value with the corresponding pixel in the source image; you can use interpolation methods to handle non-integer soure pixel coordinates.

Implement the following function(s): warp_image()
•     Prohibited Functions: cv2.warpPerspective()

•     You may use the following functions: cv2.remap(), np.meshgrid()

If you use cv2.remap, you might find the borderMode value of cv2.BORDER_TRANSPARENT useful. Also consider using bilinear interpolation in cv2.remap.

Let us now use our functions so far to replace the book cover with another image. If all above functions are implemented correctly, you should see the book cover being replaced by another image when you run the following code.

[4]: template = load_image('data/hzbook_2.jpg') original = load_image('data/hzbook_1.jpg')

 # homo is the homography matrix that maps points in template to original homo = np.array([[ 9.01221661e-01, -1.79252179e-01, 1.82000000e+02], [-1.39482959e-01, 5.49034320e-01, 9.40000000e+01],

[-1.06281109e-05, -7.07553504e-04, 1.00000000e+00]])

overlaid = warp_image(template, original, homo)

plt.figure() plt.imshow(template) plt.title('Template') plt.figure(figsize=(12,6)) plt.subplot(1,2,1) plt.imshow(original)

 

Part 3: Robust Homography Estimation using RANSAC
In this part, we upgrade our homography estimation algorithm to be robust even in the presence of outliers. Consider the following image stitching example containing 10 provided correspondences.

 [5]: im1 = load_image('data/pano_2.jpg') im2 = load_image('data/pano_3.jpg')

im1_points = np.array([[434.375, 93.625], [429.625, 190.625],

[533.625, 301.875], [452.375, 460.625],

[558.375, 188.625], [342.444, 362.596],

 

Notice that 2 of the 10 correspondences are wrong. Using just the correct correspondences yields the correct result:

[6]: im1_points_inliers = im1_points[:-2, :] # exclude the last 2 wrong␣

,→correspondences im2_points_inliers = im2_points[:-2, :] homo = compute_homography(im1_points_inliers, im2_points_inliers)

 stitched = warp_images_all([im1, im2], [homo, np.eye(3)]) plt.figure(figsize=(12, 8)) plt.imshow(stitched) print('Computed homography matrix:\n', homo)

Computed homography matrix:

[[ 1.738142e+00 9.494425e-03 -4.466937e+02] [ 2.745643e-01 1.514698e+00 -1.213769e+02]

[ 1.160079e-03 2.069049e-05 1.000000e+00]]

 

However, if we include the two outliers, notice that the alignment fails.

[7]: homo = compute_homography(im1_points, im2_points)

 stitched = warp_images_all([im2, im1], [np.eye(3), homo]) plt.figure(figsize=(12, 8)) plt.imshow(stitched);

 

The de-facto algorithm to solve this is RANdom SAmple Consensus (RANSAC). Your task is now to implement the robust homography computation which can handle outliers. RANSAC requires a error distance function to compute which correspondences are outliers. First, implement the symmetric transfer error distance measure as described in the lectures:

d  .

Implement the following function(s): compute_homography_error()
•     Prohibited Functions: cv2.findHomography()

•     You may use the following functions: np.linalg.inv(), and transform_homography() from the previous sections.

 [8]: src = np.array([[0.0, 0.0], [1.0, 0.0], [2.0, 1.0]]) dst = np.array([[1.0, 2.0], [1.5, 2.5], [2.5, 3]])

H = np.array([[0.5, 0.0, 1.0], [0.0, 0.5, 2.0],

[0.0, 0.0, 1.0]])

print('Error:', compute_homography_error(src, dst, H))

If implemented correctly, the previous code should print:

Error: [0. 1.25 2.5 ]

Now implement the RANSAC procedure covered in the lectures, using the error distance measure you just implemented:

Be sure to re-estimate the homography transformation using all the inliers at the end of the RANSAC procedure.

Implement the following function(s): compute_homography_ransac()
•     Prohibited Functions: cv2.findHomography()

•     You may use the following functions: compute_homography() function from the previous sections.

 

If you implemented the above correctly, the above images should be aligned correctly. The RANSAC algorithm detects outliers (drawn in red) and excludes them from the homography computation. The computed homography matrix should be the same as the previous sections:

[[ 1.738142e+00 9.494425e-03 -4.466937e+02] [ 2.745643e-01 1.514698e+00 -1.213769e+02]

[ 1.160079e-03 2.069049e-05 1.000000e+00]]

At this point, we now can perform automatic alignment of the image pair. Recall in the lectures that we can detect keypoints and match them between the two images. The matches generally will contain outliers but your robust homography computation can handle those and still output the correct homography.

Let’s see how our algorithm performs in such a scenario.

First, let us use OpenCV’s feature detection, extraction and matching pipeline to establish potential correspondences in the two images.

 

We can see that some of the correspondences are wrong. So how does our compute_homography_ransac() function perform?

 

 

If implemented correctly, the images should be aligned properly, with the inliers and outliers correctly detected.

 

Part 4: Adding more images
You are almost done! Now the last part is to extend the stitching to multiple images. To do so, suppose we have N images I1, . . . , IN. We perform pairwise matching between consecutive images (i.e. I1 → I2, I2 → I3, . . .). We can then set one of the images (e.g. middle one) as the reference image, and warp all the homographies to it.

Implement the function that computes the homography warping matrices for each image from pairwise homography matrices. Specifically, given H21, H32, . . . , HNN−1 which represents the homography matrices between consecutive images:

xk+1 = Hkk+1xk,

compute the absolute homography matrices H1ref , H2ref , . . . , HrefN−1 that can be used to warp all images to the reference image:

Implement the following function(s): concatenate_homographies()
[11]: # Initiate ORB detector. Detect keypoints and extracts descriptors for all␣

,→images.

orb = cv2.ORB_create() num_images = 4 all_im, all_kp, all_des = [], [], [] for i in range(num_images):

im = load_image('data/pano_{}.jpg'.format(i))

 kp, des = orb.detectAndCompute(cv2.cvtColor(im, cv2.COLOR_BGR2GRAY), None) all_im.append(im) all_kp.append(kp) all_des.append(des)

# Matches the descriptors between consecutive images via brute force matching,␣

,→and

# computes the homography matrices.

bf = cv2.BFMatcher(cv2.NORM_HAMMING, crossCheck=True) pairwise_homographies = [] for i in range(num_images-1):

matches = bf.match(all_des[i], all_des[i+1])

im1_points_orb, im2_points_orb = matches2pairs(matches, all_kp[i],␣

,→all_kp[i+1]) h_i1_i, mask = compute_homography_ransac(im1_points_orb, im2_points_orb) pairwise_homographies.append(h_i1_i)

[12]: all_h_matrices = concatenate_homographies(pairwise_homographies, ref=2)

stitched = warp_images_all(all_im, all_h_matrices) plt.figure(figsize=(16, 8)) plt.imshow(stitched);

 

If you implemented concatenate_homographies() correctly, you should get a coherent panorama. Congratulations! You have reached the end of the first assignment.