Project 3b: Panoramic Mosaic Stitching

Assigned: Friday, February 19
Due: Monday, March 1 (1:30 PM)

---

In this part of the project, you will use your feature detection and matching code from project 3a to combine a series of photographs into a 360° panorama. Your software will automatically align the photographs (determine their overlap and relative positions) and then blend the resulting photos into a single seamless panorama. You will then be able to view the resulting panorama inside an interactive Web viewer.

To start this component, you will be supplied with some test images and skeleton code that will guide you.
Taking the Pictures
Code to Write
Creating the Panorama
Debugging
Extra Credit

Submission
Panorama Links

Taking the Pictures

Each group will be checking out a panorama kit (camera, tripod, and Kaidan head). Sign up for a kit here.

Remember to bring extra batteries with you – these cameras drain batteries.

1. Take a series of photos with a digital camera mounted on a tripod. You will borrow the Kaidan head that lets you make precise rotations and the Canon PowerShot A10 camera for this purpose. Here is an explanation of how to use the equipment. Please read it before you go out to shoot. For best results, overlap each image by 50% with the previous one, and keep the camera level using the levelers on the Kaidan head.
2. Also take a series of images with a handheld camera.  You can use your own or use the Canon PowerShot A10 camera used in Step 1 above. If you are using the Canon camera, it has a “stitch assist” mode you can use to overlap your images correctly – this only works in regular landscape mode.  If you are using your own camera, you have to estimate its focal length.  The simplest way to do this is through the EXIF tags of the images, as described by Noah Snavely here.  Alternatively, you can use a camera calibration toolkit to get more precise focal length and radial distortion coefficients.  Finally, Brett Allen describes one creative way to measure rough focal length using just a book and a box.  The parameters for the class cameras are given below. You need to keep the camera zoomed out all the way for the following focal length to be valid.

Camera	resolution	focal length	k1	k2
Canon Powershot A10, tag CS30012716	480x640	678.21239 pixels	-0.21001	0.26169
Canon Powershot A10, tag CS30012717	480x640	677.50487 pixels	-0.20406	0.23276
Canon Powershot A10, tag CS30012718	480x640	676.48417 pixels	-0.20845	0.25624
Canon Powershot A10, tag CS30012927	480x640	671.16649 pixels	-0.19270	0.14168
Canon Powershot A10, tag CS30012928	480x640	674.82258 pixels	-0.21528	0.30098
Canon Powershot A10, tag CS30012929	480x640	674.79106 pixels	-0.21483	0.32286
test images	384x512	595 pixels	-0.15	0.0

3. Make sure the images are right side up (rotate the images by 90° if you took them in landscape mode), and reduce them to a more workable size (480x640 recommended). You can use external software such as PhotoShop or the Microsoft Photo Editor to do this. Or you may want to set the camera to 640x480 resolution from the start, by following the steps below:
1. Turn the mode dial on the back of the camera to one of the 3 shooting modes--auto (camera icon), manual (camera icon + M) or stitch assist (overlaid rectangles).
2. Press MENU button.
3. Press the left/right arrow to choose Resolution, then press SET.
4. Press the left/right arrow and choose S (640x480).
5. Press MENU again.

(Note: If you are using the skeleton software, save your images in (TrueVision) JPEG format (.jpg), since this is the only format the skeleton software can read. Also make sure the aspect ratio of the image (width vs. height) is either 4:3 or 3:4 (480x640 will do) which is the only aspect ratio supported by the skeleton software.)

Code to Write (80%)

1. (15%) WarpSphericalField(srcImgName, warpImgName,f,k1, k2) (file: WarpSpherical.m, routine: warpSphericalField)

Warp each image into spherical coordinates.

Compute the inverse map to warp the image by filling in the skeleton code in the warpSphericalField routine to:

1. convert the given spherical image coordinate into the corresponding planar image coordinate using the coordinate transformation equation from the lecture notes
2. apply radial distortion using the equation from the lecture notes

(Note: You will have to use the focal length f estimates for the half-resolution images provided above (you can either take pictures and save them in small files or save them in large files and reduce them afterwards) . If you use a different image size, do remember to scale f according to the image size.)

2. (15%) [M]=alignPair(locs1,locs2,match_idx, nRANSAC, RANSACthresh) (file: alignPair.m, routine: alignPair)

alignPair takes two feature sets, locs1 and locs2, and match_idx , the list of feature matches obtained from the feature detecting and matching component (described in the first part of the project),  and estimates and inter-image transform matrix M.  

alignPair uses RANSAC (RAndom SAmpling Consensus) to pull out a minimal set of feature matches (one match for this project), estimates the corresponding motion (alignment) and then invokes countInliers to count how many of the feature matches agree with the current motion estimate.  After repeated trials, the motion estimate with the largest number of inliers is used to compute a least squares estimate for the motion leastSquaresFit, which is then returned in the motion estimate M.

You will have to fill in the missing code in alignPair to:

3. Randomly select a valid matching pair and compute the translation between the two feature locations.
4. Call countInliers to count how many matches agree with this estimate.
5. Repeat the above random selection nRANSAC times and keep the estimate with the largest number of inliers.
6. Write the body of countInliers to count the number of feature matches where the SSD distance after applying the estimated transform (i.e. the distance from the match to its correct position in the image) is below the threshold. (and don’t forget to create the list of inlier ids.)
7. Write the body of leastSquaresFit, which for the simple translational case is just the average displacement between the matching feature positions.

3. (10%)[count,inliers]=countInliers(locs1,locs2,match_idx,M, RANSACthresh) (file: countInliers.m, routine: countInliers)

CountInliers computes the number of matches that have a distance below RANSACthresh is computed.  It also returns a list of inlier match ids.

4. [M]=leastSquaresFit(locs1,locs2,match_idx,inliers) (file: leastSquareFit.m, routine: leastSquareFit)

LeastSquaresFit computes a least squares estimate for the translation using all of the matches previously estimated as inliers.  It returns the resulting translation estimate in the last column of M.

5. (15%) BlendImages(ipv, blendWidth,panoImgName) (file: BlendImages.m, routine: BlendImages)

Given the warped images and their relative displacements, figure out how large the final stitched image will be and their absolute displacements in the panorama (BlendImages).

It then calls AccumulateBlend and NormalizeBlend to resample each image to its final location and blend it with its neighbors.

Crop the resulting image to make the left and right edges seam perfectly (BlendImages). The horizontal extent can be computed in the previous blending routine since the first image occurs at both the left and right end of the stitched sequence (draw the “cut” line halfway through this image).  Use a linear warp to the mosaic to remove any vertical “drift” between the first and last image.  This warp, of the form y' = y + ax, should transform the y coordinates of the mosaic such that the first image has the same y-coordinate on both the left and right end.  Calculate the value of 'a' needed to perform this transformation.

6. (10%) [acc]=AccumulateBlend(img, acc, blendWidth) (file: AccumulatedBlend.m, routine: AccumulatedBlend)

Add up the (α premultiplied) RGBα values at each pixel. Try a simple feathering function as your weighting function (see mosaics lecture slide on "feathering") (this is a simple 1-D version of the distance map described in [Szeliski & Shum]).  For extra credit, you can try other blending functions or figure out some way to compensate for exposure differences.

7. (5%) [img]=NormalizeBlend(acc) (file: NormalizeBlend.m, routine: NormalizeBlend)

Divide each pixel’s accumulated RGB by its α value. Remember to set the alpha channel of the resultant panorama to opaque!

What to Turn In (20%)

In addition to your source code, turn in a web page describing your approach and results.  In particular:

1. Who are the members in your group and how did you divide the work?
2. At least three panoramas:  (1) the test sequence, (2), one from the Kaidan head, and (3) one from a hand-held sequence.  Each panorama should be shown as (1) a low-res inlined image on the web page, (2) a link that you can click on to show the full-resolution .jpg file, AND (3) embedded in a viewer as described in the instructions above.
3. A short description of what worked well and what didn’t.  If you tried several variants or did something non-standard, please describe this as well.
4. Describe any extra credit items you did.

The web-page should be placed in the artifact directory along with all the images in JPEG format. If you are unfamiliar with HTML you can use any web-page editor such as FrontPage, Word, or Visual Studio 7.0 to make your web-page.  Here are some webpage design tips.

Creating the Panorama

1. Use the above program you wrote to warp/align/stitch images into the resulting panorama.
• Use the following to remove the radial distortion and warp the image pano1_0008.png into spherical coordinates with focal length = 595 and radial distortion coefficients k1=-0.15 and k2=0:
  
  WarpSphericalField('pano1_0008.png', 'warp08.png',595,-0.15,0.0);
• Then, use the feature detecting and matching component to compute the features in the warped images. To align two feature sets locs1 and locs2, with a set of matches in match_idx, using 200 iterations of RANSAC with an outlier threshold distance of 1 pixel:
  
  [M]=alignPair(locs1,locs2,match_idx, 200, 1);

Sample SIFT features and matches have been provided to you.

• Run the previous step for all adjacent pairs of images and save the output into a separate file pairlist.txt which may look like this:
  
      warp1.png warp2.png 213.49 -5.12
      warp2.png warp3.png 208.19 2.82
      ......
      warp9.png warp1.tga 194.76 -3.88

A helper function ipvwrite.m that takes a cell array of the image names and a cell of the output Ms from alignPair as input and outputs ipv for BlendImages is provided.

• Then stitch the images into the final panorama pano.png:
  
  BlendImages(ipv, blendWidth,'pano.png');

2. Convert your resulting image to JPEG format (Photoshop and other tools in the lab can help you with this) and paste it on a Web page along with code to run the interactive viewer. Click here for instructions on how to do this.

Debugging Guidelines

You can use the test results included in the images/ folder to check whether your program is running correctly. Comparing your output to that of the sample solution is also a good way of debugging your program.

1. Testing the warping routines:
• In the images/ folder in the skeleton code, a few example warped images are provided for test purposes. See if your program produces the same output.
2. Testing the alignment routines:
• A few example alignment results are provided in the file pairlist2/4.txt.
• To test alignPair only, try passing in an image that has been cropped with two different rectangles (and maybe rotated by a tiny amount, say 2 degrees).
3. Testing the blending routines:
• An example panorama is included in the images/ folder. Compare the resulting panorama with this image.

Extra Credit

Here is a list of suggestions for extending the program for extra credit. You are encouraged to come up with your own extensions. We're always interested in seeing new, unanticipated ways to use this program!

• (*) Although the feature-based aligner gives sub-pixel motion estimation (because of least squares), the motion vectors are rounded to integers when blending the images into the mosaic in BlendImages.cpp. Try to blend images with sub-pixel localization. 
• (*) Sometimes, there exists exposure difference between images, which results in brightness fluctuation in the final mosaic. Try to get rid of this artifact. 
• (*) Try shooting a sequence with some objects moving.  What did you do to remove “ghosted” versions of the objects?
• (*) Try a sequence in which the same person appears multiple times, as in this example.
• (**) Implement a better blending technique, e.g., pyramid blending, Poisson blending, or graph cuts.

Submission

To turn in your project, submit a single zip file to the dropbox at:

https://catalysttools.washington.edu/collectit/dropbox/iansimon/8903

The zip file should be named project3b.zip and contain all of your M-files and your assignment writeup as an HTML, PDF, or Word document. Make sure you include all the images you used for face finding and cropping (other than the provided cropped images), and your plots. Also remember the average face and eigenfaces in the first experiment. Make sure any images that need it are normalized so they can be seen clearly. Finally, make sure to fully document any extra credit on this web page.

---

Panorama Links

• Panoramas.dk: weekly archive of full-screen, high-quality panoramas worldwide
• VR Seattle: Seattle & Washington panorama
• Matt Brown's Autostitch page.
• Peru Panoramas
• A complete set of test images