mfbox_databackend.m

toolbox for spm 5 for data, model free analysis

            if (nargout>2), varargout{3} = tempvol; end
            opdone = 1;
        case 'export_mask'
            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 mask export file ...',f);
                f = fullfile(d,f);
            end
            if (isfield(tempvol,'brainmaskvol')==0)
                [ BG,tempvol ] = mfbox_databackend('getbrainmask',tempvol);
            end
            if (size(f,1)~=0)
                p = strfind(f,'.');
                if (~isempty(p)), f = f(1:p-1); end
                vol = spm_create_vol(struct('fname',[f '.img'], ...
                    'hdrfname',[f '.hdr'],'dim',tempvol.brainmaskdim,'n',1, ...
                    'dt',[spm_type('UINT8'),1],'mat',tempvol.brainmaskvol.mat, ...
                    'descrip','MFBOX mask','private',{{}}));
                dat = mfbox_transform( ...
                    tempvol.mat^(-1)*tempvol.brainmaskmat^(-1)*tempvol.brainmaskvol.mat, ...
                    tempvol.brainmaskdim,dat);
                spm_write_vol(vol,dat);
            end
            if (nargout>0), varargout{1} = f; end
            opdone = 1;
        case 'getdata'
            error(nargchk(2,3,nargin)); error(nargchk(1,1,nargout));
            seq = varargin{2};
            num = 1:seq.timesteps;
            if (nargin>2), num = varargin{3}; end
            v = cat(1,seq.private.private);
            if (length(v(1).dat.dim)==length(seq.dim))
                vn = v(num);
                vn = cat(4,vn.dat);
            else
                vn = v(1).dat(:,:,:,num);
            end
            
            if (isnumeric(vn))
                varargout{1} = reshape(single(vn),[],length(num));
            else
                if (length(num)==1)
                    varargout{1} = reshape(single(numeric(vn)),[],length(num));
                else
                    % this will result in an out of memory error for large data
                    % sets, so we do it by hand
                    
                    varargout{1} = zeros(prod(v(1).dat.dim),length(num));
                    if length(num)>10 % lets display a waitbar
                        wb = true;
                        h = waitbar(0,'Reading fMRI data files...');
                        f = 1/length(num);
                    else
                        wb = false;
                    end
                    for n=1:length(num)
                        vn = v(n);
                        vn = cat(4,vn.dat);
                        varargout{1}(:,n) = reshape(single(numeric(vn)),[],1);
                        if (wb), waitbar(n*f,h); end
                    end
                    if (wb), close(h); end
                end
            end
            opdone = 1;
        case 'init'
            error(nargchk(2,3,nargin)); error(nargchk(1,1,nargout));
            tempvol = varargin{2};
            if (nargin>2), dat = varargin{3}; else dat = struct(); end
            if (isfield(tempvol,'fname'))
                delete(tempvol.fname);
                tempvol = rmfield(tempvol,'fname');
            end
            if (isfield(tempvol,'hdrfname'))
                delete(tempvol.hdrfname);
                tempvol = rmfield(tempvol,'hdrfname');
            end
            tmp = tempname();
            if (isfield(dat,'bg'))
                if (isfield(dat,'bgmat')), mat = dat.bgmat;
                elseif (isfield(dat,'seq')) 
                    if (isfield(dat.seq,'private')), mat = dat.seq.private(1).mat;
                    else mat = eye(4);
                    end
                else mat = eye(4);
                end
                tempvol = spm_create_vol(struct('fname',[tmp '.img'], ...
                    'hdrfname',[tmp '.hdr'],'dim',size(dat.bg),'n',1, ...
                    'dt',[spm_type('UINT16'),1],'mat',mat, ...
                    'descrip','MFBOX tempvol','private',{{}}));
                warning off all
                spm_write_vol(tempvol,dat.bg);
                warning query all
                tempvol = spm_vol(tempvol.fname);
                if ((isfield(dat,'bmvolmat'))&& ...
                        (isfield(dat,'bmmat'))&&(isfield(dat,'bmvoldim')))
                    tempvol.brainmaskmat = dat.bmmat;
                    tempvol.brainmaskvol = struct('mat',dat.bmvolmat,'dim',dat.bmvoldim);
                    tempvol.brainmaskdim = dat.bmvoldim;
                end
            end
            varargout{1} = tempvol;
            opdone = 1;
        case 'getic'
            error(nargchk(4,5,nargin)); error(nargchk(1,2,nargout));
            tempvol = varargin{2};
            data = varargin{3};
            grid = varargin{4};
            if (nargin>4), mask = varargin{5};
            else mask = [];
            end
            if (~isfield(tempvol,'brainmaskvol'))
                [BG,tempvol] = mfbox_databackend('getbrainmask',tempvol);
            end
            transform = tempvol.brainmaskvol.mat^(-1)*tempvol.brainmaskmat*tempvol.mat;
            c = all(grid>=1,1);
            for i=1:size(grid,1), c = c&(grid(i,:)<=tempvol.dim(i)); end
            IC = reshape(accumarray(round(grid(:,c))',double(data(c)), ...
                tempvol.dim,@mean,NaN),tempvol.dim); % matlab 7.0 does not like single in accumarray
            if (any(any(transform-eye(4)~=0)))
                tr = transform^(-1);
                gr = tr(1:(end-1),1:(end-1))*mfbox_mkgrid(tempvol.brainmaskvol.dim)';
                for i=1:(size(tr,1)-1), gr(i,:) = gr(i,:)+tr(i,end); end
                gr = mat2cell(gr,ones(3,1),prod(tempvol.brainmaskvol.dim));
                IC = single(reshape(interpn(IC,gr{:},'nearest'),tempvol.brainmaskvol.dim));
                if (~isempty(mask))
                    mask = (mfbox_transform(tr, ...
                        tempvol.brainmaskvol.dim,mask*1)<0.75);
                    IC(mask) = single(NaN);
                end
            end
            varargout{1} = IC;
            if (nargout>1), varargout{2} = tempvol; end
            opdone = 1;
        case 'getbg'
            error(nargchk(2,2,nargin)); error(nargchk(1,2,nargout));
            tempvol = varargin{2};
            if (~isfield(tempvol,'brainmaskvol'))
                [BG,tempvol] = mfbox_databackend('getbrainmask',tempvol);
            end
            varargout{1} = mfbox_transform( ... % besser spm_sample_vol
                tempvol.mat^(-1)*tempvol.brainmaskmat^(-1)*tempvol.brainmaskvol.mat, ...
                tempvol.brainmaskvol.dim,numeric(tempvol.private.dat));
            if (nargout>1), varargout{2} = tempvol; end
            opdone = 1;
        case 'getbrainmask'
            error(nargchk(1,2,nargin)); error(nargchk(1,2,nargout));
            bmVol = spm_vol(fullfile(spm('Dir'),'apriori','brainmask.nii'));
            if (nargin>1)
                tempvol = varargin{2};
                if (~isfield(tempvol,'brainmask'))
                    a = questdlg(sprintf( ...
                        'Data registered on\nstandard brainmask?\n'), ...
                        'Get brainmask','Yes','No','Yes');
                    if (strcmp(a,'No'))
                        V = spm_coreg(tempvol,bmVol,struct('graphics',0));
                        tempvol.brainmaskmat = spm_matrix(V(:)');
                    else
                        
                        tempvol.brainmaskmat = eye(4);
                    end
                    bm = numeric(bmVol.private.dat);
                    tempvol.brainmaskdim = size(bm);
                    tempvol.brainmaskvol = bmVol;
                    tempvol.brainmask = mfbox_transform( ...
                        tempvol.brainmaskvol.mat^(-1)*tempvol.brainmaskmat*tempvol.mat, ...
                        tempvol.dim,bm);
                end
                varargout{1} = tempvol.brainmask;
                if (nargout>1), varargout{2} = tempvol; end
            else
                varargout{1} = numeric(bmVol.private.dat);
            end
            opdone = 1;
        case 'startmaskedit'
            error(nargchk(3,3,nargin)); error(nargchk(1,1,nargout));
            mask = varargin{2};
            tempvol = varargin{3};
            if (min(mask(:))==max(mask(:))), mask(1:2)=[0,1]; end        % hack
            tmp = tempname();
            maskvol = spm_create_vol(struct('fname',[tmp '.img'], ...
                'hdrfname',[tmp '.hdr'],'dim',tempvol.dim,'n',1, ...
                'dt',[spm_type('UINT16'),1],'mat',tempvol.mat, ...
                'descrip','MFBOX maskedit','private',{{}}));
            maskvol = spm_write_vol(maskvol,mask);
            tempvol.maskvol = maskvol;
            spm_orthviews('rmblobs',tempvol.volhandle);
            spm_orthviews('roi','init',tempvol.volhandle,maskvol);
            varargout{1} = tempvol;
            opdone = 1;
        case 'stopmaskedit'
            error(nargchk(2,2,nargin)); error(nargchk(1,2,nargout));
            tempvol = varargin{2};
            spm_orthviews('roi','save',tempvol.volhandle);
            varargout{1} = numeric(tempvol.maskvol.private.dat);
            delete(tempvol.maskvol.fname);
            delete(tempvol.maskvol.hdrfname);
            tempvol = rmfield(tempvol,'maskvol');
            spm_orthviews('rmblobs',tempvol.volhandle);
            if (nargout>1), varargout{2} = tempvol; end
            opdone = 1;
        case 'redraw_mask'
            error(nargchk(3,3,nargin)); error(nargchk(0,0,nargout));
            mask = varargin{2};
            tempvol = varargin{3};
            spm_orthviews('rmblobs',tempvol.volhandle);
            a = size(mask);
            c = ones(sum(mask(:)),1);
            if (size(c,1)~=0)
                g = mfbox_mkgrid(a(1:3))';
                spm_orthviews('addcolouredblobs',tempvol.volhandle, ...
                    g(:,mask(:)),c,tempvol.mat,[1,1,1]');
                spm_orthviews('redraw',tempvol.volhandle);
            end
            opdone = 1;
        case 'getposition'
            error(nargchk(2,2,nargin)); error(nargchk(0,2,nargout));
            tempvol = varargin{2};
            p = spm_orthviews('pos');
            p = tempvol.mat^(-1)*[p(:);1];
            p = p(1:3);
            if (nargout>0)
                varargout{1} = tempvol;
                if (nargout>1), varargout{2} = p; end
            end
            opdone = 1;
        case 'redraw'
            error(nargchk(2,5,nargin)); error(nargchk(0,1,nargout));
            tempvol = varargin{2};
            vols = varargin{3};
            if (nargin>3), blobs = varargin{4}; else blobs = {}; end
            if (nargin>4), pos = varargin{5}; else pos = []; end
            if (~isempty(pos))
                pos = tempvol.mat*[pos(:);1];
                pos = pos(1:3);
            end
            if (~isempty(vols))
                if (length(vols)>1), spm_image('init',vols{1}); end
                if (isfield(tempvol,'volhandle'))
                    spm_orthviews('delete',tempvol.volhandle);
                end
                tempvol.volhandle = spm_orthviews('image',vols{1},[0.0 0.45 1 0.55]);
            end
            spm_orthviews('rmblobs',tempvol.volhandle);
            for i=1:length(blobs)
                b = blobs{i};
                spm_orthviews('addcolouredblobs',tempvol.volhandle, ...
                    b{1},b{2},tempvol.mat,b{3});
            end
            if (~isempty(pos)), spm_orthviews('setcoords',pos);
            else spm_orthviews('redraw',tempvol.volhandle);
            end
            if (nargout>0), varargout{1} = tempvol; end
            opdone = 1;
    end
end

if (opdone==1)
    return;
end

switch op
    case 'load'
        error(nargchk(2,4,nargin)); error(nargchk(1,3,nargout));
        mfbss = varargin{2};
        if (nargin>2), f = varargin{3};
        else f = '';
        end
        if ((isempty(f))||(nargin<=3))
            [f,d] = uigetfile('*.mfb','Select MFBOX file ...',f);
            f = fullfile(d,f);
        end
        if ((size(f,1)~=0)&&(exist(f,'file')>0))
            v = load(f,'-mat');
            isvalid = true;
            if (isfield(v.mfbss,'version'))
                if (~strcmp(v.mfbss.version,mfbox_version()))
                    isvalid = false;
                end
            else