% PERFFT3  3D Fourier transform of Moisan's periodic image component
%
% Usage: [P, S, p, s] = perfft3(im)
%
% Argument:  im - Image to be transformed
% Returns:    P - 3D fft of periodic image component
%             S - 3D 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 = perfft3(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
% Centre for Exploration Targeting
% The University of Western Australia
% peter.kovesi at uwa edu au
% September 2012

%Joshua Morgan 2017
%Adapted to 3D filter
%All credit to original author

function [P, S, p, s] = perfft3(im)
    
    if ~isa(im, 'double'), im = double(im); end
    [rows,cols,slices] = 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,:);
    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, cz] = meshgrid(2*pi*[0:cols-1]/cols, 2*pi*[0:rows-1]/rows, 2*pi*[0:slices-1]/slices);    
    
    % Generate FFT of smooth component
    S = fftn(s)./(2*(3 - cos(cx) - cos(cy) - cos(cz))); %MAYBE CHANGE 2 TO 3?
    
    % The (1,1) element of the filter will be 0 so S(1,1) may be Inf or NaN
    S(1,1,1) = 0;          % Enforce 0 mean 

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

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