SpatialStatsFFT.m

function [T xx] = SpatialStatsFFT( A1,A2, varargin)
% Computes that spatial statistics of input signals A1 and A2.  The
% functions appropriately normalizes the statistics.  This function is
% optimized using fast Fourier transforms to expedite the computation of
% the convolution.
%
% Nonperiodic and periodic refer the the boundary conditions under which
% the material information was generated.  Experimental information will
% always have nonperiodic boundary conditions while some simulations with
% have periodic boundary conditions.  As a result, this function defaults
% to computing nonperiodic statistcs.
%
% This function can also compute statistics for partial datasets that are
% sampled on an evenly spaced pixel/voxel grid.  If some information is
% known and some is not, then one can compute the spatial statistics using
% this function.  A tutorial will be made available later.
%
% [T, xx] = f2( DATA1 ) computes the nonperiodic autocorrelation of DATA1.
% T is an array of the spatial statistics corresponding to the cell arrray
% of vector lengths in xx.
%
% [T, xx] = f2( DATA1, DATA2 )  computes the nonperiodic crosscorrelation
% of DATA1 and DATA2 where they are the head and tail of the vector, respectively.
%
% [T, xx] = f2( DATA1, DATA2, 'ARG1', val1, ..,.., 'ARGN', valN ) computes
% the crosscorrelation of DATA1 and DATA2 with a set of optional parameters
% listed below.
%
% [T, xx] = f2( DATA1, [], 'ARG1', val1, ..,.., 'ARGN', valN ) computes the
% autocorrelation of DATA1 with optional parameters
%
% Parameters ARGUMENTS
% --------------------
% ARG - class - default - description
% -----------------------------------
% 'normalize' - logical - true - Spatial statistics require computing the
%   convolutions in the numerator and denominator separately.  The numerator
%   is the number of times the statistics criteria is satisfied while the
%    denominator is the number events that were sampled.  If
%    f2( DATA1, [], 'normalize',false ) is executed then only the top
%   numerator is returned.  It is really just the convolution of DATA1.
%
% 'display' - logical - true - When display is true, a plot will appear
%    after each function call of the spatial statistics.
%
% 'cutoff' - double - Cutoff is the maximum vector size to be returned in
%   the statistics.  f2( DATA1, [], 'cutoff', 10 ) will return spatial
%   statistics for vectors whose elements are less than or equal to 10;
%
% 'periodic' - [1xd] logical - If the model information was generated by
%   simulation then it may have periodic boundary conditions.
%   f2( DATA1, [], 'periodic', true ) forces all boundaries to be periodic
%   f2( DATA1, [], 'periodic', [ true true false] ) places the condition
%   that the information is periodic in dimensions 1 and 2 and nonperiodic
%   in the last.  Conditions like this are useful for interfatial
%   simulation information.
%
% 'Mask1' - [N1xN2xN3] double - Sometimes the data returned is not complete,
%   datapoints may be unpopulated. [N1xN2xN3] is the size of the input information.
%   f2( DATA1, [], 'Mask1', M )  computes are spatial statistics for
%   populated datapoints and normalizes the function appropriately.  M is
%   a logical array that describes whether a value is populated<true> or
%   not populated<false> for the first dataset, the one at the tail of the vector .
%
% 'Mask2' - [N1xN2xN3] double -
%   f2( DATA1, [], 'Mask2', M )  computes are spatial statistics for
%   populated datapoints and normalizes the function appropriately.  M is
%   a logical array that describes whether a value is populated<true> or
%   not populated<false> for the second dataset, the one at the head of the vector .
%
% 'Mask' - [N1xN2xN3] logical -populates Mask1 and Mask2 with the same
% mask.
%
%  'shift' - logical - true - if true performs the fftshift.
%
% if varargin is a structure with all the parameter fields then the
% parameters are automatically defined and not decided.  This is useful in
% repetative applications.

%% Quickstart :: Spatial Statistics from web url
% * Download an image from the web
% * Extract States from the Image
% * Store the image as a workspace variable
% * Compute Spatial Statistics on the data.


if ischar( A1 );
    % Find potential variable names
    varnm = evalin('base', 'genvarname(''SpatialStatsVar'',who);');
    
    % Create a structure in the path that has the image data and the url
    assignin( 'base', varnm, struct( 'data', cast(  ...
                                                imread( A1 ), ... Image from the web
                                                'double' ), ...
                                      'url', A1 ) );
    
    % Save a local variable for the spatial statistics
    A1 = round( evalin( 'base', sprintf( '%s.data',varnm) )./255 );

    
    % Alert user of the new workspace variable
    disp( ...
        sprintf( 'The image data has been saved to your workspace as ``%s``.', varnm ) ...
        );
    
    % Save the variable in the workspace
    param.auto = true; 
end


%% Initialize Critical Data Elements
if numel(varargin) == 1 && isstruct( varargin )
    param = varargin;
else
    param = setparam(varargin, size(A1));
end


%%
% Decide whether the correlation is an auto or cross correlation

if exist('A2','var') && numel(A2) > 0
    if all(size(A1) == size(A2))
        if any(A1(:)~=A2(:)) param.auto = false; end
    else error('The size of the input signals are not the same.'); end
end

%% Compute numerator

if param.auto  % Autocorrelation
    if numel( param.Mask1 ) == 0
        T = convolveSSFFT( param.periodic, A1 );
    else   % Partial
        T = convolveSSFFT( param.periodic, param.Mask1.*A1 );
    end
else   % Crosscorrelation
    if numel( param.Mask1 ) == 0
        T = convolveSSFFT( param.periodic, A1, A2 );
    else  % Partial
        T = convolveSSFFT( param.periodic, param.Mask1.*A1, param.Mask2.*A2 );
    end
end


%% Compute denominator
if param.normalize
    if numel(param.Mask1)  == 0  % Complete
        if all(param.periodic)  % All Periodic
            T = T./ numel( A1 );
        else  % Partial and Nonperiodic
            % compute it directly, for the meantime use the convolution
            T(:) = T./ convolveSSFFT( param.periodic, ones(size(A1)));
        end
    else  % Partial
        if param.auto  % Auto
            T(:) = T./ convolveSSFFT( param.periodic, param.Mask1 );
        else   % Cross
            T(:) = T./ convolveSSFFT( param.periodic, param.Mask1, param.Mask2 );
        end
    end
end

%% T vector sizes
for ii = 1 : ndims(T)
    uu = 1 : floor(size(T,ii)./2);
    
    if mod( size(T,ii), 2) == 0
        xx{ii} = [uu-1,fliplr(uu)*-1];
    else
        xx{ii} = [ 0 , uu, fliplr( uu ) * -1 ];
    end
    
    incut{ii} = abs(xx{ii}) > param.cutoff(ii);
end

%% Truncate Statistics
T(incut{1},:,:) = [];
if ndims(T) >= 2
    T(:,incut{2},:) = [];
end
if ndims(T) >= 3
    T(:,:,incut{3}) = [];
end

if param.shift
    xx = arrayfun( @(x)fftshift(xx{x}(~incut{x})),1:ndims(T),'UniformOutput',false);
    T(:) = fftshift(T);
else
    xx = arrayfun( @(x)xx{x}(~incut{x}),1:ndims(T),'UniformOutput',false);
end

if param.vector 
    if nargout == 2
        switch ndims(T)
            case 1
                xx = xx{1}(:);
            case 2
                [X1,X2] = meshgrid( xx{1},xx{2});
                xx = [X1(:),X2(:)];
            case 3
                [X1,X2,X3] = meshgrid( xx{1},xx{2},xx{3});
                xx = [X1(:),X2(:),X3(:)];
        end
    end
    T = T(:);
end
%% Display the statistics
% The shifted image is the more canonical form to view the statistics in.
% In fact surface is the optimum viewing tools because then you can use
% the datatip to find probabilities for certain vectors.
%

if param.display
    % When the statistics are visualized, the outputs are
    % forced to be real, this result should be removed
    if ndims( A1 )== 2 & ~any( size(A1) == 1);
            if param.shift
                pcolor(xx{2},xx{1},real(T)); 
            else
                pcolor(fftshift(xx{2}),fftshift(xx{1}),fftshift(real(T))); 
            end
            xlabel('t_x','Fontsize',16); ylabel('t_y','Fontsize',16, 'Rotation',0); 
            hc = colorbar; shading flat; axis equal
            if param.normalize
                str = 'Probability density';
            else
                str = 'Counts';
            end
            set( get( hc, 'Ylabel'), 'String', str, 'Fontsize',16,'Rotation',270,'VerticalAlignment','Bottom');
    elseif ndims( A1 )== 2 & any( size(A1) == 1);
            if param.shift
                plot(xx{getOutput(@max,2,cellfun(@(x)numel(x),xx))},real(T),'Linewidth',3,'Color','k'); 
            else
                plot(fftshift(xx{getOutput(@max,2,cellfun(@(x)numel(x),xx))}),fftshift(real(T)),'Linewidth',3,'Color','k'); 
            end
            xlabel('t_x','Fontsize',16); ylabel('Probability','Fontsize',16); 
            ylim([0 1]); xlim([0 ceil(numel(T)./2)]); grid on
    else
        disp('Visualization is not available for 3-D statistics yet.')
    end
end


end


%%
function param = setparam( inptvar,sz )
% Set parameters for the code

param = struct('normalize',true, ...
    'display',true,...
    'cutoff', sz./2,...
    'auto',true,...
    'periodic',false*ones(1,numel(sz)),...
    'Mask1',[],...
    'Mask2',[], ...
    'shift', false,...
    'vector',false,...
    'integrate',false);


fldnm = fieldnames( param );

if numel(inptvar) > 0
    for ii = [1 : (numel(inptvar)./2)]*2-1
        if ismember( inptvar{ii}, fldnm );
            % update parameters
            param = setfield( param, inptvar{ii}, inptvar{ii+1});
        elseif strcmp( inptvar{ii}, 'Mask' );
            [ param.Mask1 param.Mask2 ] = deal( inptvar{ii+1} );
        else % error message for bad input parameters
            disp(sprintf('f2 accepts the following options:'));
            for jj = 1 : numel( fldnm )
                disp(sprintf(':: %s - type ::  %s', fldnm{jj}, class(getfield(param, fldnm{jj}))));
            end
            disp(sprintf(':: %s - type ::  %s', 'Mask', class(getfield(param, 'Mask1'))));
            error( sprintf('%s is not a valid parameter.', inptvar{ii}));
        end
    end
    
    n = {size( getfield( param, 'Mask1' )) ;  size( getfield( param, 'Mask2' ))};
    if ~( all(n{1}==0) & all(n{2}==0)) && (numel( n{1})~=numel(n{2}) || ~all( n{1}==n{2}))
        disp('test')
        if numel( getfield( param, 'Mask1') ) == 0 param.Mask1 = ones( size(getfield( param, 'Mask2')));
        elseif numel( getfield( param, 'Mask2') ) == 0 param.Mask2 = ones( size(getfield( param, 'Mask1')));
        else error( 'The size of Mask1 and Mask2 are not the same.'); end
        
    end
end
param.Mask1 = cast( param.Mask1, 'double' );
param.Mask2 = cast( param.Mask2, 'double' );

if numel( param.cutoff ) == 1 param.cutoff = param.cutoff * ones(1, numel( sz )); end
id = find(isinf(param.cutoff)); if numel(id) > 0 param.cutoff(id) = sz(id)./2; end
end


%%

%% Padding procedure for nonperiodic structures

function fA = FourierPadSSFFT( A, osz, nsz );

if all( osz == nsz )
    fA = double(A);
else
    switch numel(nsz)
        case 1
        case 2
            
            fA = [ A, zeros([osz(1),nsz(2)-osz(2)]); zeros( nsz(1)-osz(1),nsz(2))];
        case 3
            fA = cat( 3, [A, zeros( [osz(1) nsz(2)-osz(2) osz(3)]); zeros( [nsz(1)-osz(1) osz(2) osz(3)]), zeros( [nsz(1)-osz(1) nsz(2)-osz(2), osz(3)])],...  
                zeros( [ nsz(1) nsz(2), nsz(3) - osz(3)]));
    end
end

for ii = 1 :  numel(nsz)
    fA(:) = fft( fA, [], ii);
end
end


%%  Grid Point Cloud Data

function [ V xx ] = Point2Grid( Data, Lims, Out_SZ );

if ~exist( 'Lims', 'var' ) || numel( Lims ) == 0 
    Lims = [min(Data)',max(Data)'];
end

SUBS = bsxfun( @minus, Data, Lims(:,1)');
SUBS(:) = bsxfun( @rdivide, SUBS, diff(Lims,[],2)');
SUBS(:) = bsxfun( @plus, SUBS, 1./Out_SZ(:)'./2 );
SUBS(:) = round(bsxfun( @times, SUBS, Out_SZ(:)'));

%%
if numel( Out_SZ) == 1
    V = zeros(  [Out_SZ 1]);
else
    V = zeros(  Out_SZ );
end

switch numel(Out_SZ)
    case 1
        V(:) = accumarray( SUBS(:,1), ...
            ones( size(SUBS(:,1))), ...
            [ prod(Out_SZ),1], @sum);
    case 2
        V(:) = accumarray( sub2ind( Out_SZ, SUBS(:,1), SUBS(:,2)), ...
            ones( size(SUBS(:,1))), ...
            [ prod(Out_SZ),1], @sum);
    case 3
        
        V(:) = accumarray( sub2ind( Out_SZ, SUBS(:,1), SUBS(:,2), SUBS(:,3)), ...
            ones(size(SUBS(:,1))), ...
            [ prod(Out_SZ),1], @sum);
end

if nargout == 2
    for ii = 1 : numel( Out_SZ );
        dx = diff(Lims(ii,:))./(Out_SZ);
        xx{ii} = (Lims(ii,1) + dx/2):dx:(Lims(ii,2));
    end
end

end

%%

function fA = convolveSSFFT(period,A1,A2)
% convolve pads the data signals then performs the convolution of signals
% A1 and A2 using fast fourier transform algorithms.
% This is the degree to which the image is padded because we are only
% extracting half of the vector
period_mult= .6; %describe this value

% in reality you only need to pad to rmax
realvals = isreal(A1)*ones(1,2);  % supresses imaginary parts if the values in the statistics
nsz = ceil([ones(size(period)) + (1-period)*period_mult].*size( A1));


fA = FourierPadSSFFT( A1, size( A1 ), nsz );

if exist('A2','var') && numel( A2 ) > 0
    if all( nsz == ceil([ones(size(period)) + (1-period)*period_mult].*size( A2)) );    % THIS CONDITION HERE IS PROBLABLY A PROBLEM
        realvals(2) = isreal(A2);
        fA(:) = fA .* conj( FourierPadSSFFT( A2, size(A1), nsz ));
    else
        error('The sizes are incaptable.')
    end
else % Weiner-Kinchin Method to Autocorrelations
    fA(:) = abs( fA ).^2;
end

fA(:) = ifftn( fA );

if all(realvals) fA(:) = real(fA); end;


end
%%
function varargout = getOutput(func,outputNo,varargin)
% A function that selectively chooses the output of another function.
% func - function handle
% Outputnp - is the output index of interest
% varargin - the inputs to the function
    varargout = cell(max(outputNo),1);
    [varargout{:}] = func(varargin{:});
    varargout = varargout(outputNo);
end