classify_ic_meg

 Aucomatically classify ICs into brain activites and noise

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

0001 function classify_ic_meg(p)
0002 % Aucomatically classify ICs into brain activites and noise
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 method_list = {'kurtosis', 'entropy'};
0011 Nmethod = length(method_list);
0012 ica_program = 'apply_ica_meg';
0013 
0014 % Load MEG data
0015 data_file = fullfile(p.proj_root, p.meg_dirname, p.trial_dirname, 'all.info.mat');
0016 meg = vb_load_meg_data(data_file);
0017 [Nch, Nt, Ntr] = size(meg);
0018 
0019 % Load ICA result
0020 ica_dir = fullfile(p.proj_root, p.meg_dirname, p.ica_dirname);
0021 ica_result = fullfile(ica_dir, [ica_program '.mat']);
0022 load(ica_result)
0023 
0024 % Calculate ICs
0025 ic = zeros(Nch,Nt,Ntr);
0026 for tr = 1:Ntr
0027     ic(:, :, tr) = unmix*meg(:, :, tr);  % [component x sample x trial]
0028 end
0029 
0030 % Classify each IC based on the number of outliers
0031 brain_ic_each_method = cell(Nmethod, 1);
0032 for method = 1:Nmethod
0033     
0034     value = zeros(Nch,Ntr);
0035     for comp = 1:Nch
0036         for tr = 1:Ntr
0037             eval(['value(comp,tr)=' method_list{method} '(ic(comp,:,tr));']);
0038         end
0039     end
0040     mean_val = mean(value(:));% Mean value across ICs and trials
0041     std_val = std(value(:));% Standard deviation across ICs and trials
0042     
0043     fprintf('Noise ICs [%s]: \n', method_list{method});
0044     brain_ic_each_method{method} = [];
0045     for comp = 1:Nch
0046         % Count outlier trials
0047         Noutlier = 0;
0048         for tr = 1:Ntr
0049             value(comp, tr) = (value(comp,tr)-mean_val)/std_val;
0050             if abs(value(comp, tr)) > p.threshold_ourlier
0051                 Noutlier = Noutlier+1;
0052             end
0053         end
0054         
0055         % Classify IC based on the number of outliers
0056         if Noutlier/Ntr <= p.threshold_num_of_outlier
0057             brain_ic_each_method{method} = [brain_ic_each_method{method} comp];
0058         else
0059             fprintf('[%s]', num2str(comp));
0060         end
0061     end
0062     fprintf('\n');
0063 end
0064 
0065 % Determine brain and noise ICs based on classification
0066 brain_ic = brain_ic_each_method{1};
0067 for method = 2:Nmethod
0068     brain_ic = intersect(brain_ic,brain_ic_each_method{method});
0069 end
0070 noise_ic = 1:Nch;
0071 noise_ic(brain_ic) = [];
0072 
0073 % Save result
0074 output_file = fullfile(ica_dir, [mfilename '.mat']);
0075 load(ica_result, 'w', 's', 'v', 'unmix', 'mix')
0076 save(output_file,'w', 's', 'v', 'unmix', 'mix', 'brain_ic', 'noise_ic')
0077 
0078 
0079 % Inner function to calculate entropy
0080 function h = entropy(x, n)
0081 %ENTROPY Computes a first-order estimate of the entropy of a matrix.
0082 %   H = ENTROPY(X, N) returns the first-order estimate of matrix X
0083 %   with N symbols (N = 256 if omitted) in bits/symbol. The estimate
0084 %   assumes a statistically independent source characterized by the
0085 %   relative frequency of occurrence of the elements in X.
0086 
0087 %   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
0088 %   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
0089 %   $Revision: 1.4 $  $Date: 2003/10/26 18:35:35 $
0090 
0091 if nargin < 2   
0092    n = 256;                           % Default for n.
0093 end 
0094 
0095 x = double(x);                        % Make input double
0096 xh = hist(x(:), n);                   % Compute N-bin histogram
0097 xh = xh / sum(xh(:));                 % Compute probabilities
0098 
0099 % Make mask to eliminate 0's since log2(0) = -inf.
0100 i = find(xh);           
0101 
0102 h = -sum(xh(i) .* log2(xh(i)));       % Compute entropy
0103 
0104