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