glm_fillmissingvertices.m

toolbox of BVQX, This is the access between BV and matlab. It will help you to analysis data from BV

function hfile = glm_FillMissingVertices(hfile, srf, method)
% GLM::FillMissingVertices  - fill missing vertices in beta maps
%
% FORMAT:       [glm] = glm.FillMissingVertices(srf [, method])
%
% Input fields:
%
%       srf         surface object (needed for neighborhood information)
%       method      either of 'mean', {'nearest'}, 'dist-weighted'
%              where
%                   mean - interpolate missing by mean of valid neighbors
%                   dist - interpolate by an 1/distance-weighted mean
%                   nearest - set missing value to nearest valid neighbor
%
% Output fields:
%
%       glm         altered GLM structure

% Version:  v0.7b
% Build:    7083022
% Date:     Aug-30 2007, 10:21 PM CEST
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
    numel(hfile) ~= 1 || ...
    numel(srf) ~= 1 || ...
   ~isBVQXfile(hfile, 'glm') || ...
   ~isBVQXfile(srf, 'srf')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
srfc = bvqxfile_getcont(srf.L);
if bc.ProjectType ~= 2 || ...
    bc.ProjectTypeRFX == 0 || ...
    bc.NrOfStudies < 2 || ...
    bc.NrOfVertices ~= srfc.NrOfVertices
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
if nargin < 3 || ...
   ~ischar(method) || ...
    isempty(method) || ...
   ~any(strcmpi(method(:)', ...
       {'all', 'dist', 'dist-weighted', 'mean', 'mean-weighted'}))
    method = 0;
else
    switch (lower(method(:)'))
        %case {'all'}
        %    method = 3;
        case {'dist', 'dist-weighted'}
            method = 2;
        case {'mean', 'mean-weighted'}
            method = 1;
        otherwise
            method = 0;
    end
end

% get numbers, coordinates, and neighbors
numsubs = numel(bc.GLMData.Subject);
nummaps = size(bc.GLMData.Subject(1).BetaMaps, 2);
pcoords = srfc.VertexCoordinate;
neilist = srfc.Neighbors(:, 2);
myeps = eps;

% iterate over subjects
for sc = 1:numsubs
    
    % get BetaMaps
    filltry = 0;
    bmaps = bc.GLMData.Subject(sc).BetaMaps(:, :);
    
    % get indices to fill
    fillidx = find(all(bmaps == 0, 2));
    
    % while indices to fill
    while ~isempty(fillidx)
        
        % get neighborhood information on indices
        fillnei = neilist(fillidx);
        
        % for mean, don't use distance information
        if method == 1
            
            % iterate over indices
            for ic = 1:numel(fillidx)
                
                % get neighbors and their betas
                inei = fillnei{ic};
                bnei = bmaps(inei, :);
                
                % fill into missing
                try
                    knei = find(any(bnei ~= 0, 2));
                    if ~isempty(knei)
                        bmaps(fillidx(ic), :) = mean(bnei(knei, :), 1);
                    end
                catch
                    % do nothing
                end
            end
            
        % either nearest or distance weighted
        else
        
            % iterate over indices
            for ic = 1:numel(fillidx)
                
                % get neighbors and their betas
                inei = fillnei{ic};
                bnei = bmaps(inei, :);
                
                % remove invalid entries
                rnei = find(all(bnei == 0, 2));
                if ~isempty(rnei)
                    inei(rnei) = [];
                    if isempty(inei)
                        continue;
                    end
                    bnei(rnei, :) = [];
                end
                
                % get coordinate positions
                neipos = pcoords(inei, :);
                mispos = repmat(pcoords(fillidx(ic), :), [size(neipos, 1), 1]);
                
                % calculate distance
                pdist = sqrt(sum((neipos - mispos) .^ 2, 2));
                
                % nearest neighbor
                if method == 0
                    
                    % minimum position
                    [srcpos{1:2}] = min(pdist);
                    
                    % replace beta values
                    bmaps(fillidx(ic), :) = bnei(srcpos{2}, :);
                    
                % distance weighted
                else
                    
                    % set betas to weighted mean
                    pdist = 1 ./ (pdist + myeps);
                    bmaps(fillidx(ic), :) = ...
                        sum(double(bnei) .* repmat(pdist, [1, nummaps]), 1) / ...
                        sum(pdist);
                end
            end
        end
        
        % increase fill try
        filltry = filltry + 1;
        if filltry > 3
            warning( ...
                'BVQXfile:BadFileContent', ...
                'Too many beta values missing for subject %d.', ...
                sc ...
            );
            break;
        end
        
        % re-get fillidx
        fillidx = find(all(bmaps == 0, 2));
    end
    
    % set BetaMaps back
    if filltry > 0
        bc.GLMData.Subject(sc).BetaMaps = bmaps;
    end
end

% put back
bvqxfile_setcont(hfile.L, bc);