function [filtvid] = okadaFilterVideo(video)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Performs an 2D okada filter described in: https://doi.org/10.1371/journal.pone.0157595
% Inputs:
%     - video: the input video to be filtered. Can also be an image
% Output:
%     - filtvid: the filtered video
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[d1 d2 d3] = size(video);

% Reshape the video to improve the speed
linvideo = double(reshape(video, [d1*d2 d3]));
alpha = 1e8;

linvideoF(:, 1) = linvideo(:, 1);
indblock = 1:d1*d2;
indblock = reshape(indblock, [d1 d2]);

% Circularly shift the indexes so that it's easier to specify the values desired for the
% 2D okada filter in the for loop. Then convert back to a vector.
indu = reshape(circshift(indblock, -1, 2), [d1*d2 1]);
indd = reshape(circshift(indblock, 1, 2), [d1*d2 1]);
indr = reshape(circshift(indblock, 1, 1), [d1*d2 1]);
indl = reshape(circshift(indblock, -1, 1), [d1*d2 1]);

% Calculate the filtered values
for fn = 1:d3
    frame = linvideo(:, fn);
    T1n = frame(indu)+frame(indd)-2*frame;
    T1d = 2*(1+exp(-alpha.*(frame-frame(indu)).*(frame-frame(indd))));
    T2n = frame(indr)+frame(indl)-2*frame;
    T2d = 2*(1+exp(-alpha.*(frame-frame(indr)).*(frame-frame(indl))));
    linvideoF(:, fn) = frame+(T1n./T1d+T2n./T2d)./2;
end

% Reshape to the original dimensions
filtvid = reshape(linvideoF, [d1 d2 size(linvideoF, 2)]);
end