Home > functions > job > job_edit_act_dir > job_edit_act_spatial_gauss_filter.m

job_edit_act_spatial_gauss_filter

PURPOSE ^

Calculate gasussian smoothing filter of 2D cortical surface

SYNOPSIS ^

function [CW ,ixall] = job_edit_act_spatial_gauss_filter(nextDD,nextIX,R,Rmax,Index,flag,cosval,xx)

DESCRIPTION ^

 Calculate gasussian smoothing filter of 2D cortical surface 
 --- Syntax
 function [CW ,ixall] =
 vb_spatial_gauss_filter(brainfile,Rradius,Rmax,Index,flag,cosval)

 --- Input
 brainfile   : brain file 
 Rradius     : Gaussian smoothing filter radius
 Rmax        : Max radius of Gaussian smoothing filter
 Index       : Vertex index for selected area
 --- Optional Input
 flag = 0    : output region = input region
      = 1    : extend output region to neighbor of edge (Default)
 cosval : Threshold value for normal direction cos = cos(limit_angle) [-1]
 --- Output
 CW          : Smoothing Gaussian filter
 ixall       : Vertex index for smoothing filter
               smoothed_leadfield(:,Index) = leadfield(:,ixall) * CW

 Ver 1.0 written by M. Sato  2003-4-15
 Modified by Taku Yoshioka
 2005-3-29 O.Yamashita
   Bug fix : the program does not run due to 'Index' 
             after reudce_neighbor_data, when flag = 0. 
   Error check for parm.cosval.
   The comment is modified.  
 2005-3-31 O.Yamashita
   Bug fix : 'ixall' is not output correctly, when flag = 0. 
 2005-08-09 O.Yamashita
   input arguments are modified

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [CW ,ixall] = job_edit_act_spatial_gauss_filter(nextDD,nextIX,R,Rmax,Index,flag,cosval,xx)
0002 % Calculate gasussian smoothing filter of 2D cortical surface
0003 % --- Syntax
0004 % function [CW ,ixall] =
0005 % vb_spatial_gauss_filter(brainfile,Rradius,Rmax,Index,flag,cosval)
0006 %
0007 % --- Input
0008 % brainfile   : brain file
0009 % Rradius     : Gaussian smoothing filter radius
0010 % Rmax        : Max radius of Gaussian smoothing filter
0011 % Index       : Vertex index for selected area
0012 % --- Optional Input
0013 % flag = 0    : output region = input region
0014 %      = 1    : extend output region to neighbor of edge (Default)
0015 % cosval : Threshold value for normal direction cos = cos(limit_angle) [-1]
0016 % --- Output
0017 % CW          : Smoothing Gaussian filter
0018 % ixall       : Vertex index for smoothing filter
0019 %               smoothed_leadfield(:,Index) = leadfield(:,ixall) * CW
0020 %
0021 % Ver 1.0 written by M. Sato  2003-4-15
0022 % Modified by Taku Yoshioka
0023 % 2005-3-29 O.Yamashita
0024 %   Bug fix : the program does not run due to 'Index'
0025 %             after reudce_neighbor_data, when flag = 0.
0026 %   Error check for parm.cosval.
0027 %   The comment is modified.
0028 % 2005-3-31 O.Yamashita
0029 %   Bug fix : 'ixall' is not output correctly, when flag = 0.
0030 % 2005-08-09 O.Yamashita
0031 %   input arguments are modified
0032 %
0033 % Copyright (C) 2011, ATR All Rights Reserved.
0034 % License : New BSD License(see VBMEG_LICENSE.txt)
0035 
0036 Nvertex = length(Index);    % # of points in selected area
0037 Nall    = length(nextIX);    % # of points in whole brain
0038 
0039 if nargin < 7, cosval = -1; end
0040 if nargin < 6, flag = 1; end;
0041 
0042 if flag == 1 
0043     rIndex = Index;   
0044 else
0045     % Restrict neighbor vertex into input index region
0046     [nextIX,nextDD] = ...
0047         vb_reduce_neighbor_data(Index,nextIX,nextDD);
0048     rIndex = [1:Nvertex]'; %rename index
0049 end;
0050 
0051 % Threshold value for normal direction cos
0052 if cosval > 1 % avoid empty index in for loop
0053     cosval = -1;
0054 end
0055 
0056 %%%%%%% Smoothing Filter calculation %%%%%%%
0057 if    R <= 0,
0058     % No smoothing
0059     CW    = speye(Nvertex,Nvertex);
0060     ixall = Index;  
0061     return
0062 end
0063 
0064 % Calculate total number of neighbors (Nnext)
0065 Nnext   = 0;
0066 
0067 for n=1:Nvertex,
0068   inx   = find( nextDD{rIndex(n)} <= Rmax );    % neighbor within Rmax
0069   Nnext = Nnext + length(inx);
0070 end;
0071 
0072 % Initialize variable
0073 val     = zeros(Nnext,1);    % Filter value
0074 ipoint    = zeros(Nnext,1);    % Center point index
0075 inext    = zeros(Nnext,1);    % Neighbor index
0076 Nbegin  = 1;
0077 
0078 % Find neighbor within Rmax and calculate Gaussian filter value
0079 for n=1:Nvertex,
0080     i    = rIndex(n);                % center vertex index
0081     dd0 = nextDD{i};            % Neighbor distance
0082     inx = find( dd0 <= Rmax );    % Find neighbor within Rmax
0083     
0084     ix    = nextIX{i}(inx);        % Neighbor index within Rmax
0085     dd0 = dd0(inx);                % Neighbor distance
0086     
0087     % Check direction cos is larger than cosval
0088     if cosval>-1, 
0089       dcos= xx(ix,:)*xx(i,:)';
0090       inx = find( dcos > cosval );
0091       ix  = ix(inx);
0092     else
0093       inx = 1:length(ix);
0094     end
0095     
0096     % Gaussian filter value
0097     vdd = exp( - (dd0(inx)/R).^2 );    % Gaussian filter
0098     vdd = vdd/sum(vdd);            % Normalization
0099     
0100     % Total # of neighbor points
0101     Nend   = Nbegin + length(inx) - 1;
0102     Indx   = Nbegin:Nend;
0103     Nbegin = Nend + 1;
0104     
0105     ipoint(Indx) = n;        % center vertex index
0106     inext(Indx)  = ix ;     % Neighbor index
0107     val(Indx)    = vdd;        % Gaussian filter value
0108 
0109 end;
0110 
0111 ipoint = ipoint(1:Nend);
0112 inext  = inext(1:Nend) ;
0113 val    = val(1:Nend)   ;
0114     
0115 CW    = sparse( inext , ipoint , val , Nall, Nvertex) ;
0116 ixall = unique(inext);
0117 CW    = CW(ixall,:);        
0118 if flag == 0, ixall = Index; end % index before renamed.
0119

Generated on Tue 27-Aug-2013 11:46:04 by m2html © 2005