Activity 17 - Photometric Stereo

The shape of a surface is estimated and then plotted 3-dimensionally in this activity using photometric stereo.

In photometric stereo, the 3-dimensional form of an object can be determined through captured intensity profiles of that object. The intensity profile is varied through different locations of an illuminating source. A point source is considered here. The location of N point sources is defined in a matrix as follows.

A row represents a point source and each row corresponds to its x, y, and z component. Now, for N number of images of the object, the intensity at (x, y) on the surface is given by the following.

This can be written in a matrix operation shown below.

Note that g is a function of the reflectance of the surface multiplied to the normal vector at (x, y) as shown in the following equation.

Since g is initially unknown, this can be obtained using the matrix equation as follows.

Then, the normal vector of the surface is calculated through the expression below.

Note that the denominator from the equation above is the length of g computed employing Pythagorean theorem in 3-dimension. From the normal vector, the shape of the object can be determined in the following definition of its surface elevation at (x, y).

The expression above is then rewritten as follows.

The gradient of c(x, y, z) shown below is next obtained in order to relate the estimated surface normals n_x, n_y, and n_z using photometric stereo and the partial derivative of f.

Meanwhile, the relation between the surface normals and the partial derivative of f is the following.

Finally, the object's surface elevation at point (u, v) is evaluated through the line integral as follows.

A Matlab file containing four images of a surface is loaded in Scilab using the function loadmatfile. These images are shown below.

Each image above is illuminated by a point source located respectively at
V_1 = [0.085832, 0.17365, 0.98106]
V_2 = [0.085832, -0.17365, 0.98106]
V_3 = [0.17365, 0, 0.98481]
V_4 = [0.16318, -0.34202, 0.92542]

These four point source locations are then placed in matrix V with size 4 x 3 as well as the four intensity profiles in matrix I with size 4 x (128*128) or 4 x 16384. Note that each loaded intensity profile has size 128 x 128. The size of matrix g and surface normal n obtained is thus 3 x 16384. To prevent division by zero, a very small value 0.000001 is added to every element of the length or magnitude of g and to all elements of n_z. The calculated partial derivatives of f with respect to x and y are recasted into a matrix with size 128 x 128 for evaluation of the surface elevation at (u, v). The line integrals are operated using cumsum of Scilab where the limit of integration from 0 to u means the cumulative summation along the column (from 1st column to 128th column) of the partial derivative of f with respect to x. Conversely, the limit of integration from 0 to v is the cumulative summation along the row (from 1st row to 128th row) of the partial derivative of f with respect to y.
The following figure is then the plotted surface of the object in 3-dimension.

Apparently, the shape of the surface is successfully reconstructed 3-dimensionally. From the 3-D plot, it is revealed that the object is hemispherical with four divided quadrants as observed from the intensity profiles.

Since I am able to show the shape of the given surface successfully using photometric stereo, I grade myself 10/10.

I worked with Gary in completing this activity.

Appendix
The whole Scilab code utilized in this activity is as follows.

loadmatfile('Photos.mat');

//scf(0);
//imshow(I1, []);
//scf(1);
//imshow(I2, []);
//scf(2);
//imshow(I3, []);
//scf(3);
//imshow(I4, []);

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];

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];

small_number = 0.000001;
g = inv(V'*V)*V'*I;
magnitude = sqrt((g(1,:).^2) + (g(2,:).^2) + (g(3,:).^2)) + small_number;

for i = 1:3
n(i,:) = g(i,:)./magnitude;
end

n_x = n(1,:);
n_y = n(2,:);
n_z = n(3,:) + small_number;
partial_f_to_x = -n_x./n_z;
partial_f_to_y = -n_y./n_z;
partial_f_to_x = matrix(partial_f_to_x, 128, 128);
partial_f_to_y = matrix(partial_f_to_y, 128, 128);
z = cumsum(partial_f_to_x, 2) + cumsum(partial_f_to_y, 1);

//scf(4);
//plot3d(1:128, 1:128, z);