Calculate intersubject average of estimated cortical current.
(VBMEG public function)
[syntax]
vb_job_current_grandaverage(proj_root,avr_parm)
[input]
proj_root: <<string>> VBMEG project root directory.
avr_parm : <<struct>> Parameters for intersubject averaging.
--- fields of avr_parm
curr_file : <<cell>> Cortical current files.
brain_file : <<cell>> Cortical surface model files.
key : <<cell>> ID of coregistration matrix for each subjects
corresponding to brain files. In common, coregistration
matrix from the subject to a template cortical model is
specified.
reduce : <optional> <<double>> Reduction ratio of vertices of
the cortical surface model on which grandaverage is
taken.
tbrain_file : <optional> <<string>> Filename of cortical surface
model on which grandaverage is taken. This file is used
to reduce the number of vertices.
curr_file_out: <<string>> Filename of cortical current file.
---
[output]
Cortical current file (.curr.mat) on the template cortical surface
model (curr_file_out).
[example]
>> avr_parm.curr_file = {'./subjects/0001/result/0001_LR.curr.mat', ...
'./subjects/0002/result/0002_LR.curr.mat', ...
'./subjects/0003/result/0003_LR.curr.mat', ...
'./subjects/0004/result/0004_LR.curr.mat'};
>> avr_parm.brain_file ...
= {'./subjects/0001/brain/0001.brain.mat', ...
'./subjects/0002/brain/0002.brain.mat', ...
'./subjects/0003/brain/0003.brain.mat', ...
'./subjects/0004/brain/0004.brain.mat'};
>> avr_parm.key = {'fsaverage.sphere.reg','fsaverage.sphere.reg', ...
'fsaverage.sphere.reg','fsaverage.sphere.reg'};
>> avr_parm.tbrain_file = './subjects/fsaverage/fsaverage.brain.mat';
>> avr_parm.curr_file_out ...
= './subjects/fsaverage/result/fsaverage_LR.curr.mat';
>> vb_job_grandaverage(proj_root,avr_parm);
[history]
2010-12-01 Taku Yoshioka
2011-01-27 taku-y
[minor] 'reduce' parameter supported. If reduce is not given or
reduce==1.0, this function give the same result with that of the
previous version.
2011-11-15 taku-y
[minor] Bug fixed for insufficient number of current dipoles, due to
small number of Z-current dipoles not covering whole cortex.
2012-02-08 taku-y
[minor] Bug fixed for insufficient number of current dipoles, due to
small number of Z-current dipoles not covering whole cortex (first
loop).
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function vb_job_grandaverage(proj_root,avr_parm) 0002 % Calculate intersubject average of estimated cortical current. 0003 % (VBMEG public function) 0004 % 0005 % [syntax] 0006 % vb_job_current_grandaverage(proj_root,avr_parm) 0007 % 0008 % [input] 0009 % proj_root: <<string>> VBMEG project root directory. 0010 % avr_parm : <<struct>> Parameters for intersubject averaging. 0011 % --- fields of avr_parm 0012 % curr_file : <<cell>> Cortical current files. 0013 % brain_file : <<cell>> Cortical surface model files. 0014 % key : <<cell>> ID of coregistration matrix for each subjects 0015 % corresponding to brain files. In common, coregistration 0016 % matrix from the subject to a template cortical model is 0017 % specified. 0018 % reduce : <optional> <<double>> Reduction ratio of vertices of 0019 % the cortical surface model on which grandaverage is 0020 % taken. 0021 % tbrain_file : <optional> <<string>> Filename of cortical surface 0022 % model on which grandaverage is taken. This file is used 0023 % to reduce the number of vertices. 0024 % curr_file_out: <<string>> Filename of cortical current file. 0025 % --- 0026 % 0027 % [output] 0028 % Cortical current file (.curr.mat) on the template cortical surface 0029 % model (curr_file_out). 0030 % 0031 % [example] 0032 % >> avr_parm.curr_file = {'./subjects/0001/result/0001_LR.curr.mat', ... 0033 % './subjects/0002/result/0002_LR.curr.mat', ... 0034 % './subjects/0003/result/0003_LR.curr.mat', ... 0035 % './subjects/0004/result/0004_LR.curr.mat'}; 0036 % >> avr_parm.brain_file ... 0037 % = {'./subjects/0001/brain/0001.brain.mat', ... 0038 % './subjects/0002/brain/0002.brain.mat', ... 0039 % './subjects/0003/brain/0003.brain.mat', ... 0040 % './subjects/0004/brain/0004.brain.mat'}; 0041 % >> avr_parm.key = {'fsaverage.sphere.reg','fsaverage.sphere.reg', ... 0042 % 'fsaverage.sphere.reg','fsaverage.sphere.reg'}; 0043 % >> avr_parm.tbrain_file = './subjects/fsaverage/fsaverage.brain.mat'; 0044 % >> avr_parm.curr_file_out ... 0045 % = './subjects/fsaverage/result/fsaverage_LR.curr.mat'; 0046 % >> vb_job_grandaverage(proj_root,avr_parm); 0047 % 0048 % [history] 0049 % 2010-12-01 Taku Yoshioka 0050 % 2011-01-27 taku-y 0051 % [minor] 'reduce' parameter supported. If reduce is not given or 0052 % reduce==1.0, this function give the same result with that of the 0053 % previous version. 0054 % 2011-11-15 taku-y 0055 % [minor] Bug fixed for insufficient number of current dipoles, due to 0056 % small number of Z-current dipoles not covering whole cortex. 0057 % 2012-02-08 taku-y 0058 % [minor] Bug fixed for insufficient number of current dipoles, due to 0059 % small number of Z-current dipoles not covering whole cortex (first 0060 % loop). 0061 % 0062 % Copyright (C) 2011, ATR All Rights Reserved. 0063 % License : New BSD License(see VBMEG_LICENSE.txt) 0064 0065 0066 % 0067 % Input arguments 0068 % 0069 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0070 vb_struct2vars(avr_parm,{'curr_file','brain_file','key','reduce', ... 0071 'tbrain_file','curr_file_out'}); 0072 0073 if isempty(reduce), 0074 reduce = 1.0; 0075 end 0076 0077 % 0078 % Load first file 0079 % 0080 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0081 [Jinfo,J] = vb_load_current(curr_file{1},0,true); 0082 C = vb_load_cortex_pmat(brain_file{1},key{1}); 0083 0084 % added by taku-y 2012-02-08 0085 Jtmp = zeros(size(C,2),size(J,2)); 0086 Jtmp(Jinfo.ix_act_ex,:) = J; 0087 J = Jtmp; 0088 0089 if reduce==1.0, 0090 Javr = C*J; 0091 else 0092 brain_file_tmp = tbrain_file; 0093 0094 % Reduce vertex of the cortical surface model 0095 V = vb_load_cortex(brain_file_tmp); 0096 ix_act = vb_reduce_cortex(brain_file_tmp,1:size(V,1),reduce); 0097 I = length(ix_act); 0098 0099 % Calculate mean distance between reduced vertices to determine 0100 % smoothing filter size 0101 R = 0.0; 0102 N = 0; 0103 [nextDD,nextIX] = vb_load_cortex_neighbour(brain_file_tmp); 0104 for i=1:length(ix_act) 0105 [tmp,ixx] = union(nextIX{ix_act(i)},ix_act); 0106 R = R+sum(nextDD{ix_act(i)}(ixx)); 0107 N = N+length(ixx); 0108 end 0109 R = R/N; 0110 0111 % Calculate smoothing filter. It will be used for displaying purpose, 0112 % not affect further analysis. 0113 [W,ix_act_ex] ... 0114 = vb_spatial_gauss_filter(brain_file_tmp,R,2*R,ix_act,1,-1,1); 0115 W = W./repmat(sum(W,2),[1 size(W,2)]); 0116 0117 % Load average current and map onto the template cortical surface 0118 C = vb_load_cortex_pmat(brain_file{1},key{1}); 0119 I = size(C,2); 0120 T = size(J,2); 0121 Jtmp = zeros(I,T); 0122 0123 Jtmp(Jinfo.ix_act_ex,:) = J; 0124 Jtmp = C*Jtmp; % current on the template brain 0125 Zavr = Jtmp(ix_act,:); 0126 end 0127 0128 % 0129 % Load subsequent files and calculate intersubject average 0130 % 0131 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0132 for i=2:length(curr_file) 0133 [Jinfo,J] = vb_load_current(curr_file{i},0,true); 0134 C = vb_load_cortex_pmat(brain_file{i},key{i}); 0135 0136 % added by taku-y 2011-11-15 0137 Jtmp = zeros(size(C,2),size(J,2)); 0138 Jtmp(Jinfo.ix_act_ex,:) = J; 0139 J = Jtmp; 0140 0141 if reduce==1.0 0142 Javr = Javr+C*J; 0143 else 0144 % Map current onto the template cortical surface 0145 I = size(C,2); 0146 Jtmp = zeros(I,T); 0147 0148 Jtmp(Jinfo.ix_act_ex,:) = J; 0149 Jtmp = C*Jtmp; % current on the template brain 0150 Zavr = Zavr+Jtmp(ix_act,:); 0151 end 0152 end 0153 0154 % 0155 % Variables to be saved 0156 % 0157 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0158 Jinfo = vb_load_current(curr_file{1},0,true); 0159 Jinfo.Ntrial = 1; 0160 0161 if reduce==1.0 0162 NJact = size(C,1); 0163 Jinfo.ix_act = 1:NJact; 0164 Jinfo.ix_act_ex = 1:NJact; 0165 Jinfo.NJact = NJact; 0166 Jinfo.Wact = speye(NJact); 0167 0168 load(curr_file{1},'bayes_parm'); 0169 Jinfo.patch_norm = bayes_parm.patch_norm; 0170 0171 Zact = Javr./length(curr_file); 0172 else 0173 Jinfo.ix_act = ix_act; 0174 Jinfo.ix_act_ex = ix_act_ex; 0175 Jinfo.NJact = length(ix_act_ex); 0176 Jinfo.Wact = W; 0177 0178 load(curr_file{1},'bayes_parm'); 0179 Jinfo.patch_norm = bayes_parm.patch_norm; 0180 0181 Zact = Zavr./length(curr_file); 0182 end 0183 0184 % 0185 % Save result 0186 % 0187 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0188 vb_fsave(curr_file_out,'Zact','Jinfo'); 0189 0190 % 0191 % Save project file 0192 % 0193 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0194 proj_file = get_project_filename; 0195 if isempty(proj_file), 0196 return; 0197 end; 0198 0199 project_file_mgr('load', proj_file); 0200 project_file_mgr('add', 'avr_parm', avr_parm); 0201 0202 return;