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