ACTIVITY 17 - Photometric Stereo

In this activity, we estimated and then extracted the shape of an object from shadow with the use of different sources and shading models.

First, we utilized the images of of synthetic spherical surfaces that are illuminated by a far away point source.

It looks like there are no significant difference among the pictures above but the shading of the images tells much information about the surface of the object. It gives the intensity captured by the camera at point (x,y). These images are captured from the surface of the object with the sources located respectively at
This numbers were put into matrix form where each row is a source and each column is the x,y, and z component of the source.
Now we are given with I and V which is related by
We can solve for the surface normal vector by first getting g with the use of the equation
This is the reflectance of the of the object at the point normal to the surface. The surface normal vector is obtained from g divided by its magnitude
From the surface normals, we computed the elevation z = f(u,v) and the 3D plot of the shape of the object was displayed.

The surface normals (nx, ny, nz) obtained using photometric stereo are related to the partial derivative of f(x,y) as
Then, the surface elevation z at point (u,v) which is given by f(u,v) is evaluated by a line integral
Finally, the 3D plot of the shape of the object was displayed.
The shape of the object with spherical surfaces was successfully extracted and displayed.

The Scilab code used in this activity was shown below. I want to acknowledge the help of Gilbert in working on this activity.

loadmatfile('Photos.mat');

// intensity captured by camera at point (x,y)
I1 = matrix(I1, 1, size(I1, 1)*size(I1, 2));
I2 = matrix(I2, 1, size(I2, 1)*size(I2, 2));
I3 = matrix(I3, 1, size(I3, 1)*size(I3, 2));
I4 = matrix(I4, 1, size(I4, 1)*size(I4, 2));
I = [I1; I2; I3; I4];

// location of the point sources
V1 = [0.085832, 0.17365, 0.98106];
V2 = [0.085832, -0.17365, 0.98106];
V3 = [0.17365, 0, 0.98481];
V4 = [0.16318, -0.34202, 0.92542];
V = [V1; V2; V3; V4];

// calculation of surface normal vector
g = inv(V'*V)*V'*I;
magnitude = sqrt((g(1,:).^2) + (g(2,:).^2) + (g(3,:).^2))+.0001;
n = [];
for i = 1:3
n(i,:) = g(i,:)./magnitude;
end

// computation for the elevation z = f(x,y)
nx = n(1,:);
ny = n(2,:);
nz = n(3,:) + 0.0001;
dfx = -nx./nz;
dfy = -ny./nz;
z1 = matrix(dfx,128,128);
z2 = matrix(dfy,128,128);
Z1 = cumsum(z1,2); // integration from 0 to u
Z2 = cumsum(z2,1); // integration from 0 to v
z = Z1 + Z2;
scf(0);
plot3d(1:128, 1:128, z);