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)
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