% creazione di data set sintetico per photometric stereo

% funzione test e derivate parziali
f = @(x,y) 0.5*exp(x).*sin(pi*x).*sin(pi*y);
fx = @(x,y) 0.5*exp(x).*(sin(pi*x)+pi*cos(pi*x)).*sin(pi*y);
fy = @(x,y) 0.5*pi*exp(x).*sin(pi*x).*cos(pi*y);
delta = @(x,y) 0.5*exp(x).*((1-2*pi^2)*sin(pi*x)+2*pi*cos(pi*x)).*sin(pi*y);

% vettore luce
intens = .9;
vluce = [1;-1;1]; vluce = intens*vluce/norm(vluce);

% costanti varie
A = 2;			% horizontal size of image in meters (or cm, ...)
r = 101;		% horizontal size of image in pixels
s = 101;		% vertical size of image in pixels
p = r*s;		% pixels number

% grafico della funzione
h = A/(r-1);		% pixel distance
B = (s-1)*h;		% vertical size of image (square pixels)
x = linspace(-A/2,A/2,r);
y = linspace(-B/2,B/2,s);
[X Y] = meshgrid(x,y);
Z = f(X,Y);

figure(1)
%surfl(X,Y,Z)
%colormap gray
surf(X,Y,Z)
axis equal
colorbar
shading interp
pause

% calcolo gradienti e normali
P = fx(X,Y);
Q = fy(X,Y);
K = 1./sqrt(1+P.^2+Q.^2);
P = -P.*K;
Q = -Q.*K;
N = [P(:) Q(:) K(:)]';

% calcola intensita' luminosa e riordina i pixel per creare l'immagine
m = (N' * vluce);
M = reshape(m,s,r);

% mostra immagine
figure(2)
pcolor(M)
shading interp
colormap gray
axis image
axis off
pause

% texture mapping
figure(3)
hh = surfl(X,Y,Z);
axis equal
shading interp
colormap gray
set(hh,'cdata',M)

