function [f_out, fn, pixel_shifts, Dval] = ring_remove(k_in, k1, k2, fract)
% Based on E. Kellner et al, Mag. Res. Med. vol. 76, pages 1574 to 1581,
% 2016.
% k_in = input vector.
% k1, k2 as defined in Kellner's paper.
% fract defines the range of sub pixels shift [-fract, fract].

nshifts = 50;       % total shifts = 2*nshifts.
% step 1: calculate nshifts shifted versions of k_in.
% =========================================================
if nargin < 4
    fract = 0.5;
end
fract = max(min(fract, 0.5), 0.05);
shift_pixels = (-nshifts:nshifts - 1)*fract/nshifts;

k_in = k_in(:);
N = length(k_in);
vec = -N/2:N/2 - 1;
% vec = 0:N - 1;       % equivalent.
cexp_vec = exp(-2*pi*1i/N*vec(:)*shift_pixels);

ks = k_in.*cexp_vec;
fs = ifftc(ks);

% step 2: find abs diff from the right (Dp) and left (Dn) of all points.
% ======================================================================
% See explanation at WORD document.
% =====================================
diffD = abs(diff(fs));

Dn = zeros(size(fs));
Dp = zeros(size(fs));

for j = 1:N
    if j <= N - k2
        Dp(j, :) = mean(diffD(j + k1 - 1:j + k2 - 1, :));
    end
    if j >= 1 + k2
        Dn(j, :) = mean(diffD(j - k2:j - k1, :));
    end
    Dp(N - k2 + 1:end, :) = inf;
    Dn(1:k2, :) = inf;
end

D = min([Dp(:), Dn(:)], [], 2);       % select min between Dp and Dn.
D = reshape(D, size(Dp));

% step 3: find shifted vector with min TV.
% ===========================================
[~, ind] = min(D, [], 2);       % find indices of min TV.
shift_vals = shift_pixels(ind) + (1:N);       % find pixel shifts.
pixel_shifts = shift_pixels(ind);

% calculate output vector with min TV for each pixel.
fn = zeros(N, 1);
Dval = zeros(N, 1);
for j = 1:N
    fn(j) = fs(j, ind(j));
    Dval(j) = D(j, ind(j));
end
Dval = D(:, nshifts + 1) - Dval;

% step 4: linear interpolation back to 1:N locations.
% ======================================================
F = griddedInterpolant(shift_vals.', fn);
f_out = F(1:N).';
