function [N,M]= extension_proj3(dir)

mask = imread(sprintf('psmImages/%s/%s.mask.tga.ppm', dir, dir));

rows=size(mask,1);
cols=size(mask,2);

pixel_count= rows*cols;
frame_count=12;
I=zeros(pixel_count,frame_count);  
locs=zeros(pixel_count,2);
S=zeros(pixel_count,3);
N=zeros(pixel_count,3);
L=zeros(3,frame_count);
sqrtQdash=zeros(3,3);
Z=zeros(3,3);
sqrtZ=zeros(3,3);
Qdash=zeros(3,3);
A=zeros(3,3);
B=zeros(3,3);
Udash=zeros(pixel_count,3);
scap=zeros(pixel_count,3);
Vdash=zeros(3,frame_count);
lcap=zeros(3,frame_count);
Bdash=zeros(6,1);
D=zeros(6,6);
M=zeros(3,frame_count);


pixels=1;
for i=1:rows
    for j=1:cols
        if(mask(i,j)==255)
            locs(pixels,1)=i;
            locs(pixels,2)=j;
            pixels =pixels + 1;
        end
    end
end
pixel_count=pixels;


for f=1:frame_count
    img=imread(sprintf('psmImages/%s/%s.%d.tga.ppm', dir, dir, (f-1)));
    for p=1:(pixel_count-1)
        I(p,f)= img(locs(p,1),locs(p,2));
    end
end

[U,Q,V] = svd(I,0);

Qdash= Q(1:3,1:3);

%calculate square root of Q and update qdash
Qdash(1,1)=sqrt(Qdash(1,1));
Qdash(2,2)=sqrt(Qdash(2,2));
Qdash(3,3)=sqrt(Qdash(3,3));

scap = U(1:pixel_count,1:3) * Qdash;
lcap = Qdash * V(1:3,1:frame_count);

%Compute matrix A 

    %change this according to the pixel_count and the object
    sk=scap(200:205,1:3);  
    for j=1:6
        D(j,1)=sk(j,1)^2;
        D(j,2)=sk(j,2)^2;
        D(j,3)=sk(j,3)^2;
        D(j,4)=2*sk(j,1)*sk(j,2);
        D(j,5)=2*sk(j,1)*sk(j,3);
        D(j,6)=2*sk(j,2)*sk(j,3); 
    end
K=[1;1;1;1;1;1];
Bdash=K\D;
%B is symmetric
B(1,1)=Bdash(1);
B(2,2)=Bdash(2);
B(3,3)=Bdash(3);
B(1,2)=Bdash(4);
B(2,1)=Bdash(4);
B(1,3)=Bdash(5);
B(3,1)=Bdash(5);
B(2,3)=Bdash(6);
B(3,2)=Bdash(6);
[w,Z,wt]=svd(B,0);
Z(1,1)=sqrt(Z(1,1));
Z(2,2)=sqrt(Z(2,2));
Z(3,3)=sqrt(Z(3,3));
A=w * Z;

A_inv=inv(A);
S= scap * A;
L= A_inv * lcap;

%finding lighting direction matrix
M = L;   

%finding surface normals from S=R*N
 for n=1:pixel_count
     r_n=sqrt((S(n,1)^2)+(S(n,2)^2)+(S(n,3)^2));
     if(r_n~=0)
        for m=1:3
        N(n,m)=S(n,m)/r_n;
        end
     end
   
 end
 
%needle map
X= locs(:,1);
Y= locs(:,2);
NX = N(:,1);
NY = N(:,2);
NZ = N(:,3);
quiver(X,Y,NX,NY,0.2);

%RGB encoded normal map

rgb_image=zeros(340,512,3);
for i=1:(pixel_count-1)
    x=locs(i,1);
    y=locs(i,2);
    rgb_image(x,y,1) = 255*((NX(i)+1)/2);
    rgb_image(x,y,2) = 255*((NY(i)+1)/2);
    rgb_image(x,y,3) = 255*(NZ(i));
end
figure(2),imshow(uint8(rgb_image));
 
 
