1+ classdef eyeCatch
2+ % Detetcs eye ICs based solely based on on their scalpmaps.
3+ % The input can be an EEG structure, an array of scalpmaps (channwel weights and location, e.g.
4+ % from EEG.chanlocs and EEG.icawinv) or an scalpmap object from MPT.
5+ %
6+ % Example 1: (finding eye ICs in the EEG structure)
7+ %
8+ % >> eyeDetector = eyeCatch; % create an object from the class
9+ % >> [isEye similarity scalpmapObj] = eyeDetector.detectFromEEG(EEG); % detect eye ICs
10+ % >> find(isEye) % display the IC numbers for eye ICs (since isEye is a logical array)
11+ % >> scalpmapObj.plot(isEye) % plot eye ICs
12+ %
13+ % Example 2:
14+ %
15+ % >> [isEye similarity scalpmapObj] = eyeDetector.detectFromStudy(STUDY, ALLEEG); % read data from a loaded study
16+ % >> find(isEye) % display the IC numbers for eye ICs (since isEye is a logical array). The
17+ % % order of ICs is same order as in STUDY.cluster(1).comps .
18+ % >> scalpmapObj.plot(isEye) % plot eye ICs
19+ %
20+ % Written by Nima Bigdely-Shamlo, Swartz Center, INC, UCSD.
21+ % Copyright � 2012 University Of California San Diego. Distributed under BSD License.
22+
23+ properties
24+ eyeScalpmapDatabase
25+ eyeChannelWeightNormalized % precomputed normalized eye scalpmaps channel weights on interpolated 2D scalp.
26+ similarityThreshold = 0.94
27+ end ;
28+
29+ methods
30+ function obj = eyeCatch(similarityThreshold , varargin )
31+ % obj = eyeCatch(similarityThreshold)
32+ if nargin > 0
33+ obj.similarityThreshold = similarityThreshold ;
34+ end ;
35+
36+ load(' pooledEyeScalpmap.mat' % load pooledEyeScalpmap varaible
37+ load(' eyeChannelWeightNormalizedPart1.mat' % load eyeChannelWeightNormalizedPart1 variable
38+ load(' eyeChannelWeightNormalizedPart2.mat' % load eyeChannelWeightNormalizedPart1 variable
39+ eyeChannelWeightNormalized = cat(1 , eyeChannelWeightNormalizedPart1 , eyeChannelWeightNormalizedPart2 );
40+ clear eyeChannelWeightNormalizedPart1 eyeChannelWeightNormalizedPart2
41+
42+ obj.eyeScalpmapDatabase = pooledEyeScalpmap ;
43+ obj.eyeChannelWeightNormalized = eyeChannelWeightNormalized ;
44+ if size(obj .eyeChannelWeightNormalized ,1 ) ~= obj .eyeScalpmapDatabase .numberOfScalpmaps
45+ fprintf(' Precomputed weights are being recomputed.\n '
46+ % to do compute weights
47+ channelWeight = obj .eyeScalpmapDatabase .originalChannelWeight(: ,: );
48+ channelWeight(isnan(channelWeight )) = 0 ;
49+ channelWeight = channelWeight ' ;
50+
51+ channelWeightNormalized = bsxfun(@minus , channelWeight , mean(channelWeight ));
52+ channelWeightNormalized = bsxfun(@rdivide , channelWeightNormalized , std(channelWeightNormalized ));
53+ eyeChannelWeightNormalized = channelWeightNormalized ' ;
54+
55+ eyeChannelWeightNormalizedPart1 = eyeChannelWeightNormalized(1 : 2000 ,: );
56+ eyeChannelWeightNormalizedPart2 = eyeChannelWeightNormalized(2001 : end ,: );
57+
58+ % save the file (overwrite)
59+ pooledEyeScalpmap = obj .eyeScalpmapDatabase ;
60+ fullPath = which(' eyeChannelWeightNormalizedPart1.mat'
61+ try
62+ save(fullPath , ' eyeChannelWeightNormalizedPart1'
63+ fprintf(' New precomputed weights saved.\n '
64+
65+ fullPath = which(' eyeChannelWeightNormalizedPart2.mat'
66+ save(fullPath , ' eyeChannelWeightNormalizedPart2'
67+ catch
68+ fprintf(' Could not save newly precomputed weights, please check if you have write permission.\n '
69+ end ;
70+
71+
72+
73+ end ;
74+ end ;
75+
76+ function [isEye similarity ] = detectFromInterpolatedChannelWeight(obj , channelWeight , similarityThreshold , varargin )
77+ % [isEye similarity] = detectFromInterpolatedChannelWeight(obj, channelWeight, similarityThreshold)
78+
79+ if nargin < 3
80+ similarityThreshold = obj .similarityThreshold ;
81+ end ;
82+
83+ channelWeight(isnan(channelWeight )) = 0 ;
84+ channelWeight = channelWeight ' ;
85+
86+ channelWeightNormalized = bsxfun(@minus , channelWeight , mean(channelWeight ));
87+ channelWeightNormalized = bsxfun(@rdivide , channelWeightNormalized , std(channelWeightNormalized ));
88+ channelWeightNormalized = channelWeightNormalized ' ;
89+
90+ similarity = max(abs(obj .eyeChannelWeightNormalized * channelWeightNormalized ' )) / (size(obj .eyeChannelWeightNormalized ,2 ));
91+ isEye = similarity > similarityThreshold ;
92+ end ;
93+
94+ function [isEye similarity ] = detectFromScalpmapObj(obj , scalpmapObj , similarityThreshold , varargin )
95+ % [isEye similarity] = detectFromScalpmapObj(obj, scalpmapObj, similarityThreshold)
96+
97+ if nargin < 3
98+ similarityThreshold = obj .similarityThreshold ;
99+ end ;
100+
101+ [isEye similarity ] = detectFromInterpolatedChannelWeight(obj , scalpmapObj .originalChannelWeight(: ,: ), similarityThreshold );
102+ end ;
103+
104+ function [isEye similarity scalpmapObj ] = detectFromEEG(obj , EEG , similarityThreshold , varargin )
105+ % [isEye similarity scalpmapObj] = detectFromEEG(obj, EEG, similarityThreshold)
106+ % EEG data MUST have been highpassed at 1 Hz.
107+
108+ if nargin < 3
109+ similarityThreshold = obj .similarityThreshold ;
110+ end ;
111+
112+ if isempty(EEG .icachansind )
113+ EEG.icachansind= 1 : length(EEG .chanlocs );
114+ end ;
115+
116+ [isEye , similarity , scalpmapObj ] = detectFromChannelWeight(obj , EEG .icawinv , EEG .chanlocs(EEG .icachansind ), similarityThreshold );
117+ isEye = findEyeMovementICsByPowerRatio(obj , EEG , isEye , similarity );
118+ end ;
119+
120+ function [isEye similarity scalpmapObj ] = detectFromChannelWeight(obj , channelWeight , channelLocation , similarityThreshold , varargin )
121+ % [isEye similarity] = detectFromChannelWeight(obj, channelWeight, channelLocation, similarityThreshold)
122+ if nargin < 4
123+ similarityThreshold = obj .similarityThreshold ;
124+ end ;
125+
126+ scalpmapObj = scalpmap ;
127+ scalpmapObj = scalpmapObj .addFromChannels(channelWeight , channelLocation );
128+ [isEye similarity ] = detectFromScalpmapObj(obj , scalpmapObj , similarityThreshold );
129+ end ;
130+
131+ function [isEye similarity scalpmapObj ] = detectFromStudy(obj , STUDY , ALLEEG , similarityThreshold , varargin )
132+ scalpmapObj = scalpmapOfStudy(STUDY , ALLEEG , [], ' normalizePolarity' false );
133+ [isEye similarity ] = detectFromScalpmapObj(obj , scalpmapObj , similarityThreshold );
134+ end ;
135+
136+ function isEye = findEyeMovementICsByPowerRatio(obj , EEG , isEye , eyeCatchSimilarity )
137+ % isEye = findEyeMovementICsByPowerRatio(EEG, isEye, eyeCatchSimilarity)
138+ % EEG data MUST have been highpassed at 1 Hz.
139+
140+ criticalFreq = 3 ;
141+ powerRatioThreshold = 100 ;
142+ timeRatioOfPowerRatioTooHigh = 0.01 ;
143+
144+ icNumbers = find(eyeCatchSimilarity > 0.85 & ~isEye );
145+
146+ wname = ' cmor1-1.5'
147+ T = 1 / EEG .srate ;
148+ frequencyRange = [1 15 ];
149+ numberOfFrequencies = 20 ;
150+ [scales , freqs ] = freq2scales(frequencyRange(1 ), frequencyRange(2 ), numberOfFrequencies , wname , T );
151+
152+ if isempty(EEG .icaact )
153+ if isempty(EEG .icachansind )
154+ EEG.icachansind = 1 : size(EEG .data ,1 );
155+ end ;
156+ EEG.icaact = EEG .icaweights * EEG .icasphere * EEG .data(EEG .icachansind ,: );
157+ end ;
158+
159+ powerRatioTooHigh = false(length(icNumbers ), size(EEG .data ,2 ));
160+
161+ for i= 1 : length(icNumbers )
162+
163+ tfdecomposition = cwt(EEG .icaact(icNumbers(i ),: )' ,scales , wname );
164+
165+ tfdecomposition = abs(tfdecomposition );
166+ powerRatio = sum(tfdecomposition(freqs < criticalFreq ,: ).^2 ) ./ sum(tfdecomposition(freqs >= criticalFreq ,: ).^2 );
167+ powerRatioTooHigh(i ,: ) = powerRatio > powerRatioThreshold ;
168+ end ;
169+
170+ eyeMovementICs = mean(powerRatioTooHigh ,2 ) > timeRatioOfPowerRatioTooHigh ; % more than 1% of time have a very high low-frequency activity
171+ if any(eyeMovementICs )
172+ fprintf(' Found %d additional eye movement ICs: %s\n ' eyeMovementICs ), ...
173+ strjoin_adjoiner_first(' , ' num2str , find(eyeMovementICs ), ' UniformOutput' false )));
174+ end ;
175+ isEye(icNumbers(eyeMovementICs )) = true ;
176+ end ;
177+ end
178+ end
0 commit comments