function [ normals ] = FindNormals( lightingVectors, sphereFileName, numberOfFiles )

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

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

% 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 + ...

% Matrices for system of equations
A = sparse(3 * height * width * numberOfFiles, 3 * height * width);
b = sparse(3 * height * width * numberOfFiles);
x = sparse(3 * height * width);

bRow = 1;
for img=1:numberOfFiles
    
    % read image
    image = imread(strcat(sphereFileName,'.',num2str(img),'.bmp'));
    for i=1:height
       for j=1:width
           
           % use grayscale pixel value for I
           I = 0.3*image(height,width,1) + 0.59*image(height,width,2) + 0.11*image(height,width,3);
           
           for element=1:3
                % fill A
                colPos = ((i-1) * width * 3) + (j-1) * 3 + 1; % 1 indexing is stupid
                A(bRow, colPos) = 2 * lightingVectors(img,element) * lightingVectors(img,1) * I^2;
                A(bRow, colPos+1) = 2 * lightingVectors(img,element) * lightingVectors(img,2) * I^2;
                A(bRow, colPos+2) = 2 * lightingVectors(img,element) * lightingVectors(img,3) * I^2;
           
                % fill b
                b(bRow) = 2 * lightingVectors(img,element) * I^3;
           
                bRow = bRow + 1;
           end
       end
    end
end

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

normals = x;

end