mfbox_databackend.m

toolbox for spm 5 for data, model free analysis

function varargout=mfbox_databackend(varargin)
% in/output functions
%
% Usage:
%  [output,...]=mfbox_databackend(command,input,...)
%
%  available commands:
%          export           - input:X,private,filename,force output:filename
%                             export data in spm format
%          export_mask      - input:private,data,filename,force output:filename
%                             export mask in spm format
%          export_result    - input:private,data,filename,force output:filename
%                             export components in spm format
%          getatlas         - input:privat output:atlas,private
%                             get 3 atlas components in std space
%          getbg            - input:private output:background,private
%                             get background in std space
%          getbrainmask     - input:private,force output:brainmask,private
%          getdata          - input:sequence,num output:data
%                             get timesteps in sequence data
%          getic            - input:private,data,grid,mask output:ic,private
%                             get ic in std space
%          getposition      - input:private output:private,position
%          import           - input:filename,force output:sequence,filename
%                             import data in spm format
%          import_bg        - input:private,filename,force output:background,private
%          import_design    - input:filename,force output:design,filename
%          import_mask      - input:filename,private output:mask,filename,pivate
%          import_reference - input:private,filename,force output:reference,filename,private
%          import_timemask  - input:filename,force output:timemask,filename
%          init             - input:private,init_data output:private
%                             initialize private data possibly using data from load
%          load             - input:mfbss,filename,force output:mfbss,filename,init_data
%                             load mfbss structure
%          loadcompare      - input:mfbss_compare,filename,force output:mfbss_compare,filename
%                             load mfbss_compare structure
%          load_params      - input:filename,force output:params
%                             load parameter file
%          redraw           - input:private,vols,blobs,pos output:private
%          redraw_mask      - input:mask,privat output:
%          save             - input:mfbss,private,filename,force output:filename
%                             save mfbss structure
%          savecompare      - input:mfbss_compare,filename,force output:filename
%                             save mfbss_compare structure
%          save_params      - input:params,filename,force output:filename
%                             save parameter file
%          startmaskedit    - input:mask,private output:private
%          stopmaskedit     - input:private output:mask,private
%
%  private  - private data from init
%  force    - interactive switch
%  X        - (NxT) data
%  mfbss    - mfbss structure (see mfbox_init)
%  mask     - data mask (XxYxZ) with sum(mask(:))==N
%  grid     - 3d positions (3xN) of the data values
%  sequence - sequence from import
%
% Copyright by Peter Gruber, Fabian J. Theis and Ingo R. Keck
% 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

if (nargin<1), return; end

opdone = 0;

op = varargin{1};

if (~isempty(which('spm')))

    switch op
        case 'import' % import data sequence, optional argument is the filename
            error(nargchk(1,3,nargin)); error(nargchk(1,2,nargout));
            if (nargin>1), f = varargin{2};
            else f = [];
            end
            if ((isempty(f))||(nargin<=2))
                f = spm_select(Inf,'image','Select fMRI sequence ...',f);
            end
            if (size(f,1)==1)
                if (strcmp(f(1,end-5:end),'.nii,1')), f = f(1,1:end-2); end
            end
            if (size(f,1)~=0)
                varargout{1} = mfbox_read_spmfiles(f);
                if (nargout>1), varargout{2} = f; end
            else
                varargout{1} = struct('num',0,'dim',[0,0],'timesteps',0);
                if (nargout>1), varargout{2} = f; end
            end
            opdone = 1;
        case 'export' % exports (preprocessed) data sequence, first argument is the data, optional the filename
            error(nargchk(3,5,nargin)); error(nargchk(0,1,nargout));
            dat = varargin{2};
            tempvol = varargin{3};
            if (nargin>3), f = varargin{4};
            else f = '';
            end
            if ((isempty(f))||(nargin<=4))
                [fn,d] = uiputfile('*.nii','Select file to save preprocessed data ...');
                f = fullfile(d,fn);
            end
            if (~isempty(fn))
                if (iscell(dat)~=1), X = dat; mask = [];
                elseif (length(dat)>1), X = dat{1}; mask = dat{2}; grid = dat{3};
                else X = dat{1}; mask = [];
                end
                if (isstruct(X))
                    oX = X;
                    X = zeros(sum(mask(:)),n,'single');
                    for i=1:n
                        v = reshape(mfbox_databackend('getdata',oX,i),[],1);
                        X(:,i) = v(mask);
                    end
                end
                if (~isempty(mask))
                    s = size(X);
                    s = [size(mask),s(end)];
                    X = reshape(X,[],s(end));
                    Y = zeros(prod(s(1:(end-1))),s(end),'single');
                    g = mat2cell(round(grid),ones(1,size(grid,1)),size(grid,2));
                    nmask = sub2ind(size(mask),g{:});
                    if (any(find(mask(:))~=nmask(:)))
                        for i=1:s(end)
                            Y(:,i) = reshape(single(accumarray(round(grid)', ...
                                double(X(:,i)),s(1:(end-1)),@mean)),[],1); % matlab 7.0 does not like single in accumarray
                            Y(~mask,i) = NaN;
                        end
                    else
                        for i=1:s(end)
                            Y(mask,:) = X;
                            Y(~mask,i) = NaN;
                        end
                    end
                    mfbox_write_spmfiles(f,Y,tempvol);
                else
                    mfbox_write_spmfiles(f,X,tempvol);
                end
            end
            if (nargout>0), varargout{1} = f; end
            opdone = 1;
        case 'import_reference'
            error(nargchk(2,4,nargin)); error(nargchk(1,3,nargout));
            tempvol = varargin{2};
            if (nargin>2), f = varargin{3};
            else f = [];
            end
            if ((isempty(f))||(nargin<=3))
                f = spm_select(1,'image','Select reference ...',f);
            end
            ref = [];
            if (size(f,1)~=0)
                ref = spm_vol(f(1,:));
                a = questdlg(sprintf( ...
                    'Reference registered on\nstandard brainmask or data?\n'), ...
                    'Import reference','Brainmask','Data','Brainmask');
                dat = numeric(ref.private.dat);
                if (length(size(dat))==4)
                    dat = squeeze(dat(:,:,:,1));
                end
                if (strcmp(a,'Brainmask'))
                    if (~isfield(tempvol,'brainmaskvol'))
                        [BG,tempvol] = mfbox_databackend('getbrainmask',tempvol);
                    end
                    ref = (mfbox_transform( ...
                        ref.mat^(-1)*tempvol.brainmaskmat*tempvol.mat, ...
                        tempvol.dim,dat)>0.5);
                else
                    ref = (mfbox_transform(ref.mat^(-1)*tempvol.mat, ...
                        tempvol.dim,dat)>0.5);
                end
            end
            varargout{1} = ref;
            if (nargout>1), varargout{2} = f; end
            if (nargout>2), varargout{3} = tempvol; end
            opdone = 1;
        case 'import_bg'
            error(nargchk(2,4,nargin)); error(nargchk(1,2,nargout));
            tempvol = varargin{2};
            if (nargin>2), f = varargin{3};
            else f = [];
            end
            if ((isempty(f))||(nargin<=3))
                f = spm_select(1,'image','Select background ...',f);
            end
            ref = [];
            if (size(f,1)~=0)
                ref = spm_vol(f(1,:));
                a = questdlg(sprintf( ...
                    'Background registered on\nstandard brainmask or data?\n'), ...
                    'Import background','Brainmask','Data','Brainmask');
                if (strcmp(a,'Brainmask'))
                    if (~isfield(tempvol,'brainmaskvol'))
                        [BG,tempvol] = mfbox_databackend('getbrainmask',tempvol);
                    end
                    ref = mfbox_transform( ...
                        ref.mat^(-1)*tempvol.brainmaskmat*tempvol.mat, ...
                        tempvol.dim,numeric(ref.private.dat));
                else
                    ref = mfbox_transform(ref.mat^(-1)*tempvol.mat, ...
                        tempvol.dim,numeric(ref.private.dat));
                end
            end
            ref = reshape(ref,tempvol.dim);
            varargout{1} = ref;
            if (nargout>1), varargout{2} = tempvol; end
            opdone = 1;
        case 'export_result'
            error(nargchk(3,5,nargin)); error(nargchk(0,1,nargout));
            tempvol = varargin{2};
            dat = varargin{3};
            if (nargin>3), f = varargin{4};
            else f = '';
            end
            if ((isempty(f))||(nargin<=4))
                %[f,d] = uiputfile('*.img','Select result export file ...',f);
                [f,d] = uiputfile('*.nii','Select result export file ...',f);
                f = fullfile(d,f);
            end
            if (size(f,1)~=0)
                p = strfind(f,'.');
                if (~isempty(p)), f = f(1:p-1); end
                dat = int16(((2^15)-1)*dat/max(max(dat(:)),-min(dat(:))));
                % in former times we exported .img and .hdr:
                %vol = spm_create_vol(struct('fname',[f '.img'], ...
                %    'hdrfname',[f '.hdr'],'dim',tempvol.dim,'n',1, ...
                %    'dt',[spm_type('INT16'),1],'mat',tempvol.mat, ...
                %    'descrip','MFBOX result','private',{{}}));
                % now is time for .nii. 
                % spm_create_vol simply looks for the fname extention to
                % decide...
                vol = spm_create_vol(struct('fname',[f '.nii'], ...
                    'dim',tempvol.dim,'n',1, ...
                    'dt',[spm_type('INT16'),1],'mat',tempvol.mat, ...
                    'descrip','MFBOX result','private',{{}}));
                spm_write_vol(vol,dat);
            end
            if (nargout>0), varargout{1} = f; end
            opdone = 1;
        case 'import_mask'
            error(nargchk(2,3,nargin)); error(nargchk(1,3,nargout));
            f = varargin{2};
            if (nargin>2), tempvol = varargin{3};
            else tempvol = [];
            end
            if (isempty(f))
                f = spm_select(1,'image','Select mask ...',f);
                a = questdlg(sprintf( ...
                    'Mask registered on\nstandard brainmask or data?\n'), ...
                    'Import reference','Brainmask','Data','Brainmask');
            elseif (iscell(f))
                a = f{2};
                f = f{1};
            else
                a = questdlg(sprintf( ...
                    'Mask registered on\nstandard brainmask or data?\n'), ...
                    'Import reference','Brainmask','Data','Brainmask');
            end
            if (size(f,1)~=0)
                f = f(1,:);
                bmVol = spm_vol(f);
                if (~isempty(tempvol))
                    if (strcmp(a,'Brainmask'))
                        if (~isfield(tempvol,'brainmaskvol'))
                            [BG,tempvol] = mfbox_databackend('getbrainmask',tempvol);
                        end
                        mask = mfbox_transform( ...
                            bmVol.mat^(-1)*tempvol.brainmaskmat*tempvol.mat, ...
                            tempvol.dim,numeric(bmVol.private.dat));
                    else
                        mask = mfbox_transform( ...
                            bmVol.mat^(-1)*tempvol.mat, ...
                            tempvol.dim,numeric(bmVol.private.dat));
                    end
                else
                    mask = numeric(bmVol.private.dat);
                end
            else
                mask = true(tempvol.dim);
            end
            varargout{1} = mask;
            if (nargout>1), varargout{2} = {f,a}; end