function relighting(myPath, myLightRes, myBlocksW, myBlocksH, myBlockSize, myRank, myRank2)
% Process first and second stages of PCA for the input set of images and
% given parameters.
%
% ATTENTION : This function requires a lot of memory, take care with the
% parameters that you choose.
% 
% INPUT
% Use as:   relighting(myPath, myLightRes, myBlocksW, myBlocksH,
% myBlockSize, myRank, myRanks2)
% 
% Where:
%    myPath : Path to the set of input images.
%             Example: 'c:/senrique/relighting/nicoleb6x256/final/'
%             Path must include the final slash to work properly.
%             Input images should be jpg files with the name:
%             'file(f)(uu)(vv).jpg' -without parenthesis-
%             .
%             Where:
%               f : Face, 0: +X, 1: -X, 2: +Y, 3: -Y, 4: +Z, 5: -Z.
%                   +X points to the right, +Y points to the sky, +Z points
%                   to you.
%               uu, vv : 00..myLightRes-1
%                        if f=0 then uu is from +Y to -Y,
%                                    vv is from +Z to -Z
%                           f=1 then uu is from +Y to -Y,
%                                    vv is from -Z to +Z
%                           f=2 then uu is from -Z to +Z,
%                                    vv is from -X to +X
%                           f=3 then uu is from +Z to -Z,
%                                    vv is from -X to +X
%                           f=4 then uu is from +Y to -Y,
%                                    vv is from -X to +X
%                           f=5 then uu is from -Y to +Y,
%                                    vv is from -X to +X
%               Example filename: file10314.jpg
%    myLightRes : Each face was divided in a grid of myLightRes x
%                 myLightRes to generate the corresponding picture with a
%                 point light source. There are 6 x myLightRes x myLightRes
%                 image files.
%                 Example: 16
%    myBlocksW : Each image is composed by myBlocksE of width.
%                Example: 32
%    myBlocksH : Each image is composed by myBlocksH of height.
%                Example: 32
%    myBlockSize : Each block is composed of myBlockSize x myBlockSize
%                  pixels.
%                  Example: 16
%                  (myBlocksW * myBlockSize) should be the image width.
%                  (myBlocksH * myBlockSize) should be the image height.
%                  In this example, width = height = 32 x 16 = 512
%    myRank : Rank to use in the first PCA computation.
%             Example: 10
%    myRank2 : Rank to use in the second PCA computation.
%              Example: 200
%
% OUTPUT
% The following files will be generated in the working directory:
%
% - Files for input in the Relighting Framework (Real-Time):
%   Er.txt
%   Eg.txt
%   Eb.txt
%   Ur.txt
%   Ug.txt
%   Ub.txt
%   Vr.txt
%   Vg.txt
%   Vb.txt
%   --> see LSD paper for details about the meaning of these files.
% 
% - Additional files:
%   status.mat : Some variables used for the calculation and for debugging
%                purposes, like initial and final time, elapsed time, etc.
%   Ir.mat, Ig.mat, Ib.mat : All the original images packed in a big
%                            matrix.
%   For each .txt file generated, the corresponding .mat file.
%
% ---
% Sebastian Enrique - senrique@cs.columbia.edu
% Columbia University - January 2004

% FULL CUBE MAP
% clear;

% CONSTANT DEFINITIONS
%myPath = 'c:/senrique/relighting/nicoleb6x256/final/';  % Path for the images
%myLightRes = 16;
%myBlocks = 32;                       % Blocks square = qty images
%myBlockSize = 16;                    % Blocks size square in pixels
%myRank = 10;
%myRank2 = 200;

m = myBlocksW * myBlocksH;
n = myLightRes * myLightRes;% * 6;
nfc = n * 6;
p = myBlockSize * myBlockSize;
b = myRank;
q = myRank2;

myMessage = sprintf('STARTING!')
myClockStart = fix(clock)
save status;

% =========================================================================
%  READ IMAGES DATA: RED
% -------------------------------------------------------------------------
Ir = uint8([]);
Ir(m,p,nfc) = 0;
i1 = 0;
for e=1:6;
    % For each image in the set
    for h=1:myLightRes;
        for j=1:myLightRes;
            myMessage = sprintf('RED - FACE %.1u - Processing blocks for image [%u] file%.1u%.2u%.2u.jpg',e-1,i1,e-1,h-1,j-1)
            i1 = i1 + 1;
            filename = sprintf([myPath,'file%.1u%.2u%.2u.jpg'],e-1,h-1,j-1);
            image = imread(filename,'jpg');
            i2 = 0;
            % for each block in the image
            for f=1:myBlocksH;
                for g=1:myBlocksW;
                    i2 = i2 + 1;
                    % For each pixel in the block
                    i3 = 0;
                    for i=1:myBlockSize;
                        for k=1:myBlockSize;
                            i3 = i3 + 1;
                            Ir(i2,i3,i1) = image((f-1)*myBlockSize+i,(g-1)*myBlockSize+k,1);
                        end;
                    end;
                end;
            end;
        end;
    end;
end;

save Ir Ir;
clear;
pack;
% =========================================================================

% =========================================================================
%  READ IMAGES DATA: GREEN
% -------------------------------------------------------------------------
load status;
Ig = uint8([]);
Ig(m,p,nfc) = 0;
i1 = 0;
for e=1:6;
    % For each image in the set
    for h=1:myLightRes;
        for j=1:myLightRes;
            myMessage = sprintf('GREEN - FACE %.1u - Processing blocks for image [%u] file%.1u%.2u%.2u.jpg',e-1,i1,e-1,h-1,j-1)
            i1 = i1 + 1;
            filename = sprintf([myPath,'file%.1u%.2u%.2u.jpg'],e-1,h-1,j-1);
            image = imread(filename,'jpg');
            i2 = 0;
            % for each block in the image
            for f=1:myBlocksH;
                for g=1:myBlocksW;
                    i2 = i2 + 1;
                    % For each pixel in the block
                    i3 = 0;
                    for i=1:myBlockSize;
                        for k=1:myBlockSize;
                            i3 = i3 + 1;
                            Ig(i2,i3,i1) = image((f-1)*myBlockSize+i,(g-1)*myBlockSize+k,2);
                        end;
                    end;
                end;
            end;
        end;
    end;
end;

save Ig Ig;
clear;
pack;
% =========================================================================

% =========================================================================
%  READ IMAGES DATA: BLUE
% -------------------------------------------------------------------------
load status;
Ib = uint8([]);
Ib(m,p,nfc) = 0;
i1 = 0;
for e=1:6;
    % For each image in the set
    for h=1:myLightRes;
        for j=1:myLightRes;
            myMessage = sprintf('BLUE - FACE %.1u - Processing blocks for image [%u] file%.1u%.2u%.2u.jpg',e-1,i1,e-1,h-1,j-1)
            i1 = i1 + 1;
            filename = sprintf([myPath,'file%.1u%.2u%.2u.jpg'],e-1,h-1,j-1);
            image = imread(filename,'jpg');
            i2 = 0;
            % for each block in the image
            for f=1:myBlocksH;
                for g=1:myBlocksW;
                    i2 = i2 + 1;
                    % For each pixel in the block
                    i3 = 0;
                    for i=1:myBlockSize;
                        for k=1:myBlockSize;
                            i3 = i3 + 1;
                            Ib(i2,i3,i1) = image((f-1)*myBlockSize+i,(g-1)*myBlockSize+k,3);
                        end;
                    end;
                end;
            end;
        end;
    end;
end;

save Ib Ib;
clear;
pack;
% =========================================================================

% =========================================================================
%  PROCESS SVD 1st STAGE FOR RED
% -------------------------------------------------------------------------
load status;

Lr = [];
Er = [];

load Ir;
Er(p*m, myRank) = 0;
Lr(myRank*m, nfc) = 0;

i2 = 0;
% for each block in the image
for f=1:myBlocksH;
    for g=1:myBlocksW;
        myMessage = sprintf('RED: Processing SVD for blocks [%.2u,%.2u]',f-1,g-1)
        i2 = i2 + 1;
        
        Ijr = double(reshape(Ir(i2,:,:),p,nfc)) / 255.0;
       
        [Ejr, Sjr, Cjr] = svd(Ijr,0);
        clear Ijr;
        Ejr = Ejr(1:p, 1:myRank);
        Sjr = Sjr(1:myRank, 1:myRank);
        Cjr = Cjr(1:nfc, 1:myRank);
        Ljr = Sjr * Cjr';
        clear Cjr;
        clear Sjr;
 
        Er( (i2-1)*p + 1 : (i2-1)*p + p, 1:myRank) = Ejr(:,:);
        clear Ejr;

        Lr( (i2-1)*myRank + 1 : (i2-1)*myRank + myRank, 1:nfc ) = Ljr(:,:);
    end;
    pack;
end;

save Lr Lr;
save Er Er;

clear;
pack;
% =========================================================================

% =========================================================================
%  PROCESS SVD 1st STAGE FOR GREEN
% -------------------------------------------------------------------------
load status;

Lg = [];
Eg = [];

load Ig;
Eg(p*m, myRank) = 0;
Lg(myRank*m, nfc) = 0;

i2 = 0;
% for each block in the image
for f=1:myBlocksH;
    for g=1:myBlocksW;
        myMessage = sprintf('GREEN: Processing SVD for blocks [%.2u,%.2u]',f-1,g-1)
        i2 = i2 + 1;
        
        Ijg = double(reshape(Ig(i2,:,:),p,nfc)) / 255.0;

        [Ejg, Sjg, Cjg] = svd(Ijg,0);
        clear Ijg;
        Ejg = Ejg(1:p, 1:myRank);
        Sjg = Sjg(1:myRank, 1:myRank);
        Cjg = Cjg(1:nfc, 1:myRank);
        Ljg = Sjg * Cjg';
        clear Cjg;
        clear Sjg;
 
        Eg( (i2-1)*p + 1 : (i2-1)*p + p, 1:myRank) = Ejg(:,:);
        clear Ejg;

        Lg( (i2-1)*myRank + 1 : (i2-1)*myRank + myRank, 1:nfc ) = Ljg(:,:);
    end;
    pack;
end;

save Lg Lg;
save Eg Eg;

clear;
pack;
% =========================================================================

% =========================================================================
%  PROCESS SVD 1st STAGE FOR BLUE
% -------------------------------------------------------------------------
load status;

Lb = [];
Eb = [];

load Ib;
Eb(p*m, myRank) = 0;
Lb(myRank*m, nfc) = 0;

i2 = 0;
% for each block in the image
for f=1:myBlocksH;
    for g=1:myBlocksW;
        myMessage = sprintf('BLUE: Processing SVD for blocks [%.2u,%.2u]',f-1,g-1)
        i2 = i2 + 1;
        
        Ijb = double(reshape(Ib(i2,:,:),p,nfc)) / 255.0;

        [Ejb, Sjb, Cjb] = svd(Ijb,0);
        clear Ijb;
        Ejb = Ejb(1:p, 1:myRank);
        Sjb = Sjb(1:myRank, 1:myRank);
        Cjb = Cjb(1:nfc, 1:myRank);
        Ljb = Sjb * Cjb';
        clear Cjb;
        clear Sjb;
 
        Eb( (i2-1)*p + 1 : (i2-1)*p + p, 1:myRank) = Ejb(:,:);
        clear Ejb;

        Lb( (i2-1)*myRank + 1 : (i2-1)*myRank + myRank, 1:nfc ) = Ljb(:,:);
    end;
    pack;
end;

save Lb Lb;
save Eb Eb;

clear;
pack;
% =========================================================================


% =========================================================================
%  PROCESS SVD 2nd STAGE FOR RED
% -------------------------------------------------------------------------
load status;
myMessage = sprintf('Processing 2nd Stage RED SVD')
load Lr;
[Ur, Vlr, Vrr] = svd(Lr,0);
Ur = Ur(1:m*b, 1:q);
Vlr = Vlr(1:q, 1:q);
Vrr = Vrr(1:nfc, 1:q);
Vr = Vlr * Vrr';
save Ur Ur;
save Vr Vr;
clear;
pack;
% =========================================================================


% =========================================================================
%  PROCESS SVD 2nd STAGE FOR GREEN
% -------------------------------------------------------------------------
load status;
myMessage = sprintf('Processing 2nd Stage GREEN SVD')
load Lg;
[Ug, Vlg, Vrg] = svd(Lg,0);
Ug = Ug(1:m*b, 1:q);
Vlg = Vlg(1:q, 1:q);
Vrg = Vrg(1:nfc, 1:q);
Vg = Vlg * Vrg';
save Ug Ug;
save Vg Vg;
clear;
pack;
% =========================================================================

% =========================================================================
%  PROCESS SVD 2nd STAGE FOR BLUE
% -------------------------------------------------------------------------
load status;
myMessage = sprintf('Processing 2nd Stage BLUE SVD')
load Lb;
[Ub, Vlb, Vrb] = svd(Lb,0);
Ub = Ub(1:m*b, 1:q);
Vlb = Vlb(1:q, 1:q);
Vrb = Vrb(1:nfc, 1:q);
Vb = Vlb * Vrb';
save Ub Ub;
save Vb Vb;
clear;
pack;
% =========================================================================

% =========================================================================
%  LOAD FINAL RESULTS AND SAVE MATRICES IN TXT FILES
% -------------------------------------------------------------------------
myMessage = sprintf('Saving Final Results')
load Er;
save Er.txt Er -ASCII;
clear;

load Eg;
save Eg.txt Eg -ASCII;
clear;

load Eb;
save Eb.txt Eb -ASCII;
clear;

load Ur;
save Ur.txt Ur -ASCII;
clear;

load Ug;
save Ug.txt Ug -ASCII;
clear;

load Ub;
save Ub.txt Ub -ASCII;
clear;

load Vr;
save Vr.txt Vr -ASCII;
clear;

load Vg;
save Vg.txt Vg -ASCII;
clear;

load Vb;
save Vb.txt Vb -ASCII;
clear;
% =========================================================================

load status;
myClockFinish = fix(clock)
myProcessingTime = myClockFinish - myClockStart
save status;
myMessage = sprintf('FINISHED!')