mfbox_sobimats.m

toolbox for spm 5 for data, model free analysis

function [ci,M,A]=mfbox_sobimats(varargin)
% measures the spatial decorrelation using multidimensional delays
%
% usage: [ci,M,A]=mfbox_sobimats(X[,Y],[maskcolor[,num[,step[,normalize]]]])
%
% Copyright by Peter Gruber
% Signal Processing & Information Theory group
% Institute of Biophysics, University of Regensburg, Germany
% Homepage: http://research.fabian.theis.name
%           http://www-aglang.uni-regensburg.de
%
% This file is free software, subject to the 
% GNU GENERAL PUBLIC LICENSE, see gpl.txt

usemask = false;
normalize = false;
num = 10;
step = 1;
maskcolor = NaN;

if (nargin<1), return; end

xmask = []; ymask = [];
if (~iscell(varargin{1}))
    X = varargin{1};
else
    X = varargin{1}{1};
    xmask = varargin{1}{2};
end
Y = [];
if (nargin>1)
    c = 1;
    if (iscell(varargin{2}))
        Y = varargin{2}{1};
        ymask = varargin{2}{2};
        c = 2;
    elseif ((length(xmask)==0)&&(length(varargin{2})==1))
        maskcolor = varargin{2}; c = 2;
    elseif (length(varargin{2})>1)
        Y = varargin{2};
        if (nargin>2)
            if ((length(xmask)==0)&&(length(varargin{3})==1))
                maskcolor = varargin{3}; c = 3;
            end
        end
    end
    if (nargin>c), if (length(varargin{c+1})==1), num = varargin{c+1}; c = c+1; end, end
    if (nargin>c), if (length(varargin{c+1})==1), step = varargin{c+1}; c = c+1; end, end
    if (nargin>c), if (length(varargin{c+1})==1), normalize = logical(varargin{c+1}); end, end
end

s = size(X);
if (~isempty(xmask))
    if (~isempty(ymask))
        if (length(size(xmask))==length(size(ymask)))
            if (any(size(xmask)~=size(ymask)))
                error('incompatible masks');
            end
        else
            error('incompatible masks');
        end
    end
    s = [size(xmask),s(end)];
end
X = reshape(X,[],s(end)); 

if ((~isempty(Y))&&(length(size(X))==length(size(Y))))
    t = size(Y);
    if (any(s(1:(end-1))~=t(1:(end-1)))), error('incompatible data verctors'); end
    Y = reshape(Y,[],t(end));
end

m = s(end);
s = s(s(1:(end-1))>1);
if (length(s)<2), rs = [s(1),1]; else rs = s; end
N = prod(rs);
s = [s,m];

if (~isempty(xmask))
    if (isempty(Y))
        mask = reshape(xmask,rs);
    elseif (~isempty(ymask))
        mask = reshape(xmask&ymask,rs);
        X = X(mask(xmask),:); Y = Y(mask(ymask),:);
    elseif (~isnan(maskcolor))
        mask = reshape(xmask&all(Y~=maskcolor),rs);
        X = X(mask(xmask),:); Y = Y(mask,:);
    else
        mask = reshape(xmask,rs);
        Y = Y(mask,:);
    end
elseif (~isempty(ymask))
    if (~isnan(maskcolor))
        mask = reshape(all(X~=maskcolor)&ymask,rs);
        X = X(mask,:); Y = Y(mask(ymask),:);
    else
        mask = reshape(ymask,rs);
        X = X(mask,:);
    end
elseif (~isnan(maskcolor))
    if (isempty(Y))
        mask = reshape(all(X~=maskcolor),rs);
        X = X(mask,:);
    else
        mask = reshape(all(X~=maskcolor)&all(Y~=maskcolor),rs);
        X = X(mask,:); Y = Y(mask,:);
    end
else
    mask = logical(ones(rs));
end

Xmean = mean(X,1);
for i=1:size(X,2), X(:,i) = X(:,i)-Xmean(i); end

if (~isempty(Y))
    Ymean = mean(Y,1);
    for i=1:size(Y,2), Y(:,i) = Y(:,i)-Ymean(i); end
end

if (nargout>2)
    if (~isempty(Y)), A = zeros(m,size(Y,1),num);
    else, A = zeros(m,m,num);
    end
end

linoffs = mfbox_mkgrid(3*ones(1,length(s)-1))-2;
linoffs = linoffs((linoffs(:,1)>=0)&(sum(linoffs.^2,2)~=0),:);

nl = size(linoffs,1);
if (~isfinite(num)), num = 2*nl; end
maxoffs = step*ceil(num/nl)*max(abs(linoffs));
if (length(s)==2), fullmask = logical(zeros(2*[maxoffs,0]+rs));
else, fullmask = logical(zeros(2*maxoffs+size(mask)));
end
ref = {};
for i=1:length(maxoffs), ref{i} = maxoffs(i)+(1:size(mask,i)); end
fullmask(ref{:}) = mask;

ci = 0;
if (~isempty(Y))
    M = zeros(size(X,2),size(Y,2));
    RX = sqrt(sum(X.^2,1)/N);
    RY = sqrt(sum(Y.^2,1)/N);
    for k=1:num
        lo = step*linoffs(mod(k-1,nl)+1,:)*ceil(k/nl);
        nm = circshift(fullmask,lo)&fullmask;
        om = circshift(fullmask,-lo)&fullmask;
        nm = nm(fullmask); om = om(fullmask);
        P = X(om,:)'*Y(nm,:)/sum(om(:));
        if (normalize)
            dPx = sqrt(sum(X(om,:).^2,1));
            dPy = sqrt(sum(Y(nm,:).^2,1));
            P = sum(om(:))*P./(dPx'*dPy);
        else
            P = P./(RX'*RY);
        end
        if (nargout>2), A(:,:,k) = P; end
        P = P.^2;
        M = M+P; ci = ci+sum(P(:));
    end
    ci = ci/num; M = M/num;
else
    M = zeros(size(X,2));
    RX = sqrt(sum(X.^2,1)/N);
    for k=1:num
        lo = step*linoffs(mod(k-1,nl)+1,:)*ceil(k/nl);
        nm = circshift(fullmask,lo)&fullmask;
        om = circshift(fullmask,-lo)&fullmask;
        nm = nm(fullmask); om = om(fullmask);
        P = X(om,:)'*X(nm,:)/sum(om(:));
        if (normalize)
            dP = sqrt(diag(P));
            P = P./(dP*dP');
        else
            P = P./(RX'*RX);
        end
        if (nargout>2), A(:,:,k) = P+P'; end
        P = P.^2;
        M = M+P; ci = ci+sum(P(logical(1-eye(m))));
    end
    ci = ci/num; M = M/num;
end

return