correct_classify_ic_meg

 Manually correct classification of ICs

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

0001 function correct_classify_ic_meg(p)
0002 % Manually correct classification of ICs
0003 %
0004 % Copyright (C) 2011, ATR All Rights Reserved.
0005 % License : New BSD License(see VBMEG_LICENSE.txt)
0006 
0007 disp(mfilename);
0008 
0009 % Classification program
0010 classify_ic_prog = 'classify_ic_meg';
0011 
0012 % Load MEG data
0013 data_file = fullfile(p.proj_root, p.meg_dirname, p.trial_dirname, 'all.info.mat');
0014 [meg, ~, time_info] = vb_load_meg_data(data_file);
0015 [Nch, Nt, Ntr] = size(meg);
0016 samplingrate = time_info.sample_frequency;
0017 time = time_info.time;
0018 
0019 % Load positions of channels
0020 [pos, channel_info] = vb_load_channel(data_file);
0021 ix_axial = find(channel_info.Type == 2);% axial sensor
0022 pos = pos(ix_axial, :);
0023 
0024 % Load ICA result
0025 ica_dir = fullfile(p.proj_root, p.meg_dirname, p.ica_dirname);
0026 ica_result = fullfile(ica_dir, [classify_ic_prog '.mat']);
0027 load(ica_result,'w','s','v','unmix','mix','brain_ic','noise_ic')
0028 
0029 % Calculate ICs
0030 ic = zeros(Nch, Nt, Ntr);
0031 for tr = 1:Ntr
0032     ic(:,:,tr) = unmix*meg(:,:,tr);  % [component x sample x trial]
0033 end
0034 
0035 % Calculate stimulu-triggered average and power spectrum
0036 mic = mean(ic, 3);% Average each IC across trials
0037 [pow, f] = vb_fftpower(ic, samplingrate);% Calculate power spectrum
0038 mp = mean(log(pow), 3);% Average power spectrum across trials
0039 
0040 % Obtain classification
0041 classification = -ones(1, Nch);
0042 classification(brain_ic) = 1;% 1: brain_ic, -1: noise_ic
0043 
0044 % Check and correct classification
0045 disp('Brain/noise ICs are represented by blue/red lines.')
0046 disp('Input IC number to correct or just push enter to go ahead.')
0047 disp('If you want to exit correction, input 0.')
0048 h = figure;
0049 cc = 1;
0050 while cc <= Nch
0051     clf(h, 'reset')
0052     set(0, 'CurrentFigure', h)
0053     k = 1;% Position of figure
0054     for c = cc:min([cc+4,Nch]);
0055         if classification(c) == 1
0056             color = 'b';
0057         else
0058             color = 'r';
0059         end
0060 
0061         % Plot spatial pattern of mixing vector
0062         subplot(5, 3, 1+3*(k-1))
0063         vb_plot_sensor_2d(pos, mix(ix_axial,c));
0064         vb_plot_sensor_2d_head_plot_add(gca);
0065         axis square off
0066         title(['IC ' num2str(c)])
0067         
0068         % Plot stimulus-triggered average of components
0069         subplot(5, 3, 2+3*(k-1))
0070         plot(time, mic(c, :), color);
0071         xlim([-0.5 1])
0072         if k == 1
0073             title('Stimulus-triggered average')
0074         elseif k == 5
0075             xlabel('Time [s]')
0076         end
0077         
0078         % Plot power spectrum
0079         subplot(5, 3, 3+3*(k-1))
0080         plot(f, mp(c, :), color);
0081         xlim([0 50])
0082         if k == 1
0083             title('Power spectrum (log)')
0084         elseif k == 5
0085             xlabel('Frequency [Hz]')
0086         end
0087         k = k+1;
0088     end
0089     ax1 = axes('Position', [0 0 1 1], 'Visible', 'off');
0090     text1 = {'Blue : Brain'; 'Red : Noise'};
0091     axes(ax1)
0092     text(0.93, 0.87, text1)
0093     drawnow
0094 
0095     % Obtain iput
0096     a = input('Input IC number to correct:');
0097     if isempty(a)% Go ahead
0098         cc = cc+5;
0099     elseif a == 0% Exit
0100         break
0101     else% Correct classification
0102         classification(a) = -classification(a);
0103     end
0104 end
0105 
0106 % Save result
0107 brain_ic = find(classification == 1);
0108 noise_ic = find(classification == -1);
0109 output_file = fullfile(ica_dir, [mfilename '.mat']);
0110 load(ica_result, 'w', 's', 'v', 'unmix', 'mix')
0111 save(output_file, 'w', 's', 'v', 'unmix', 'mix', 'brain_ic', 'noise_ic')
0112 
0113 
0114 
0115