function [ normals, albedos ] = FindNormals( lightingVectors, fileName, numberOfFiles )

image = imread(strcat(fileName,'.1.bmp'));
mask = imread(strcat(fileName,'.mask.bmp'));
imgSize = size(image);
height = imgSize(1);
width = imgSize(2);

% g = kd * N
normals = zeros(height,width,3);
albedos = zeros(height,width);

% Find the normal by setting up a system of equations and solving
% Minimizing the following objective function by differentiating with
% respect to g and setting equal 0
%
% Objective function: Q = Sigma_i( (I_i^2 - I_i * L_i * transpose(g))^2 )
%                     Q = (I_i^2 - Ii* (L_i1*g_i1 + L_i2*g_i2 +
%                     L_i3*g_i3))^2 + ...

% setup matrices for solving the linear system
A = sparse([],[],[],3 * numberOfFiles * height * width, 3 * height * width,9*numberOfFiles*height*width);
b = zeros(3 * numberOfFiles * height * width, 1);

bRow = 1;

for img=1:numberOfFiles

    aCol = 1;
    % read image
    image = imread(strcat(fileName,'.',num2str(img),'.bmp'));
    for i=1:height
       for j=1:width

            if mask(i,j,1) > 0

                % use grayscale pixel value for I
                % Might need to change range to 1-256 from 0-255
                I = 0.3*double(image(i,j,1)+1) + 0.59*double(image(i,j,2)+1) + 0.11*double(image(i,j,3)+1);

                for element=1:3

                    % fill A
                    A(bRow, aCol) = 2 * I^2 * lightingVectors(img,element) * lightingVectors(img,1);
                    A(bRow, aCol+1) = 2 * I^2 * lightingVectors(img,element) * lightingVectors(img,2);
                    A(bRow, aCol+2) = 2 * I^2 * lightingVectors(img,element) * lightingVectors(img,3); 

                    % fill b
                    b(bRow) = 2 * I^3 * lightingVectors(img,element);
                    bRow = bRow + 1;
                end
            else
                bRow = bRow + 3;
            end
            aCol = aCol + 3;
            
        end
    end
end

% Solve System using <strike>magic</strike> math
x = A\b;

display('Done with first bit.')

xIndex = 1;
for i=1:height
    for j=1:width
        for k=1:3
            normals(i,j,k) = x(xIndex);
            xIndex = xIndex + 1;
        end
        g = [normals(i,j,1),normals(i,j,2),normals(i,j,3)];
        albedos(i,j) = norm(g);
        if albedos(i,j) ~= 0
            normals(i,j,:) = normals(i,j,:)/albedos(i,j);
        end
    end
end

end