show_source_current_meg

 Show source current estimated from MEG data

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

0001 function show_source_current_meg(p)
0002 % Show source current estimated from MEG data
0003 %
0004 % Copyright (C) 2011, ATR All Rights Reserved.
0005 % License : New BSD License(see VBMEG_LICENSE.txt)
0006 
0007 disp(mfilename);
0008 
0009 % Set parameters
0010 time_to_show = [-0.1 0.4];% sec
0011 time_of_interest = [0 0.24];% sec
0012 angle(:,:,1) = [-90 0;90 0];% Angle to see brain (auditory)
0013 angle(:,:,2) = [-90 0;90 0];% Angle to see brain (sensorimotor)
0014 angle(:,:,3) = [-50 0;50 0];% Angle to see brain (visual)
0015 
0016 % Make directory to save figures
0017 save_dir = fullfile(p.fig_root, mfilename, p.sub);
0018 if exist(save_dir, 'dir') ~= 7
0019     vb_mkdir(save_dir);
0020 end
0021 
0022 % Load time information
0023 current_file = fullfile(p.proj_root, p.meg_dirname, p.current_dirname, [p.task_list{1} '.curr.mat']);
0024 Jinfo = vb_load_current(current_file);
0025 time = Jinfo.Tmsec/1000;% sec
0026 [~, from_show] = min(abs(time-time_to_show(1)));
0027 [~, to_show] = min(abs(time-time_to_show(2)));
0028 time = time(from_show: to_show);
0029 
0030 % Load brain model on which current is shown
0031 brain_file = fullfile(p.proj_root, p.brain_dirname, [p.struct_name, '.brain.mat']);
0032 [V, F, xx, inf_C] = vb_load_cortex(brain_file, 'Inflate');
0033 plot_parm = vb_set_plot_parm;
0034 
0035 for ta = 1:length(p.task_list)
0036     
0037     % Load current
0038     current_file = fullfile(p.proj_root, p.meg_dirname, p.current_dirname, [p.task_list{ta} '.curr.mat']);
0039     [~, current] = vb_load_current(current_file);
0040     
0041     % Average source currents across trials
0042     y = mean(current(:, from_show:to_show, :), 3);
0043     
0044     % Calculate absolute value of current
0045     [~, from_power] = min(abs(time-time_of_interest(1)));
0046     [~, to_power] = min(abs(time-time_of_interest(2)));
0047     ay=mean(abs(y(:,from_power:to_power)),2);
0048     
0049     % Plot time series of averaged MEG data
0050     h = figure;
0051     max_y = max(y(:));
0052     min_y = min(y(:));
0053     subplot(2, 1, 1)
0054     plot(time, y)
0055     xlabel('Time [s]')
0056     ylabel('Current [Am/m^2]')
0057     axis([time(1) time(end) min_y max_y])
0058     
0059     % Plot spatial pattern of absolute value of current
0060     for an = 1:size(angle, 1)
0061         subplot(2, 2, 2+an)
0062         vb_plot_cortex(plot_parm, V, F, inf_C, ay, max(ay)*0.3, max(ay));
0063         axis equal off
0064         view(angle(an, :, ta))
0065         %zoom(1.5)
0066         if an==1
0067             title(['Absolute value of current (' num2str(time_of_interest(1)) ' - ' num2str(time_of_interest(2)) ' sec)'])
0068         end
0069     end
0070     
0071     % Save figure
0072     vb_savefig_as_shown(h, fullfile(save_dir, p.task_list{ta}))
0073 end
0074 
0075