PREPARE LEADFIELD MATRICES
Processes handled in this function as follows
1. Area (optional)
2. Reduce (optional)
3. Spatial smoothing matrix (optional)
4. Load leadfield for all the sessions
5. If 'area_remove' is specified, this area is removed from area
--- Syntax
[Basis,ix, W, ix_ex,L] = vb_prepare_leadfield(basisfile, brainfile, ...
areafile, area_key, reduce, Rfilt, expand_spatial_filter, patch_norm)
[Basis,ix, W, ix_ex,L] = vb_prepare_leadfield(basisfile, brainfile, ...
areafile, area_key, reduce, Rfilt, expand_spatial_filter, ...
patch_norm, area_remove, basis_smoother)
--- Input
See the field names of "bayes_parm"
smoother : basis is smoothed using
= 'std' : standard brain(nextDD, nextIX, xx, xxA)
= 'subj' : Individual brain(subj.nextDD, nextIX, xx, xxA)
--- Output
Basis{Nsession}: Lead field (Nsensor x Ndipole)
W : Smoothing filter or distributed current filter
ix_area : Dipole indices
ix_area_ex : Dipole indices for distributed current. If this is
not empty, W is the distributed current filter.
L : Number of components of basis
--- History
2005-08-16
Okito Yamashita modified
2005-08-22
O. Yamashita ver.30b
2006-03-03 Taku Yoshioka
Support patch size normalization
2006-04-07 Kawawaki, Dai.
2006-08-25 M. Sato
Remove area option is added
2006-11-11 M. Sato
Rmax of vb_spatial_gauss_filter is specified by 'Rfilt'
2008-11-28 Taku Yoshioka
Use vb_disp() for displaying message
2017-03-15 rhayashi
Added spatial_smoother
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function [sbasis, L, W, ix_area, ix_area_ex] = vb_prepare_leadfield( ... 0002 basisfile, brainfile, areafile, area_key, reduce, Rfilt, ... 0003 expand_spatial_filter, patch_norm, area_remove, spatial_smoother) 0004 % PREPARE LEADFIELD MATRICES 0005 % 0006 % Processes handled in this function as follows 0007 % 1. Area (optional) 0008 % 2. Reduce (optional) 0009 % 3. Spatial smoothing matrix (optional) 0010 % 4. Load leadfield for all the sessions 0011 % 5. If 'area_remove' is specified, this area is removed from area 0012 % 0013 % --- Syntax 0014 % [Basis,ix, W, ix_ex,L] = vb_prepare_leadfield(basisfile, brainfile, ... 0015 % areafile, area_key, reduce, Rfilt, expand_spatial_filter, patch_norm) 0016 % [Basis,ix, W, ix_ex,L] = vb_prepare_leadfield(basisfile, brainfile, ... 0017 % areafile, area_key, reduce, Rfilt, expand_spatial_filter, ... 0018 % patch_norm, area_remove, basis_smoother) 0019 % 0020 % --- Input 0021 % See the field names of "bayes_parm" 0022 % smoother : basis is smoothed using 0023 % = 'std' : standard brain(nextDD, nextIX, xx, xxA) 0024 % = 'subj' : Individual brain(subj.nextDD, nextIX, xx, xxA) 0025 % --- Output 0026 % Basis{Nsession}: Lead field (Nsensor x Ndipole) 0027 % W : Smoothing filter or distributed current filter 0028 % ix_area : Dipole indices 0029 % ix_area_ex : Dipole indices for distributed current. If this is 0030 % not empty, W is the distributed current filter. 0031 % L : Number of components of basis 0032 % 0033 % --- History 0034 % 2005-08-16 0035 % Okito Yamashita modified 0036 % 2005-08-22 0037 % O. Yamashita ver.30b 0038 % 2006-03-03 Taku Yoshioka 0039 % Support patch size normalization 0040 % 2006-04-07 Kawawaki, Dai. 0041 % 2006-08-25 M. Sato 0042 % Remove area option is added 0043 % 2006-11-11 M. Sato 0044 % Rmax of vb_spatial_gauss_filter is specified by 'Rfilt' 0045 % 2008-11-28 Taku Yoshioka 0046 % Use vb_disp() for displaying message 0047 % 2017-03-15 rhayashi 0048 % Added spatial_smoother 0049 % Copyright (C) 2011, ATR All Rights Reserved. 0050 % License : New BSD License(see VBMEG_LICENSE.txt) 0051 0052 const = vb_define_verbose; 0053 VERBOSE_LEVEL_NOTICE = const.VERBOSE_LEVEL_NOTICE; 0054 0055 % Leadfield file names for all sessions 0056 if iscell(basisfile), 0057 Nfile = length(basisfile); 0058 else 0059 tmp = basisfile; 0060 Nfile = 1; 0061 basisfile = {tmp}; 0062 end; 0063 0064 vb_disp(['Number of basis files (= Number of sessions): ' ... 0065 num2str(Nfile)],VERBOSE_LEVEL_NOTICE); 0066 0067 % number of dipole vertex 0068 [basis,L]=vb_load_basis(basisfile{1}); 0069 Nj = size(basis,1); 0070 Nv = Nj/L; 0071 0072 %%%% Area for forward model 0073 if ~isempty(areafile) & ~isempty(area_key) 0074 Area = vb_get_area(areafile, area_key); 0075 ix_area = Area.Iextract; 0076 vb_disp(['Area ID: ' area_key],VERBOSE_LEVEL_NOTICE) 0077 vb_disp(['Number of vertices: ' num2str(length(ix_area))], ... 0078 VERBOSE_LEVEL_NOTICE); 0079 else 0080 ix_area = [1:Nv]'; %modified by Kawawaki, Dai at 2006/04/07 0081 vb_disp(['Area ID: All vertices'],VERBOSE_LEVEL_NOTICE) 0082 vb_disp(['Number of vertices: ' num2str(length(ix_area))], ... 0083 VERBOSE_LEVEL_NOTICE); 0084 end 0085 0086 %%%% Remove Area 0087 if exist('area_remove','var') & ~isempty(areafile) & ~isempty(area_remove) 0088 Area = vb_get_area(areafile, area_remove); 0089 ix_remove = Area.Iextract; 0090 % Remove 'ix_remove' from 'ix_area' 0091 ix_area = vb_setdiff2(ix_area, ix_remove); 0092 end 0093 0094 %%%% Reduce vertex 0095 if ~isempty(reduce) & reduce ~=1, 0096 ix_area = vb_reduce_cortex(brainfile, ix_area, reduce); 0097 vb_disp(['Number of reduced vertices: ' ... 0098 num2str(length(ix_area))], VERBOSE_LEVEL_NOTICE); 0099 end 0100 0101 %%%% Spatial Smoothing 0102 if ~isempty(Rfilt) & Rfilt(1) > 0 0103 R = Rfilt(1); 0104 0105 if length(Rfilt) == 2, 0106 Rmax = Rfilt(2); 0107 else 0108 Rmax = 2*R; 0109 end 0110 0111 vb_disp('--- Spatial smoothing filter calculation', ... 0112 VERBOSE_LEVEL_NOTICE); 0113 vb_disp(['R = ' num2str(R,'%2.2e') ', Rmax= ' ... 0114 num2str(Rmax,'%2.2e')],VERBOSE_LEVEL_NOTICE); 0115 0116 [W,ix_area_ex] = vb_spatial_gauss_filter(... 0117 brainfile, R, Rmax, ix_area, expand_spatial_filter, [], [], spatial_smoother); 0118 else 0119 W = speye(length(ix_area), length(ix_area)); 0120 ix_area_ex = ix_area; 0121 end 0122 0123 vb_disp(['Number of vertices = ' num2str(length(ix_area),'%5d')], ... 0124 VERBOSE_LEVEL_NOTICE); 0125 vb_disp(['Number of vertices in expanded area = ' ... 0126 num2str(length(ix_area_ex))],VERBOSE_LEVEL_NOTICE); 0127 0128 %%%% Patch size 0129 if patch_norm==ON, 0130 if strcmpi(spatial_smoother, 'std') 0131 [tmp1,tmp2,tmp3,tmp4,xxA] = vb_load_cortex(brainfile); 0132 else 0133 [tmp1,tmp2,tmp3,tmp4,xxA] = vb_load_cortex(brainfile, 'subj'); 0134 end 0135 clear tmp1 tmp2 tmp3 tmp4 0136 end 0137 0138 %%%% Load lead field matrix 0139 sbasis = cell(Nfile,1); 0140 0141 for n = 1:Nfile 0142 % Load leadfield 0143 vb_disp(['Basis file for session ' num2str(n) ': ' basisfile{n}], ... 0144 VERBOSE_LEVEL_NOTICE); 0145 [basis,L]=vb_load_basis(basisfile{n}); 0146 basis = basis'; 0147 0148 % Patch norm normalization constant 0149 if patch_norm==ON, 0150 S = repmat(xxA(ix_area_ex)',[size(basis,1) 1]); 0151 else 0152 S = ones(size(basis(:,ix_area_ex))); 0153 end 0154 0155 % Spatial Smoothing for leadfield 0156 switch L 0157 case 1, 0158 basis = (basis(:,ix_area_ex).*S)*W; 0159 case 2, 0160 basis = [(basis(:,ix_area_ex).*S)*W (basis(:,Nv+ix_area_ex).*S)*W]; 0161 case 3, 0162 basis = [(basis(:,ix_area_ex).*S)*W (basis(:,Nv+ix_area_ex).*S)*W ... 0163 (basis(:,2*Nv+ix_area_ex).*S)*W]; 0164 end 0165 0166 sbasis{n} = basis; % smoothed basis 0167 end