job_plot_cmap_load_current

 Load cortical current or time-frequency data. 

 [syntax]
 job_plot_cmap_load_current(inst_id,currfile,ix_trial,verbose)

 [history]
 2008-09-01 Taku Yoshioka
 2009-07-30 Taku Yoshioka
  Multiple trial current data has been supported. 
 2010-03-08 rhayashi
  Read single trial data at a one time
  when multiple trial current data is inputted
 2010-12-01 taku-y
  [internal change] 'Jinfo' replaces 'MEGinfo'.
 2010-12-06 taku-y
  [internal change] 'Jinfo' of version 0.9 or later preferentially
  used. 
 2010-12-08 taku-y
  [internal change] vb_version_cmp used.
 2010-12-11 taku-y
  [enhancement] Time-frequency data supported. The filetype is
  determined by file extention ('.curr.mat' or '.tf.mat'). 
 2011-01-29 taku-y
  [enhancement] Activity map plot supported. 
 2011-02-02 taku-y
  [minor] When loading cortical current it can be averaged (as do
  job_plot_meg). 
 2011-02-16 taku-y
  [enhancement] Cortical area plot supported. 
 2012-01-26 sako
  [minor change] n_trial -> "n_trial" was replaced to "ix_trial"
 2012-02-09 taku-y
  [debug] Support for .curr.mat file including both Jinfo and MEGinfo
  (old format)

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function job_plot_cmap_load_current(inst_id,currfile,ix_trial,verbose)
0002 % Load cortical current or time-frequency data.
0003 %
0004 % [syntax]
0005 % job_plot_cmap_load_current(inst_id,currfile,ix_trial,verbose)
0006 %
0007 % [history]
0008 % 2008-09-01 Taku Yoshioka
0009 % 2009-07-30 Taku Yoshioka
0010 %  Multiple trial current data has been supported.
0011 % 2010-03-08 rhayashi
0012 %  Read single trial data at a one time
0013 %  when multiple trial current data is inputted
0014 % 2010-12-01 taku-y
0015 %  [internal change] 'Jinfo' replaces 'MEGinfo'.
0016 % 2010-12-06 taku-y
0017 %  [internal change] 'Jinfo' of version 0.9 or later preferentially
0018 %  used.
0019 % 2010-12-08 taku-y
0020 %  [internal change] vb_version_cmp used.
0021 % 2010-12-11 taku-y
0022 %  [enhancement] Time-frequency data supported. The filetype is
0023 %  determined by file extention ('.curr.mat' or '.tf.mat').
0024 % 2011-01-29 taku-y
0025 %  [enhancement] Activity map plot supported.
0026 % 2011-02-02 taku-y
0027 %  [minor] When loading cortical current it can be averaged (as do
0028 %  job_plot_meg).
0029 % 2011-02-16 taku-y
0030 %  [enhancement] Cortical area plot supported.
0031 % 2012-01-26 sako
0032 %  [minor change] n_trial -> "n_trial" was replaced to "ix_trial"
0033 % 2012-02-09 taku-y
0034 %  [debug] Support for .curr.mat file including both Jinfo and MEGinfo
0035 %  (old format)
0036 %
0037 % Copyright (C) 2011, ATR All Rights Reserved.
0038 % License : New BSD License(see VBMEG_LICENSE.txt)
0039 
0040 global vbmeg_inst; 
0041 
0042 error('This function will become obsolete. (Feb 20 2012 taku-y)')
0043 
0044 %
0045 % Input parameters
0046 %
0047 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0048 if nargin<4, verbose = true; end
0049 if nargin<3, ix_trial= 1; end
0050 
0051 %
0052 % Global variables
0053 %
0054 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0055 h_fig = vbmeg_inst.plotcmap{inst_id}.h_fig;
0056 ed11  = vbmeg_inst.plotcmap{inst_id}.ed11;
0057 ed12  = vbmeg_inst.plotcmap{inst_id}.ed12;
0058 
0059 set(h_fig,'Pointer','watch');
0060 drawnow;
0061 
0062 %
0063 % Load estimated current
0064 %
0065 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0066 if ~exist(currfile,'file'),
0067   warndlg(['File ' currfile ' does not exist.']);
0068   return;
0069 end
0070 
0071 drawnow;
0072 
0073 if findstr(currfile,'.curr.mat'),  % Load cortical current file
0074   % Clear time-frequency variables
0075   vbmeg_inst.plotcmap{inst_id}.tffile     = [];
0076   vbmeg_inst.plotcmap{inst_id}.data       = [];
0077   vbmeg_inst.plotcmap{inst_id}.TFmax      = [];
0078   vbmeg_inst.plotcmap{inst_id}.TFmin      = [];
0079   vbmeg_inst.plotcmap{inst_id}.TFmax_plot = [];
0080   vbmeg_inst.plotcmap{inst_id}.TFmin_plot = [];
0081   vbmeg_inst.plotcmap{inst_id}.foi        = [];
0082   vbmeg_inst.plotcmap{inst_id}.TFinfo     = [];
0083   
0084   % Filename
0085   vbmeg_inst.plotcmap{inst_id}.currfile = currfile;
0086   
0087   % If cortical current file has multiple trials, ask if you load
0088   % averaged data or not.
0089   Jinfo = vb_load_current_info(currfile);
0090   if Jinfo.Ntrial>1, 
0091     button = questdlg(['The specified cortical current file has ' ...
0092                        'multiple trials. Do you want to load' ...
0093                        'trial-average data or the original one?'], ...
0094                       'Load cortical current data','Average','Raw', ...
0095                       'Average');
0096     if strcmp(button,'Average'), 
0097       ave_mode = true;
0098     else
0099       ave_mode = false;
0100     end
0101   else
0102     ave_mode = true;
0103   end
0104   vbmeg_inst.plotcmap{inst_id}.ave_mode = ave_mode;
0105 
0106   % load current data
0107   [tmp,Jact] = vb_load_current(currfile,0,ave_mode,ix_trial);
0108   if ~isempty(vbmeg_inst.plotcmap{inst_id}.V), 
0109     if max(Jinfo.ix_act_ex)>size(vbmeg_inst.plotcmap{inst_id}.V,1), 
0110       warndlg(['Index of dipole currents exceeds the number of ' ...
0111                'vertices.']);
0112       return;
0113     end
0114   end
0115 
0116   % Dipole current direction is flipped to inward
0117   vbmeg_inst.plotcmap{inst_id}.Jact = -1*Jact;
0118   clear Jact;
0119   
0120   % plotcmap.ix_act was obsoleted
0121   if ~isfield(Jinfo,'version') || ...
0122         vb_version_cmp(Jinfo.version,'<','0.9-0.a.0'),
0123     if isfield(Jinfo,'ix_act_ex'), 
0124       vbmeg_inst.plotcmap{inst_id}.ix_act = Jinfo.ix_act_ex;
0125     else
0126       vbmeg_inst.plotcmap{inst_id}.ix_act = Jinfo.ix_act;
0127     end
0128   end
0129 
0130   % Background current is not supported.
0131   vbmeg_inst.plotcmap{inst_id}.Jbck = [];
0132   vbmeg_inst.plotcmap{inst_id}.ix_bck = [];
0133 
0134   % Temporal information
0135   load(currfile,'bayes_parm','MEGinfo');
0136   
0137   if isfield(Jinfo,'version') & ...
0138         vb_version_cmp(Jinfo.version,'>=','0.9-0.a.0'), 
0139     T = Jinfo.Tmsec;
0140     vbmeg_inst.plotcmap{inst_id}.Jinfo = Jinfo;
0141   elseif exist('MEGinfo'), 
0142     T = (1000/MEGinfo.SampleFreq) ...
0143         *(Jinfo.Tsample-vb_meginfo_get_pre_trigger(MEGinfo));
0144     if exist('Jinfo'), 
0145       vbmeg_inst.plotcmap{inst_id}.Jinfo = Jinfo;
0146     else
0147       vbmeg_inst.plotcmap{inst_id}.Jinfo = MEGinfo;
0148     end
0149     %vbmeg_inst.plotcmap{inst_id}.Tmsec = T;
0150     vbmeg_inst.plotcmap{inst_id}.Jinfo.version = '0.7';
0151     vbmeg_inst.plotcmap{inst_id}.Tstart = Jinfo.Tsample(1);
0152   else
0153     T = Jinfo.Tmsec;
0154     %vbmeg_inst.plotcmap{inst_id}.Tmsec = T;
0155     vbmeg_inst.plotcmap{inst_id}.Jinfo = Jinfo;
0156     vbmeg_inst.plotcmap{inst_id}.Jinfo.version = '0.7';
0157     vbmeg_inst.plotcmap{inst_id}.Tstart = Jinfo.Tsample(1);
0158   end
0159 
0160   vbmeg_inst.plotcmap{inst_id}.toi = [T(1) T(2)];
0161 
0162   % Other parameters
0163   if isfield(Jinfo,'patch_norm'), 
0164     patch_norm = Jinfo.patch_norm;
0165   elseif isfield(Jinfo,'bayes_parm') ...
0166         & isfield(Jinfo.bayes_parm,'patch_norm'), 
0167     patch_norm = Jinfo.bayes_parm.patch_norm
0168     vbmeg_inst.plotcmap{inst_id}.Jinfo.patch_norm ...
0169         = Jinfo.bayes_parm.patch_norm;
0170   elseif exist('bayes_parm'), 
0171     % following line will be deleted after checking function:
0172     % vbmeg_inst.plotcmap{inst_id}.bayes_parm = bayes_parm;
0173     vbmeg_inst.plotcmap{inst_id}.Jinfo.patch_norm = ...
0174         bayes_parm.patch_norm;
0175     patch_norm = bayes_parm.patch_norm;
0176   else
0177     patch_norm = false;
0178   end
0179 
0180   if ~isempty(ix_trial)
0181     vbmeg_inst.plotcmap{inst_id}.ix_trial = ix_trial;
0182   else
0183     vbmeg_inst.plotcmap{inst_id}.ix_trial = 1;
0184   end
0185   set(vbmeg_inst.plotcmap{inst_id}.ed1,'String', ...
0186     num2str(vbmeg_inst.plotcmap{inst_id}.ix_trial));
0187 
0188   % Physical unit of current
0189   if patch_norm, % [pAm/mm^2]
0190     vbmeg_inst.plotcmap{inst_id}.Jact ...
0191         = vbmeg_inst.plotcmap{inst_id}.Jact*1e6; % (=1e12*1e-6);
0192     vbmeg_inst.plotcmap{inst_id}.Jbck ...
0193         = vbmeg_inst.plotcmap{inst_id}.Jbck*1e6; % (=1e12*1e-6);
0194     %set(vbmeg_inst.plotcmap{inst_id}.Jmin_unit_text,'String', ...
0195     %                  '[pAm/mm^2]');
0196     %set(vbmeg_inst.plotcmap{inst_id}.Jmax_unit_text,'String', ...
0197     %                  '[pAm/mm^2]');
0198   else                          % [nAm]
0199     vbmeg_inst.plotcmap{inst_id}.Jact ...
0200         = vbmeg_inst.plotcmap{inst_id}.Jact*1e9;
0201     vbmeg_inst.plotcmap{inst_id}.Jbck ...
0202         = vbmeg_inst.plotcmap{inst_id}.Jbck*1e9;
0203     %set(vbmeg_inst.plotcmap{inst_id}.Jmin_unit_text,'String', ...
0204     %                  '[nAm]');
0205     %set(vbmeg_inst.plotcmap{inst_id}.Jmax_unit_text,'String', ...
0206     %                  '[nAm]');
0207   end
0208 
0209   % Scale of current amplitude
0210   if size(vbmeg_inst.plotcmap{inst_id}.Jact,3) > 1
0211     ix_trial = vbmeg_inst.plotcmap{inst_id}.ix_trial;
0212     Jmax = max(max(abs(vbmeg_inst.plotcmap{inst_id}.Jact(:,:,ix_trial))));
0213     Jmin = min(min(abs(vbmeg_inst.plotcmap{inst_id}.Jact(:,:,ix_trial))));
0214   else
0215     Jmax = max(max(abs(vbmeg_inst.plotcmap{inst_id}.Jact(:,:))));
0216     Jmin = min(min(abs(vbmeg_inst.plotcmap{inst_id}.Jact(:,:))));
0217   end
0218   vbmeg_inst.plotcmap{inst_id}.Jmax = Jmax;
0219   vbmeg_inst.plotcmap{inst_id}.Jmin = Jmin;
0220   vbmeg_inst.plotcmap{inst_id}.Jmax_plot = Jmax;
0221   vbmeg_inst.plotcmap{inst_id}.Jmin_plot = Jmax*0.1;
0222   %set(vbmeg_inst.plotcmap{inst_id}.Jmax_plot_edit,'String',...
0223   %                  num2str(Jmax));
0224   %set(vbmeg_inst.plotcmap{inst_id}.Jmin_plot_edit,'String',...
0225   %                  num2str(0.1*Jmax));
0226   set(vbmeg_inst.plotcmap{inst_id}.ed10,'String','');
0227 elseif findstr(currfile,'.tf.mat'), % Load time frequency file
0228   % Clear cortical current variables
0229   vbmeg_inst.plotcmap{inst_id}.currfile  = [];
0230   vbmeg_inst.plotcmap{inst_id}.data      = [];
0231   vbmeg_inst.plotcmap{inst_id}.Jmax      = [];
0232   vbmeg_inst.plotcmap{inst_id}.Jmin      = [];
0233   vbmeg_inst.plotcmap{inst_id}.Jmax_plot = [];
0234   vbmeg_inst.plotcmap{inst_id}.Jmin_plot = [];
0235   vbmeg_inst.plotcmap{inst_id}.foi       = [];
0236   vbmeg_inst.plotcmap{inst_id}.Jinfo     = [];
0237   
0238   % Filename
0239   vbmeg_inst.plotcmap{inst_id}.tffile = currfile;
0240   
0241   % Time-frequency data
0242   load(currfile,'data','TFinfo');
0243   if ~isfield(TFinfo,'freq_scale'), % TFinfo format under development
0244     TFinfo.freq_scale = 'log';
0245   end
0246   vbmeg_inst.plotcmap{inst_id}.data   = data;
0247   vbmeg_inst.plotcmap{inst_id}.TFinfo = TFinfo;
0248   clear data;
0249   
0250   % Plotting threshold
0251   TFmax = max(max(max(abs(vbmeg_inst.plotcmap{inst_id}.data))));
0252   TFmin = min(min(min(vbmeg_inst.plotcmap{inst_id}.data)));
0253   vbmeg_inst.plotcmap{inst_id}.TFmax      = TFmax;
0254   vbmeg_inst.plotcmap{inst_id}.TFmin      = TFmin;
0255   vbmeg_inst.plotcmap{inst_id}.TFmax_plot = TFmax;
0256   vbmeg_inst.plotcmap{inst_id}.TFmin_plot = 0.1*TFmax;
0257   %set(vbmeg_inst.plotcmap{inst_id}.Jmax_plot_edit,'String',...
0258   %                  num2str(TFmax));
0259   %set(vbmeg_inst.plotcmap{inst_id}.Jmin_plot_edit,'String',...
0260   %                  num2str(0.1*TFmax));
0261   set(vbmeg_inst.plotcmap{inst_id}.ed10,'String',...
0262                     num2str(TFmax));
0263   
0264   % Time-frequency information
0265   T = TFinfo.Tmsec;
0266   F = TFinfo.freq;
0267   vbmeg_inst.plotcmap{inst_id}.toi = [T(1) T(2)];
0268   vbmeg_inst.plotcmap{inst_id}.foi = [F(1) F(2)];
0269 elseif findstr(currfile,'.act.mat'), % Load activity map file
0270   % Filename
0271   vbmeg_inst.plotcmap{inst_id}.act_file = currfile;
0272   
0273   % Update listbox
0274   keyset = vb_get_keyset_act(currfile);
0275   set(vbmeg_inst.plotcmap{inst_id}.lb1,'String',keyset);
0276   set(vbmeg_inst.plotcmap{inst_id}.lb1,'Value',1);
0277   
0278   % Update activity map
0279   job_plot_cmap_update_xxP(inst_id);
0280 else % Load cortical area file
0281   % Filename
0282   vbmeg_inst.plotcmap{inst_id}.area_file = currfile;
0283   
0284   % Update listbox
0285   keyset = vb_get_keyset_area(currfile);
0286   set(vbmeg_inst.plotcmap{inst_id}.lb2,'String',keyset);
0287   set(vbmeg_inst.plotcmap{inst_id}.lb2,'Value',1);
0288   
0289   % Update cortical area
0290   job_plot_cmap_update_area(inst_id);
0291 end
0292 
0293 % Update edittext of color scaling
0294 job_plot_cmap_update_colorbarscale(inst_id);
0295 
0296 %
0297 % Change figure title/act_file/area_file edit box
0298 %
0299 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0300 [tmp,fname,ext] = fileparts(currfile);
0301 if ~isempty(findstr(currfile,'.curr.mat')) ...
0302       || ~isempty(findstr(currfile,'.tf.mat')), 
0303   set(h_fig,'Name',['job_plot_currentmap: ' fname ext]);
0304 elseif findstr(currfile,'.act.mat'), 
0305   set(ed11,'String',[fname ext]);
0306 else
0307   set(ed12,'String',[fname ext]);
0308 end
0309 
0310 %
0311 % Load cortical surface model
0312 %
0313 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0314 if exist('bayes_parm'), 
0315   brainfile = bayes_parm.brainfile;
0316 else
0317   brainfile = '';
0318 end
0319 
0320 % Check if a cortical surface model has been loaded
0321 brainfile_old = vbmeg_inst.plotcmap{inst_id}.brainfile; 
0322 if ~isempty(brainfile_old) & ~isempty(brainfile), 
0323   if verbose == true
0324     qstring = ['Cortical surface model ''' brainfile_old ' has already ' ...
0325                'been loaded. Do you replace this model by ' ...
0326                brainfile '?'];
0327     button = questdlg(qstring,'Overwrite cortical surface model',...
0328                         'Yes','No','Yes');
0329   else
0330     button = 'No';
0331   end
0332   
0333   if strcmp(button,'No'), 
0334     job_plot_cmap_update_filterradius(inst_id);
0335     job_plot_cmap_update_timecourse(inst_id); 
0336     job_plot_cmap_update_spatialpattern(inst_id);
0337     return;
0338   end
0339 end
0340     
0341 if exist(brainfile,'file'),   
0342   vbmeg_inst.plotcmap{inst_id}.brainfile = brainfile;
0343   
0344   % Cortical surface model
0345   [Vinf,F,xx,inf_C,xxA] = vb_load_cortex(brainfile,'Inflate');
0346   V = vb_load_cortex(brainfile);
0347   vbmeg_inst.plotcmap{inst_id}.V = V;
0348   vbmeg_inst.plotcmap{inst_id}.Vinf = Vinf;
0349   vbmeg_inst.plotcmap{inst_id}.F = F;
0350   vbmeg_inst.plotcmap{inst_id}.inf_C = inf_C;
0351   vbmeg_inst.plotcmap{inst_id}.xxA = xxA;
0352 
0353   % Standard brain coordinate
0354   load(brainfile,'Vtal','Vmni');
0355   if ~exist('Vtal','var'), Vtal = zeros(size(V)); end
0356   if ~exist('Vmni','var'), Vmni = zeros(size(V)); end
0357   vbmeg_inst.plotcmap{inst_id}.Vtal = Vtal*1e3;
0358   vbmeg_inst.plotcmap{inst_id}.Vmni = Vmni*1e3;
0359 
0360   % Load neighbour data for filter calculation
0361   [nextDD,nextIX] = vb_load_cortex_neighbour(brainfile); 
0362   vbmeg_inst.plotcmap{inst_id}.nextDD = nextDD;
0363   vbmeg_inst.plotcmap{inst_id}.nextIX = nextIX;
0364 end;
0365 %elseif verbose,
0366 %  warndlg(['File ''' brainfile ''' does not exist. You have to ' ...
0367 %           'manually choose the cortical surface file ' ...
0368 %           '(.brain.mat).']);
0369 %  return;
0370 %end
0371 
0372 %
0373 % Update plots
0374 %
0375 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0376 job_plot_cmap_update_filterradius(inst_id);
0377 job_plot_cmap_update_timecourse(inst_id); 
0378 job_plot_cmap_update_spatialpattern(inst_id);
0379 set(h_fig,'Pointer','arrow');