%COMP 256 Assignment 1
%Adrian Ilie

%create matrices
albedo = zeros(height,width);
normal = zeros(height,width,3);
p = zeros(height,width);
q = zeros(height,width);

%compute the elements
for x = 1:height
    x
    for y = 1:width
        iv = imarray(x,y,:);
        iv = iv(:);%make it a vector
        im = diag(iv);
        g = pinv(im*lightdirs')*(im*iv);
%        g=(im*lightdirs')\(im*iv);
        n=norm(g);
        albedo(x,y) = n;
        if(n~=0)
            normal(x,y,:) = g/n;
            if (normal(x,y,3)~=0)
                p(x,y) = normal(x,y,1)/normal(x,y,3);
                q(x,y) = normal(x,y,2)/normal(x,y,3);
            end
        end
    end
end

p(find(p>10))=10;
p(find(p<-10))=-10;
q(find(q>10))=10;
q(find(q<-10))=-10;

%blur derivatives with a Gaussian
% gs=4;%Gaussian size
% sigma=2.5;
% f=fspecial('gaussian',[2*gs+1 2*gs+1],sigma);
% p=conv2(p,f,'full');
% p=p(1+gs:height+gs,1+gs:width+gs);
% q=conv2(q,f,'full');
% q=q(1+gs:height+gs,1+gs:width+gs);

savepfmraw('1.pfm',p);
savepfmraw('2.pfm',q);

%cap the albedo to white
albedo(find(albedo(:)>1))=1;
figure;
imagesc(albedo);colormap(gray(256));
savepfmraw('a.pfm',albedo);

% %find derivarives
% [fxy,fyx] = myCrossDerivs(p,q);
% 
% %eliminate regions where derivatives do not agree
% sq=(fxy-fyx);
% sq=sq.*sq;
% ind=find(sq>1);%1 is the threshold
% p(ind)=0;
% figure;
% imagesc(p);colormap(gray(256));
% savepfmraw('1a.pfm',p);
% q(ind)=0;
% figure;
% imagesc(q);colormap(gray(256));
% savepfmraw('2a.pfm',q);