% Refocus in Fourier Domain

function IOut = RefocusFSP(FFTLF,theta,intpm)

% FFTLF has (nAngles x nAngles) View Images
% m = F1Y, n = F1X
m = size(FFTLF,1);  n = size(FFTLF,2);

% theta is slice coordinate vector [theta, theta]
% theta > 0: focus to short, theta < 0: focus to back
Theta1 = theta(1);   Theta2 = theta(2);

% The Center of Each View Image
% cTem = [F12Y F12X cAngles cAngles]
% cTem(1) = F12Y, cTem(2) = F12X, cTem(3) = cTem(4) = cAngles
cTem = floor(size(FFTLF)/2) + 1;


% Find the coordinates of 2D slice

% XX: (F1Y x F1X),  YY: (F1Y x F1X)
% XX = [1 ... n] x m, YY = [1 ... m]' x n (': transposed vector)
[XX,YY] = meshgrid(1:n,1:m);

cc = (XX - cTem(2)) / n;    dd = (YY - cTem(1)) / m;    % Fraction from Center
% cTem(2) = F12X x ones(F1Y,F1X),  cTem(1) = F12Y x ones(F1Y,F1X)
% cc: (F1Y x F1X),  dd: (F1Y x F1X)
% cc = [-0.5 ~ 0 ~ 0.5] x m,  dd = [-0.5 ~ 0 ~ 0.5]' x n (': transposed vector)
cc = Theta2*cc + cTem(4);   dd = Theta1*dd + cTem(3);   % Coordinates of 2D Slice

% Original Focus: cc = cTem(4) * ones(m,n),  dd = cTem(3) * ones(m,n)


% Resample 4D light field along the 2D slice
% Function "interpn" extracts 2D Slice from FFTLF
% FFTLF: 4D (F1Y x F1X x nAngles x nAngles),  v: 2D (F1Y x F1X)
if (intpm==1)
    v = interpn(FFTLF,YY,XX,dd,cc,'linear');   % Bilinear Interpolation
elseif (intpm==2)
    v = interpn(FFTLF,YY,XX,dd,cc,'cubic');    % Bicubic Interpolation
elseif (intpm==3)
    v = interpn(FFTLF,YY,XX,dd,cc,'spline');   % Spline Interpolation   % NG 
elseif (intpm==4)
    v = interpn(FFTLF,YY,XX,dd,cc,'nearest');  % Nearest Interpolation  % Bad
end


% 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;


% k = v(ceil(m/2),ceil(n/2));  k = floor(k/20);


% Take inverse 2D FFT to get the image
IOut = real(ifft2(ifftshift(v)));
