Homework 2: Photometric Stereo Solution

Code instructions




The goal of this assignment is to implement shape from shading as described in class. This is also described in shape from shading section (Sec 2.2) in Forsyth and Ponce book (pdf link for this section).




The data directory consists of 64 images each of four subjects from the Yale Face database. The light source directions are encoded in the file names. The code consists of several .m functions. We provide loadFaceImages.m, which loads the images and light source directions. loadFaceImages.m returns the images for the 64 light source directions and an ambient image (i.e., image taken with all the light sources turned off).




Your task is to add some code to the top-level script, evalCode.m, and to fill in the code in prepareData.m, photometriStereo.m, getSurface.m, as explained below. The remaining files are utilities to load the input data and display the output. When testing your code you can change parameters in evalCode.m to change which image set and surface method is used. You can start by setting the subjectName='debug' which creates images from a toy scene. You can debug your code against this before you try the faces. You need to implement:




Preprocess the data: subtract the ambient image from each image in the light source stack, set any negative values to zero, rescale the resulting intensities to between 0 and 1 (they are originally between 0 and 255).



Hint: these operations can be done without using any loops whatsoever. You may want to look into MATLAB's bsxfun function.




Estimate the albedo and surface normals. For this, you need to complete photometricStereo.m, which is a function taking as input the image stack corresponding to the different light source directions and the matrix of the light source directions, and returning an albedo image and surface normal estimates. The
latter should be stored in a three-dimensional matrix. That is, if your original image dimensions are h x w, the surface normal matrix should be h x w x 3, where the third dimension corresponds to the x-, y-, and z-components of the normals. To solve for the albedo and the normals, you will need to set up a linear system as shown in the lecture slides.




Hints:




To get the least-squares solution of a linear system, use MATLAB's backslash operator. That is, the solution to Ax = b is given by x = A\b.




If you directly implement the formulation in our lecture slides, you will have to loop over every image pixel and separately solve a linear system in each iteration. There is a way to get all the solutions at once by stacking the unknown g vectors for every pixel into a 3 x n matrix and getting all the solutions with a single application of the backslash operator. (This is a suggestion, not a requirement.)




You will most likely need to reshape your data in various ways before and after solving the linear system. Useful MATLAB functions for this include reshape and cat.




You may also need to use element-wise operations. For example, for two equal-size matrices X and Y, X .* Y multiplies corresponding elements, and X.^2 squares every element. As before, bsxfun can be a useful function here.







Compute the surface height map by integration. The method is shown in the slides, except that instead of continuous integration of the partial derivatives over a path, you will simply be summing their discrete values. Your code implementing the integration should go in the getSurface.m file. As stated in the slide, to get the best results, you should compute integrals over multiple paths and average the results. You should implement the following variants of integration:



Integrating first the rows, then the columns. That is, your path first goes along the same row as the pixel along the top, and then goes vertically down to the pixel. It is possible to implement this without nested loops using the cumsum function.
Integrating first along the columns, then the rows.
Average of the first two options.
Average of multiple random paths. For this, it is fine to use nested loops. You should determine the number of paths experimentally.



Display the results using functions displayOutput and plotSurfaceNormals. To make things easier, running evalCode.m will call everything automatically.










Report instructions




For full credit, your report should include a section for each of the following questions:




Briefly describe your algorithms for preprocessing, albedo, surface normals, and integration. Include the formulas you have used.



For all four subjects, display your estimated albedos, surface normals, and height maps, with all four integration method. Use displayOutput and plotSurfaceNormals. When inserting results images into your report, you should resize/compress them appropriately to keep the file size manageable -- but make sure that the correctness and quality of your output can be clearly and easily judged. For the height maps, take screenshots of the 3D surfaces from a viewpoint that makes the structure as clear as possible. Feel free to include more than one screenshots from each height map.



Discuss the differences between the four different integration methods.



You should choose only one subject. Specify which subject you choose to evaluate, and discuss which method produces the best results and why.
b. You should also compare the running times of the different approaches on the subject of your choice. For timing, you can use tic and toc functions.




Discuss how the Yale Face data violates the assumptions of the shape-from-shading method covered in the slides. What features of the data can contribute to errors in the results? Feel free to include specific input images to illustrate your points.






Submission instructions




Turn in the following files:




prepareData.m







photometricStereo.m




getSurface.m




report.pdf







Submit your report on Gradescope and your code on Canvas.







Acknowledgements




This homework is based on a similar one made by Lana Lazebnik at UIUC.