Shape from shading Pentland approach

Pentland used the linear approximation of the reflectance function in terms of the surface gradient, and applied a Fourier transform to the linear function to get a closed form solution for the depth at each point. (Pentland takes the Fourier transform of both sides of the equation). (Pentland, A., “Shape Information From Shading: A Theory About Human Perception,” Computer Vision., Second International Conference on , vol., no., pp.404-413, 5-8 Dec 1988.)

Algorithm of shape from shading using Pentland approach as follows:

clear all;
close all;
% read the image
E = imread('Nike.jpg');

% making sure that it is a grayscale image
if size(E,3)==3
    E = rgb2gray(E);
end
E = double(E);
% normalizing the image to have maximum of one
E = E ./ max(E(:));
% first compute the surface albedo and illumination direction
[albedo,I,slant,tilt] = estimate_albedo_illumination (E);

% compute the fourier transform of the image
Fe = fft2(E);

% wx and wy
[M,N] = size(E);
[x,y] = meshgrid(1:N,1:M);
wx = (2.* pi .* x) ./ M;
wy = (2.* pi .* y) ./ N;

% Using the estimated illumination direction
Fz =Fe./(-1i.*wx.*cos(tilt).*sin(slant)-1i.*wy.*sin(tilt).*sin(slant));

% Compute the inverse Fourier transform to recover the surface.
Z = abs(ifft2(Fz));

% visualizing the result
figure;
surfl(Z);
shading interp;
colormap gray(256);
lighting phong;

“estimate_albedo_illumination” function:

function [albedo,I,slant,tilt] = estimate_albedo_illumination (E)
% E: the image
% normalizing the image to have maximum of one
%E = E ./ max(E(:));
% compute the average of the image brightness
Mu1 = mean(E(:));
% compute the average of the image brightness square
Mu2 = mean(mean(E.^2));

% now lets compute the image's spatial gradient in x and y directions
[Ex,Ey] = gradient(E);
% %imshow(Ex)
% %imshow(Ey)
% figure(1),
% subplot(1,2,1), imshow(Ex)
% subplot(1,2,2), imshow(Ey)

% normalize the gradients to be unit vectors
Exy = sqrt(Ex.^2 + Ey.^2);
nEx = Ex ./(Exy + eps); % to avoid dividing by zero
nEy = Ey ./(Exy + eps);
% % imshow(nEx)
% % imshow(nEy)
% figure(2),
% subplot(1,2,1), imshow(nEx)
% subplot(1,2,2), imshow(nEy)

% computing the average of the normalized gradients
avgEx = mean(nEx(:));
avgEy = mean(nEy(:));

% now lets estimate the surface albedo
gamma = sqrt((6 *(pi^2)* Mu2) - (48 * (Mu1^2)));
albedo = gamma/pi;
% estimating the slant
slant = acos((4*Mu1)/gamma);
% estimating the tilt
tilt = atan(avgEy/avgEx);
if tilt < 0
    tilt = tilt + pi;
end
% the illumination direction will be:
I = [cos(tilt)*sin(slant) sin(tilt)*sin(slant) cos(slant)];
end

Leave a Reply

Your email address will not be published. Required fields are marked *