% PERFFT2  2D Fourier transform of Moisan's periodic image component
%
% Usage: [P, S, p, s] = perfft2(im)
%
% Argument:  im - Image to be transformed
% Returns:    P - 2D fft of periodic image component
%             S - 2D fft of smooth component
%             p - Periodic component (spatial domain)
%             s - Smooth component (spatial domain)
%
% Moisan's "Periodic plus Smooth Image Decomposition" decomposes an image 
% into two components
%        im = p + s
% where s is the 'smooth' component with mean 0 and p is the 'periodic'
% component which has no sharp discontinuities when one moves cyclically across
% the image boundaries.  
%
% This wonderful decomposition is very useful when one wants to obtain an FFT of
% an image with minimal artifacts introduced from the boundary discontinuities.
% The image p gathers most of the image information but avoids periodization
% artifacts.
%
% The typical use of this function is to obtain a 'periodic only' fft of an
% image 
%   >>  P = perfft2(im);
%
% Displaying the amplitude spectrum of P will yield a clean spectrum without the
% typical vertical-horizontal 'cross' arising from the image boundaries that you
% would normally see.
%
% Note if you are using the function to perform filtering in the frequency
% domain you may want to retain s (the smooth component in the spatial domain)
% and add it back to the filtered result at the end.  
%
% The computational cost of obtaining the 'periodic only' FFT involves taking an
% additional FFT.
%
%
% Reference: 
% This code is adapted from Lionel Moisan's Scilab function 'perdecomp.sci' 
% "Periodic plus Smooth Image Decomposition" 07/2012 available at
%
%   http://www.mi.parisdescartes.fr/~moisan/p+s
%
% Paper:
% L. Moisan, "Periodic plus Smooth Image Decomposition", Journal of
% Mathematical Imaging and Vision, vol 39:2, pp. 161-179, 2011.

% Peter Kovesi
% www.peterkovesi.com
% September 2012

function [P, S, p, s] = perfft2(im)
    
    if ~isa(im, 'double'), im = double(im); end
    [rows,cols] = size(im);
    
    % Compute the boundary image which is equal to the image discontinuity
    % values across the boundaries at the edges and is 0 elsewhere
    s = zeros(size(im));
    s(1,:)   = im(1,:) - im(end,:);
    s(end,:) = -s(1,:);
    s(:,1)   = s(:,1)   + im(:,1) - im(:,end);
    s(:,end) = s(:,end) - im(:,1) + im(:,end);
    
    % Generate grid upon which to compute the filter for the boundary image in
    % the frequency domain.  Note that cos() is cyclic hence the grid values can
    % range from 0 .. 2*pi rather than 0 .. pi and then pi .. 0
    [cx, cy] = meshgrid(2*pi*[0:cols-1]/cols, 2*pi*[0:rows-1]/rows);    
    
    % Generate FFT of smooth component
    S = fft2(s)./(2*(2 - cos(cx) - cos(cy)));
    
    % The (1,1) element of the filter will be 0 so S(1,1) may be Inf or NaN
    S(1,1) = 0;          % Enforce 0 mean 

    P = fft2(im) - S;    % FFT of periodic component

    if nargout > 2       % Generate spatial domain results 
        s = real(ifft2(S)); 
        p = im - s;         
    end
    
