% Refocus in Fourier Domain
% Low-Pass Filter

function [IOut1,IOut2] = LPF_FSP(FFTLF,theta,str)

[m,n,p,q] = size(FFTLF);    % m = F1Y, n = F1X, p = q = nAngles
Theta1 = theta(1);
Theta2 = theta(2);

cTem = floor(size(FFTLF)/2) + 1;


% Find the coordinates of 2D slice
[XX,YY] = meshgrid(1:n,1:m);
cc = (XX - cTem(2))/size(FFTLF,2);
cc = Theta2*cc + cTem(4);
dd = (YY - cTem(1))/size(FFTLF,1);
dd = Theta1*dd + cTem(3);


% Resample 4D light field along the 2D slice
v = interpn(FFTLF,YY,XX,dd,cc,'cubic');


% Set NaN values to zero
idx = find(isnan(v)==1);
disp(sprintf('Number of Nans in sampling = %d',size(idx,1)));
v(isnan(v)) = 0;


% Image Filtering

% Filtering Matrix f1 (Low-Pass Filter, Etc)
% f1 = fspecial('average',[3 3]);       strf = sprintf('F_Average_%s',str);
% f1 = fspecial('gaussian',[3 3],0.5);  strf = sprintf('F_Gaussian_%s',str);
% f1 = [1,1,1,1,1;1,2,2,2,1;1,2,3,2,1;1,2,2,2,1;1,1,1,1,1] ./ 35;   strf = sprintf('F_Kiya_%s',str);
f1 = [-1,-1,-1;-1,9,-1;-1,-1,-1];  strf = sprintf('F_Sharp_%s',str);  % Sharpening Filter
[a,b] = size(f1);

%{
% High-Emphasis Filter  % Failure
strf = sprintf('E_%s',strf);
k = 1;
for u = 1:a
    for v = 1:b
        f1(u,v) = k + 1 - k * f1(u,v);
    end
end
%}

% Filtering Matrix f (mxn)
f = zeros(m,n);  f((m-a)/2+1:(m+a)/2,(n-b)/2+1:(n+b)/2) = f1;
% fvtool(f);
% wvtool(f);
% imhist(f);
F = fftshift(fft2(f));

% Filtering
v = F .* v; % Multiply in Fourier Domain


% Take inverse 2D FFT to get the image
IOut1 = real(ifft2(ifftshift(v)));
IOut1 = fftshift(IOut1);
IOut2 = strf;
