vb_step_filter

 Calculate step function filter on cortical surface 
 --- Syntax
 function [CW ,ixall] = spatial_average_filter(nextDD,nextIX,R,Index)

 --- Input
 nextDD      : Distance to neighbour vertices
 nextIX      : Indices of neighbour vertices
 R           : Smoothing filter radius
 Index       : Vertex index for selected area

 --- Optional Input
 flag = 0    : output region = input region
      = 1    : extend output region to neighbor of edge (Default)

 --- Output
 CW          : Smoothing Gaussian filter
 ixall       : Vertex index for smoothing filter
               smoothed_leadfield(:,Index) = leadfield(:,ixall) * CW

 --- History
 2007-02-23 Taku Yoshioka

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function [CW ,ixall] = vb_step_filter(nextDD,nextIX,R,Index)
0002 % Calculate step function filter on cortical surface
0003 % --- Syntax
0004 % function [CW ,ixall] = spatial_average_filter(nextDD,nextIX,R,Index)
0005 %
0006 % --- Input
0007 % nextDD      : Distance to neighbour vertices
0008 % nextIX      : Indices of neighbour vertices
0009 % R           : Smoothing filter radius
0010 % Index       : Vertex index for selected area
0011 %
0012 % --- Optional Input
0013 % flag = 0    : output region = input region
0014 %      = 1    : extend output region to neighbor of edge (Default)
0015 %
0016 % --- Output
0017 % CW          : Smoothing Gaussian filter
0018 % ixall       : Vertex index for smoothing filter
0019 %               smoothed_leadfield(:,Index) = leadfield(:,ixall) * CW
0020 %
0021 % --- History
0022 % 2007-02-23 Taku Yoshioka
0023 %
0024 % Copyright (C) 2011, ATR All Rights Reserved.
0025 % License : New BSD License(see VBMEG_LICENSE.txt)
0026 
0027 Nvertex = length(Index);    % # of points in selected area
0028 Nall    = length(nextIX);    % # of points in whole brain
0029 rIndex = Index;   
0030 
0031 % Calculate total number of neighbors (Nnext)
0032 Nnext   = 0;
0033 
0034 for n=1:Nvertex,
0035   inx   = find( nextDD{rIndex(n)} <= R );    % neighbor within Rmax
0036   Nnext = Nnext + length(inx);
0037 end;
0038 
0039 % Initialize variable
0040 val     = zeros(Nnext,1);    % Filter value
0041 ipoint    = zeros(Nnext,1);    % Center point index
0042 inext    = zeros(Nnext,1);    % Neighbor index
0043 Nbegin  = 1;
0044 
0045 % Find neighbor within Rmax and calculate Gaussian filter value
0046 for n=1:Nvertex,
0047     i    = rIndex(n);        % center vertex index
0048     dd0 = nextDD{i};        % Neighbor distance
0049     inx = find( dd0 <= R );            % Find neighbor within Rmax
0050     
0051     ix = nextIX{i}(inx);        % Neighbor index within Rmax
0052     dd0 = dd0(inx);            % Neighbor distance
0053     
0054     % Filter value
0055     vdd = ones(length(inx),1);      % Spatial filter
0056     vdd = vdd/sum(vdd);        % Normalization
0057     
0058     % Total # of neighbor points
0059     Nend   = Nbegin + length(inx) - 1;
0060     Indx   = Nbegin:Nend;
0061     Nbegin = Nend + 1;
0062     
0063     ipoint(Indx) = n;        % center vertex index
0064     inext(Indx)  = ix ;     % Neighbor index
0065     val(Indx)    = vdd;        % Gaussian filter value
0066 
0067 end;
0068 
0069 ipoint = ipoint(1:Nend);
0070 inext  = inext(1:Nend) ;
0071 val    = val(1:Nend)   ;
0072     
0073 CW    = sparse( inext , ipoint , val , Nall, Nvertex) ;
0074 ixall = unique(inext);
0075 CW    = CW(ixall,:);