## Photometric Stereo

Tushar Khot
tushar@cs.wisc.edu

Since Matlab cant handle tga files, I convert TGA files into BMP files on the fly using ImageMagick's convert commandline tool.
Hence this code will need linux machine with "convert" command available.
• Calibration
For calibrating the direction of the light using the chrome ball, I
1. Determined the center of the ball by averaging all pixels in mask with value=255, say [xc yc]
2. Similarly determined the center of the bright spot, say [xl yl]
3. Determined the radius by averaging the distance between points with value in mask greater than 0 but less than 255(boundary points). Say r.
4. N = [xl-xc yl-yc sqrt(r^2 - (xl-xc)^2 - (yl-yc)^2)] normalized to norm=1. R=[0 0 -1]
5. L = 2*N*R' - R
Code: lightDirection.m
• Normals
For detemining the normals, I
1. For each pixel in ever image, combine the values for each channel.
2. Try to minimize the objective function, (w*I - w*L*g) where w = I.
3. I even attempted using the hat function as the weight just like the one used in HDR, but that had almost no impact on the output.
4. N = (1/norm(g))*g. k = norm(g).
Code: getNormals.m
• Color Albedo
For obtaining the color albedo, I used the equation I.J/J.J as given on the project page. This seems to work well.
Code: getAlbedo.m
• Least Square Fitting
For solving this problem, I
1. Created an index of pixel co-ordinates where mask value is non-zero.
2. Created two sparse matrices A and b
3. A and b contain rows corresponding to gradient in each direction as per the equations on the project page. Since every row in A has only two non-zero entries, the matrix is very sparse.
4. A has one row to fix one pixel to a height > 0 so that the entire surface appears to be above 0.

Code: computeZ.m
• Extra
There were multiple small changes that I attempted.
1. Used a hat weighing function for the normal computation. Causes some peaks - look at the eyes.

2. Used no weighing function for the normal computation. Still doesn't get rid off the peaks.

3. Used a hat weighing function for surface estimation. So every row in the matrix A and b was scaled by the intensity. Creates wierd peaks in the image.

4. Used the intensity itself as the weighing function. This also creates peaks in the image.

5. To reduce the peaks in the surface, I added another minimization term to minimize the difference between the adjacent points i.e. minimize the sum of squared differences between neighbouring pixels. This still keeps A sparse.
This works well to reduce noise. This is what I have used in all the output images.
6. I used MATLAB's builtin "pcg" function to solve the surface estimation equation with A'Az = A'b.
The result was not very different from the output given by \.
NOTE: They both use the squared difference terms mentioned in the previous point. Also PCG with 100 iterations gave a very noisy answer, so I used 1000 iterations.  With \ With pcg (1000 iterations)

7. Sometimes, the Z-values gradually increase to high values. Hence, I tried adding a constraint to minimize the norm of the Z values similar to HDR.
But this doesn't help much because I was raising the entire surface. Also without that condition, even if the norm is small, 10^-3 would appear as a peak if all values are 10^-3. The norm in this case would still be low.
• How to Run
0. Prereq:
Copy the psmImages folder into the directory that contains the code.
Make sure convert is a valid command or add /usr/bin/ to \$PATH.
1. stereo e.g. stereo buddha
2. Various parameters that can be changed from stereo.m
• usePCG - Use PCG for solving surface estimation.
• useSquaredDiff - Use the squared differences between adjacent Z values.
• useWeightsSurf - Use weights for estimating the surface in computeZ.m
• useWeightsNormal - Use weights for estimating the normal in getNormals.m

• Output
 DataSet RGB-encoded Normals(absolute values) Albedo-map Surface Buddha Owl Horse Rock Gray Cat

References:
1. Woodham, Optical Engineering, 1980, Photometric Stereo.