vb_prepare_leadfield

 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)

 --- Input
 See the field names of "bayes_parm"

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

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