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