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

counter = 0;
for i=1:height
   for j=1:width
       counter = counter + 1;
       
        % setup matrices for solving the linear system
        A = zeros(3 * numberOfFiles, 3);
        b = zeros(3 * numberOfFiles, 1);

        bRow = 1;
        for img=1:numberOfFiles
            
            % This is where setup used to be performed

            % read image
            image = imread(strcat(sphereFileName,'.',num2str(img),'.bmp'));

            % 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, 1) = 2 * I^2 * lightingVectors(img,element) * lightingVectors(img,1);
                A(bRow, 2) = 2 * I^2 * lightingVectors(img,element) * lightingVectors(img,2);
                A(bRow, 3) = 2 * I^2 * lightingVectors(img,element) * lightingVectors(img,3); 
                
                % fill b
                b(bRow) = 2 * I^3 * lightingVectors(img,element);
                
                bRow = bRow + 1;
            end
        end
        
        % Solve System using <strike>magic</strike> math
        x = A\b;

        g(i,j,:) = x;
    end
end

normals = g;

end