% Implementation of Light Field Image Processing (Fourier Domain)

%
% Initialization of MATLAB (Turn Off in Time Measurement)
close all;  % Delete all figures
clear all;  % Delete all parameters
clc;        % Clear a command window
%

% Part 1: Config the Parameters

% The physical mask placed close to the sensor has 4 harmonics,
% therefore we will have 9 angular samples in the light field.

% imread: Input JPEG & PNG Image,   readpfm_color: Input Color PFM Image
% nAnlges : The Number of Angles (Odd Number),   cAngles : The Center of Angles


% Parameter 1

% The floor Number of multiple view images
t = 4;   % 3x3 Images
% t = 12;  % 5x5 Images
% t = 24;  % 7x7 Images
% t = 40;  % 9x9 Images
% t = 60;  % 11x11 Images


% Parameter 2

% Read 2D Image
% Input Image has to have Odd Size
% PFM = Portable Float Map
disp('Reading Input Image...');

% Color Image "Toy and Book": mxn = 2101x2827, Angles = 11, F1Y = 191, F1X = 257
% phi1 = 120;  phi2 = 120;  bst = 2 or 4;
file_name_in = 'ToyAndBook_c.jpg'; img = double(imread(file_name_in)); num = 'FSP_TAB_C'; GRAY_SCALE = 0; nAngles = 11;

% Color Image "Toy and Woman": mxn = 1575x2133, Angles = 9, F1Y = 175, F1X = 237
% phi1 = 150;  phi2 = 200;  bst = 4;
% file_name_in = 'ToyAndWoman_c.jpg'; img = double(imread(file_name_in)); num = 'FSP_TAW_C'; GRAY_SCALE = 0; nAngles = 9;

% Color Image "Toys": mxn = 1925x2585, Angles = 11, F1Y = 175, F1X = 235
% phi1 = 150;  phi2 = 360;  bst = 4;
% file_name_in = 'Toys_c.jpg'; img = double(imread(file_name_in)); num = 'FSP_TOY_C'; GRAY_SCALE = 0; nAngles = 11;

% GrayScale Image "Toy and Book": mxn = 2101x2827, Angles = 11, F1Y = 191, F1X = 257
% phi1 = 120;  phi2 = 120;  bst = 2 or 4;
% file_name_in = 'ToyAndBook_g.jpg'; img = double(imread(file_name_in)); num = 'FSP_TAB_G'; GRAY_SCALE = 1; nAngles = 11;

% GrayScale Image "Toy and Woman": mxn = 1575x2133, Angles = 9, F1Y = 175, F1X = 237
% phi1 = 150;  phi2 = 200;  bst = 4;
% file_name_in = 'ToyAndWoman_g.jpg'; img = double(imread(file_name_in)); num = 'FSP_TAW_G'; GRAY_SCALE = 1; nAngles = 9;

% GrayScale Image "Toys": mxn = 1925x2585, Angles = 11, F1Y = 175, F1X = 235
% phi1 = 150;  phi2 = 360;  bst = 4;
% file_name_in = 'Toys_g.jpg'; img = double(imread(file_name_in)); num = 'FSP_TOY_G'; GRAY_SCALE = 1; nAngles = 11;

cAngles = ceil(nAngles/2);

% Set GRAY_SCALE to 0 for color images. Will need more memory.
if(GRAY_SCALE)
   img = img(:,:,1);  % take red channel only
end

% Make image odd size
% img = img(1:end-1,1:end-6,:);

% find size of input image - m:height, n:width, CH:color channel
[m,n,CH] = size(img);  % Grayscale: CH = 1, Color: CH = 3

% The fundamental frequency of the cosine in the mask in pixels
% Cosine Frequency in Pixels from Calibration Image (Y:width, X:height)
F1Y  = floor(m/nAngles);  F1X  = floor(n/nAngles);  % Size of Each View Image
F12Y = floor(F1Y/2);      F12X = floor(F1X/2);      % Half Size of Each View Image


% Parameter 3

% PhaseShift due to Mask In-Plane Translation wrt Sensor
% Reduce Image Noise
% phi1 = 0;    phi2 = 0;    % Plain
phi1 = 120;  phi2 = 120;  % "Toy and Book"
% phi1 = 150;  phi2 = 200;  % "Toy and Woman"
% phi1 = 150;  phi2 = 360;  % "Toys"

% Compute Spectral Tile Centers, Peak Strengths and Phase.
for i=1:nAngles
    for j=1:nAngles
        CentY(i,j) = ceil(m/2) + (i-cAngles)*F1Y;  % 88~1488  (1~1575)
        CentX(i,j) = ceil(n/2) + (j-cAngles)*F1X;  % 119~2015 (1~2133)
        % Mat(i,j) = exp(-sqrt(-1)*((phi1*pi/180)*(i-cAngles) + (phi2*pi/180)*(j-cAngles)));
    end
end

% Mask Interpolation Unused
Mat = ones(nAngles,nAngles); % Spectrum Component


% Parameter 4

% Boost DC Components of FFT of 2D Mask Function
% 20 is because we cannot have negative values in the mask
% So strenght of DC component is 20 times that of harmonics
% bst = 1;    % Plain
bst = 2;    % "Toy and Book"
% bst = 4;    % "Toy and Woman" & "Toys"
% bst = 10;   % Bad
% bst = 20;   % Too Bad
Mat(cAngles,cAngles) = Mat(cAngles,cAngles) * bst;


% Parameter 5

% Method of Mask Interpolation
% intp = 'linear';   intpm = 1;  % Linear Interpolation
intp = 'cubic';    intpm = 2;  % Cubic Interpolation
% intp = 'spline';   intpm = 3;  % Spline Interpolation   % NG
% intp = 'nearest';  intpm = 4;  % Nearest Interpolation  % Bad



% Part 2: Convert 2D Image to 4D Light-Field

% Beginning of 4D light field computation
% do for all color channel

% Spectrum Image
imgs = zeros(m,n,CH);

for ch = 1:CH

    disp('=================================');
    disp(sprintf('Processing channel %d',ch));

    % Find FFT of image
    disp('Computing FFT of 2D image');
    f = fftshift(fft2(img(:,:,ch))); % f: row = m, column = n
    
    % Spectrum Image (phi1=0,phi2=0,bst=1)
    fs = log10(abs(f));
    imgs(:,:,ch) = fs;
        
    % Rearrange Tiles of 2D FFT into 4D Planes to obtain FFT of 4D Light-Field
    disp('Rearranging 2D FFT into 4D');
    
    for i = 1:nAngles
        for j = 1:nAngles
            idx1 = CentY(i,j)-F12Y:CentY(i,j)+F12Y;
            if (size(idx1,2)>F1Y)
                idx1(1)=[];
            end
            idx2 = CentX(i,j)-F12X:CentX(i,j)+F12X;
            if (size(idx2,2)>F1X)
                idx2(1)=[];
            end
            FFT_LF(:,:,i,j) = f(idx1,idx2) / Mat(i,j);
        end
    end

    clear f

    %
    % Mask Interpolation
    k = sqrt(-1); % imaginary number
    for i = 1:nAngles
        for j = 1:nAngles
            shift = (phi1*pi/180)*(i-cAngles) + (phi2*pi/180)*(j-cAngles);
            FFT_LF(:,:,i,j) = FFT_LF(:,:,i,j) * exp(k*shift);
        end
    end
    %

    disp('Computing inverse 4D FFT');
    LF = ifftn(ifftshift(FFT_LF)); % Compute Light-Field by 4D Inverse FFT

    clear FFT_LF

    if(ch==1)
        LF_R = LF;
    elseif(ch==2)
        LF_G = LF;
    elseif(ch==3)
        LF_B = LF;
    end

    clear LF

end

clear img

% Now we have 4D light field
disp('Light Field Computed. Done...');
disp('==========================================');

%{
% Write the Light-Filed Views
disp(sprintf('Processing channel %d',ch));
disp('Computing Light-Field View images');
disp('==========================================');

LF_count = 0;
for i = 1:nAngles
    for j = 1:nAngles
        if(GRAY_SCALE)
            LF_img = zeros(F1Y,F1X);
            LF_img(:,:) = LF_R(:,:,i,j);
        else
            LF_img = zeros(F1Y,F1X,3);
            for ch = 1:CH
                if(ch==1)
                    LF_img(:,:,1) = LF_R(:,:,i,j);
                elseif(ch==2)
                    LF_img(:,:,2) = LF_G(:,:,i,j);
                elseif(ch==3)
                    LF_img(:,:,3) = LF_B(:,:,i,j);
                end
            end
        end
        LF_count = LF_count + 1;
        LF_str = sprintf('%s_LF_DCx%02d_P1=%03d_P2=%03d_%03d.png',num,bst,phi1,phi2,LF_count);
        LF_Write(LF_img,LF_str);
    end
end
disp('==========================================');
%}



% Part 3: Convert 4D Light-Field to 2D Image

% Digital Refocusing Code
% Take a 2D slice of 4D light field
% For refocusing, we only need the FFT of light field, not the light field

disp('Synthesizing Refocused Images by taking 2D slice of 4D Light Field');

if(GRAY_SCALE)
    FFT_LF_R = fftshift(fftn(LF_R));
    clear LF_R
else
    FFT_LF_R = fftshift(fftn(LF_R));
    clear LF_R

    FFT_LF_G = fftshift(fftn(LF_G));
    clear LF_G

    FFT_LF_B = fftshift(fftn(LF_B));
    clear LF_B
end

% Height and Width of refocused image
% Caution: [H, W] = size(FFT_LF_R)  means  H = F1Y, W = F1X * nAngles * nAngles
H = size(FFT_LF_R,1);  W = size(FFT_LF_R,2);

count = 0;
for theta = -t:t

    count = count + 1;
    
    % Name of Output Image
    str = sprintf('%s_DCx%02d_P1=%03d_P2=%03d_%s_%03d.png',num,bst,phi1,phi2,intp,count);
    
    disp('===============================================');
    disp(sprintf('Calculating New ReFocused Image: theta = %d',theta));
    
    if(GRAY_SCALE)
        img_out = RefocusFSP(FFT_LF_R,[theta,theta],intpm);
    else
        img_out = zeros(H,W,3);
        img_out(:,:,1) = RefocusFSP(FFT_LF_R,[theta,theta],intpm);
        img_out(:,:,2) = RefocusFSP(FFT_LF_G,[theta,theta],intpm);
        img_out(:,:,3) = RefocusFSP(FFT_LF_B,[theta,theta],intpm);
    end
    
    % Scale RefocusedImage in [0,255]
    img_out = img_out - min(img_out(:)); % Minimum Number is 0
    img_out = img_out/max(img_out(:));   % Maximum Number is 1
    img_out = uint8(255*img_out);        % values are 0 ~ 255
    imwrite(img_out,str);
    
    % Histgram(str);
    
    disp(sprintf('Refocused image written as %s',str));
    
    % figure(1);  image(img_out);         % Display Image
    % figure(1);  imagesc(img_out);       % Display Image
    % imtool(str);                       % imtool = imshow
    
    %{
    % Image Filtering (Low-Pass)
    if(GRAY_SCALE)
        [img_out_F,strf] = LPF_FSP(FFT_LF_R,[theta,theta],str);
    else
        img_out_F = zeros(H,W,3);
        [img_out_F(:,:,1),strf] = LPF_FSP(FFT_LF_R,[theta,theta],str);
        [img_out_F(:,:,2),strf] = LPF_FSP(FFT_LF_G,[theta,theta],str);
        [img_out_F(:,:,3),strf] = LPF_FSP(FFT_LF_B,[theta,theta],str);
    end
    
    % Scale FilteredImage in [0,255]
    img_out_F = img_out_F - min(img_out_F(:)); % Minimum Number is 0
    img_out_F = img_out_F/max(img_out_F(:));   % Maximum Number is 1
    img_out_F = uint8(255*img_out_F);          % values are 0 ~ 255
    imwrite(img_out_F,strf);                   % Image Filtering
    
    % Histgram(strf);
    
    disp(sprintf('Filtered image written as %s',strf));
    
    % figure(2); image(img_out_F);         % Display Image
    % figure(2); imagesc(img_out_F);       % Display Image
    % imtool(strf);                        % imtool = imshow
    %}
    
end

%
% Spectrum Color Image
strs = sprintf('Spectrum_%s_DCx%02d_P1=%03d_P2=%03d.png',num,bst,phi1,phi2);
imgs = imgs - min(imgs(:));  % Minimum Number is 0
imgs = imgs/max(imgs(:));    % Maximum Number is 1
imgs = uint8(255*imgs);      % values are 0 ~ 255
imwrite(imgs,strs);

disp('==========================================');
disp(sprintf('Spectrum image written as %s',strs));
%

disp('==========================================');
