sfs.m

matlab source code for shape from shading

function [zmap]=sfs(image)
iter=10;
Sx=0.001;
Sy=0.001;
Sz=1.0;
[szy,szx]=size(image);
Zn=zeros(szy,szx);
Zn1=zeros(szy,szx);
Si1=zeros(szy,szx);
Si=zeros(szy,szx);
Wn = 0.001*0.001; %(regularization constraint)
for i=1:szy
    for j=1:szx
        Zn1(i,j)=0.0;
        Si1(i,j)=1.0;
    end
end

%     if Sx == 0 & Sy == 0
%         Sx=0.001;
%         Sy=0.001;
%     end

Ps = Sx/Sz;
Qs = Sy/Sz;
PQs = 1 + Ps*Ps + Qs*Qs;
E = double(image);
d= max(E);
d = max(d);
for I=1:iter
    E=double(image);

        for i=1:szy
        for j=1:szx
 if (j-1<1) | (i-1 <1) | (j==szx) | (i==szy)
                p=0.0;
                q=0.0;
            else
                p = Zn1(i,j)-Zn1(i,j-1);
                q = Zn1(i,j)-Zn1(i -1,j);
            end
            pq = 1 + p*p + q*q;
%             PQs = 1 + Ps*Ps + Qs*Qs;
            E(i,j)=E(i,j)/d;
            fZ = -1.0*(E(i,j) - max(0.0,(1 + p*Ps + q*Qs)/(sqrt(pq)*sqrt(PQs))));
            dfZ = -1.0 *max(1.0,((Ps+Qs)/(sqrt(pq)*sqrt(PQs)) - ( (p+q)*(1 + p*Ps + q*Qs)/(sqrt(pq*pq*pq)*sqrt(PQs))))); 
            
%              if dfZ==0
%                  dfZ=-1.0;
%              end
            Y = fZ + dfZ*Zn1(i,j);
            
            k = Si1(i,j)*dfZ/(Wn+dfZ*Si1(i,j)*dfZ);
            
            Si(i,j) = (1 - k*dfZ)*Si1(i,j);
            
            Zn(i,j) = Zn1(i,j) + k*(Y-dfZ*Zn1(i,j));
        end
    end
    
    for i=1:szy
        for j =1:szx

Zn1(i,j)=Zn(i,j);
            Si1(i,j)=Si(i,j);
        end
    end
end

zmap = Zn1;
smap = Si1;
return