convert current unit(density/amplitude)
[syntax]
J_out = vb_convert_current_unit(proj_root, ...
current_file, to_density...
[,new_current_file]);
[input]
proj_root : VBMEG result file directory .
current_file : Cortical current file (.curr.mat).
* relative path from proj_root
to_density : unit after conversion
= true : convert to density (Am/m^2).
= false : convert to amplitude (Am).
new_current_file : [Optional] save current file.
* relative path from proj_root
[output]
J_out :
to_density = true
J_out: current density for 'J-current' (Am/m^2).
to_density = false
J_out: current amplitude for 'J-current' (Am).
--- 'J-current' is always returned
-----------------------------------------------------------
--- Smoothing leadfield
G = basis*S*W (S: area of vertex, W: smoothing filter)
--- MEG - current relation
MEG = G * Zact; % Nsensor x TimeSample
= basis*S*Wact*Zact
= basis*S*Jact_density (patch_norm: ON)
= basis*Jact_amplitude (patch_norm: OFF)
--- conversion between J_amplitude and J_density
Jact_amplitude = Jact_density.*S
Jact_density = Jact_amplitude./S
--- S should be applied to Jact
[history]
2014-10-2 Masa-aki Sato
[note]
It should be noted that in general converted current density may be
different from the current density obtained by setting patch_norm = 1
using vb_job_vb(current variance estimation) because of inversion
process.
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function [J_out] = vb_convert_current_unit(proj_root, ... 0002 current_file, ... 0003 to_density, new_current_file) 0004 % convert current unit(density/amplitude) 0005 % 0006 % [syntax] 0007 % J_out = vb_convert_current_unit(proj_root, ... 0008 % current_file, to_density... 0009 % [,new_current_file]); 0010 % 0011 % [input] 0012 % proj_root : VBMEG result file directory . 0013 % current_file : Cortical current file (.curr.mat). 0014 % * relative path from proj_root 0015 % to_density : unit after conversion 0016 % = true : convert to density (Am/m^2). 0017 % = false : convert to amplitude (Am). 0018 % new_current_file : [Optional] save current file. 0019 % * relative path from proj_root 0020 % 0021 % [output] 0022 % J_out : 0023 % to_density = true 0024 % J_out: current density for 'J-current' (Am/m^2). 0025 % to_density = false 0026 % J_out: current amplitude for 'J-current' (Am). 0027 % 0028 %--- 'J-current' is always returned 0029 %----------------------------------------------------------- 0030 %--- Smoothing leadfield 0031 % G = basis*S*W (S: area of vertex, W: smoothing filter) 0032 %--- MEG - current relation 0033 % MEG = G * Zact; % Nsensor x TimeSample 0034 % = basis*S*Wact*Zact 0035 % = basis*S*Jact_density (patch_norm: ON) 0036 % = basis*Jact_amplitude (patch_norm: OFF) 0037 %--- conversion between J_amplitude and J_density 0038 % Jact_amplitude = Jact_density.*S 0039 % Jact_density = Jact_amplitude./S 0040 %--- S should be applied to Jact 0041 % 0042 % [history] 0043 % 2014-10-2 Masa-aki Sato 0044 % 0045 % [note] 0046 % It should be noted that in general converted current density may be 0047 % different from the current density obtained by setting patch_norm = 1 0048 % using vb_job_vb(current variance estimation) because of inversion 0049 % process. 0050 % 0051 % Copyright (C) 2011, ATR All Rights Reserved. 0052 % License : New BSD License(see VBMEG_LICENSE.txt) 0053 0054 % 0055 % --- Previous check 0056 % 0057 if nargin < 3 0058 error('Please check function usage.'); 0059 end 0060 if ~exist('new_current_file', 'var') 0061 new_current_file = ''; 0062 end 0063 save_flag = true; 0064 if isempty(new_current_file) 0065 save_flag = false; 0066 end 0067 % 0068 % --- Main Procedure 0069 % 0070 fprintf('starting unit conversion...'); 0071 current_file = fullfile(proj_root, current_file); 0072 0073 [Jinfo,Jact] = vb_load_current(current_file); 0074 % Jinfo.patch_norm: ON/OFF (Am/m^2) / (Am). 0075 % Jinfo.curr_type : 0: 'J-current' or 1: 'Z-current' 0076 % Jinfo.Wact : Smoothing Gaussian filter 0077 % Jact = Wact * Zact 0078 % Jinfo.ix_act : Vertex indices of Z-current. 0079 % Jinfo.ix_act_ex : Vertex indices of J-current. 0080 0081 0082 % 0083 % --- Create J_out 0084 % 0085 switch to_density 0086 case true 0087 % J_out: current density for 'J-current' (Am/m^2). 0088 if Jinfo.patch_norm==ON 0089 J_out = Jact; 0090 fprintf('\nNo need to convert unit because the specified file contains data in density.\n'); 0091 if save_flag 0092 save_flag = false; 0093 end 0094 else 0095 [Nact,Nsample,Ntrials] = size(Jact); 0096 % Jact : current timecourse 0097 % Jact(Nact,Nsample,Ntrials) 0098 % Nact = Lact * Nvact, Nvact = length(ix_act) 0099 % Nsample : # of time sample, 0100 % Ntrials : # of trials in all session] 0101 0102 ix_act_ex = Jinfo.ix_act_ex; 0103 Lact = Nact/length(ix_act_ex); 0104 0105 load(current_file,'bayes_parm'); 0106 0107 brainfile = fullfile(proj_root, bayes_parm.brainfile); 0108 0109 % Patch norm normalization 0110 xxA = vb_load_area_density(brainfile); 0111 0112 % Jact_density = Jact_amplitude./S 0113 S = repmat(xxA(ix_act_ex),[Lact, Nsample, Ntrials]); 0114 J_out = Jact./S; 0115 end 0116 0117 patch_norm = ON; 0118 0119 case false 0120 % J_out: current amplitude for 'J-current' (Am). 0121 if Jinfo.patch_norm==OFF 0122 J_out = Jact; 0123 fprintf('\nNo need to convert unit because the specified file contains data in amplitude.\n'); 0124 if save_flag 0125 save_flag = false; 0126 end 0127 else 0128 [Nact,Nsample,Ntrials] = size(Jact); 0129 % Jact : current timecourse 0130 % Jact(Nact,Nsample,Ntrials) 0131 % Nact = Lact * Nvact, Nvact = length(ix_act) 0132 % Nsample : # of time sample, 0133 % Ntrials : # of trials in all session] 0134 0135 ix_act_ex = Jinfo.ix_act_ex; 0136 Lact = Nact/length(ix_act_ex); 0137 0138 load(current_file,'bayes_parm'); 0139 0140 brainfile = fullfile(proj_root, bayes_parm.brainfile); 0141 0142 % Patch norm normalization 0143 xxA = vb_load_area_density(brainfile); 0144 0145 % Jact_amplitude = Jact_density.*S 0146 S = repmat(xxA(ix_act_ex),[Lact, Nsample, Ntrials]); 0147 J_out = Jact.*S; 0148 0149 patch_norm = OFF; 0150 end 0151 end 0152 disp('done.'); 0153 0154 if save_flag 0155 % 0156 % --- new current directory file creation 0157 % 0158 new_file = fullfile(proj_root, new_current_file); 0159 inner_new_current_file_dir_create(new_file); 0160 s = load(current_file); 0161 0162 0163 % Header modification 0164 s.Jinfo.curr_type = 'J-current'; 0165 s.Jinfo.patch_norm = patch_norm; 0166 0167 ix_act = s.Jinfo.ix_act_ex; 0168 0169 s.Jinfo = rmfield(s.Jinfo, 'ix_act_ex'); 0170 s.Jinfo.ix_act = ix_act; 0171 s.Jinfo.NJact = length(ix_act); 0172 s.Jinfo.Wact = []; 0173 0174 % Current timecourse 0175 if isfield(s, 'Zact') 0176 s = rmfield(s, 'Zact'); 0177 end 0178 s.Jact = J_out; 0179 0180 % save 0181 disp('Now saving...'); 0182 vb_save_struct(new_file, s); 0183 disp('done.'); 0184 end 0185 0186 if nargout 0187 if to_density 0188 unit = 'density'; 0189 else 0190 unit = 'amplitude'; 0191 end 0192 disp(['returning value: unit in ' unit, '.']); 0193 end 0194 0195 function inner_new_current_file_dir_create(new_current_file) 0196 0197 [p] = vb_get_file_parts(new_current_file); 0198 if exist(p, 'dir') ~= 7 0199 vb_mkdir(p); 0200 end