function [img_cleaned]=striperemoval(img, varargin)
%STRIPEREMOVAL  removes interference fringes/stripes from an image.
%  [img_cleaned]=STRIPEREMOVAL(img, varargin)
%     img - image with stripes to be removed
%     varargin - not used, it is here for compatibility
% By Eugeniy E. Mikhailov 2014/11/25
% 
% Be aware of "slow" spatial stripes 
% i.e. stripes with characteristic width comparable to the Gaussian beam.
% Their Fourier components will overlap with Gaussian.

% Fringes appear like strong isolated peaks in the Fourier transformed image.
% What is even nicer they appear symmetrically with respect to origin but not symmetrical to 
% the fold with respect to X or Y axes. See figure below
%
%                 ^
%                 |
%               * |
%                 |
%                oOo
%   -------------O0O-------------->
%                oOo <- Gaussian beam fft
%                 |
%                 | * <- stripes fft peak
%                 |
%
%
% So it easy to clean them, especially if their spatial frequency is high.
%

img_fourier=fftshift(fft2(img)); % move to Fourier space

imgfa= abs(img_fourier); % amplitude of Fourier components
imgfa= imgfa/max(imgfa(:)); % normalized amplitudes

[Ny, Nx] = size (img);     % image size

mask = img_fourier.*0;
for xi = 1:Nx
    % strictly speaking one of dimensions can be only half sampled due to symmetry of the mask
    for yi = 1:Ny
        xr = Nx - xi + 1; % reflected x
        yr = Ny - yi + 1; % reflected y
        mt = imgfa( yi, xi ); % test point
        ms = imgfa( yr, xr ); % symmetrical point with respect to origin 
        mr = imgfa( yi, xr ); % symmetrical point with respect to Y axis (right)
        md = imgfa( yr, xi ); % symmetrical point with respect to Y axis (down)

        m = ( mt + ms ) - ( mr + md ); % high value to the points shown at figure
        m = m/max([mt, ms, mr, md]);
        mask(yi,xi) = m;
    end
end

threshold = max(mask(:))/2;
mask = mask < threshold; % points due to stripes need to be removed

img_cleaned = abs( ifft2(ifftshift(img_fourier.*mask)) );