clear all;
shadow_image = cell(1,7);
% Import the images here
% Thinking about flipping the images up_side_down
% into the ordinary x_y plane coordinate

[Height, Width] = size(shadow_image{1});

% Build your own V here
% V = 255*[  1      0      0;
%            0      1      0;
%            1      0      0;
%            0      1      0;
%            1      0      0;
%            0      1      0;
%            1      0      0;            ];

shape_r = 1; % Grid of reconstruction, to save time set it 2~4 during debuging 

% Normal calculation
for x =  1 : Width/shape_r
    for y =  1 : Height/shape_r
        for k = 1 : 7
            i(k,1) = shadow_image{k}(shape_r*y,shape_r*x);
        end
        I = diag(i);
        Ii = I*i;
        IV = I*V;
        if(det(IV'*IV)~=0)          
            temp = inv(IV'*IV)*IV'*Ii;
            g(y,x,:) = temp;
            if(abs(norm( temp))<10)
                albedo(y,x) = norm( temp);                     
                n(y,x,:) = g(y,x,:)./albedo(y,x);
                p(y,x) = -n(y,x,1)/n(y,x,3);
                q(y,x) = -n(y,x,2)/n(y,x,3);   
            else
                temp = [0;0;0];
                g(y,x,:) = temp;
                albedo(y,x) = 0;     
                n(y,x,:) = temp;
                p(y,x) = 0;   q(y,x) = 0;      
            end
        else
            temp = [0;0;0];
            g(y,x,:) = temp;
            albedo(y,x) = 0;     
            n(y,x,:) = temp;
            p(y,x) = 0;   q(y,x) = 0;        
        end
    end
end

for x =  2 : Width/shape_r
    for y =  2 : Height/shape_r
        if(abs((p(y,x)-p(y-1,x))-(q(y,x)-q(y,x-1)))>5)
            p(y,x)=0;
            q(y,x)=0;
        end
    end
end