[422372]: / functions / miscfunc / loc_subsets.m

Download this file

227 lines (196 with data), 9.1 kB

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
% LOC_SUBSETS - Separate channels into maximally evenly-spaced subsets.
% This is achieved by exchanging channels between subsets so as to
% increase the sum of average of distances within each channel subset.
% Usage:
% >> subset = loc_subsets(chanlocs, nchans); % select an evenly spaced nchans
% >> [subsets subidx pos] = loc_subsets(chanlocs, nchans, plotobj, plotchans, keepchans);
%
% Inputs:
%
% chanlocs - EEGLAB dataset channel locations structure (e.g., EEG.chanlocs)
%
% nchans - (1,N) matrix containing number of channels that should be in each
% of N channel subsets respectively. if total number of channels in
% all subsets is less than the number of EEG channels in the chanlocs
% structure, N+1 subsets are created. The last subset includes the
% remaining channels.
%
% Optional inputs:
%
% plotobj - (true|false|1|0) plot the time course of the objective function
% {default: false|0)
% plotchans - (true|false|1|0) make 3-D plots showing the channel locations of
% the subsets {default: false|0)
% keepchans - (cell array) channels that has to be kept in each set and
% not undergo optimization. You can use this option to make
% sure certain channels will be assigned to each set. For
% example, to keep channels 1:10 to the first subset and
% channels 20:30 to the second, use keepchans = {[1:10], [20:30]}.
% To only keeps channels 1:5 in the first set, use
% keepchans = {1:5}. {default: {}}
%
% Outputs:
%
% subsets - {1,N} or {1,N+1} cell array containing channel indices for each
% requested channel subset.
% subidx - (1, EEG.nbchans) vector giving the index of the subset associated
% with each channel.
% pos - (3,N) matrix, columns containing the cartesian (x,y,z) channel
% positions plotted.
% Example:
%
% % Create three sub-montages of a 256-channel montage containing 32, 60, and 100
% % channels respectively.The 64 remaining channels will be put into a fourth subset.
% % Also visualize time course of optimization and the channel subset head locations.
%
% >> subset = loc_subsets(EEG.chanlocs, [32 60 100], true, true);
%
% Author: Nima Bigdely Shamlo, SCCN/INC/UCSD, 2007
% Copyright (C) 2007 Nima Bigdely Shamlo, SCCN/INC/UCSD, nima@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.
% if plotOptimization|plotSubsets in line 130 removed by nima 3/6/2007
% line 133, removed num2str removed by nima 3/6/2007
function [subset idx pos] = loc_subsets(chanlocs, numberOfChannelsInSubset, plotOptimization, plotSubsets, mandatoryChannelsForSet);
if sum(numberOfChannelsInSubset)> length(chanlocs)
error('Total channels in requested subsets larger than number of EEG channels.');
end
if min(numberOfChannelsInSubset) < 2
error('Number of channels in the requested subsets must be >= 2.');
end
rand('state',0);
if nargin < 5
mandatoryChannelsForSet = {};
end
if nargin < 4
plotSubsets = false;
end
if nargin < 3
plotOptimization = false;
end
pos=[cell2mat({chanlocs.X}); cell2mat({chanlocs.Y}); cell2mat({chanlocs.Z});];
dist = squareform(pdist(pos'));
nChans = length(chanlocs);
idx = ones(nChans,1);
setId = cell(1, length(numberOfChannelsInSubset)); % cell array containing channels in each set
remainingChannels = 1:nChans; % channels that are to be assigned to subsets
% channels that have to stay in their original subset (as in
% mandatoryChannelsForSet) and should not be re-assigned
allMandatoryChannels = cell2mat(mandatoryChannelsForSet);
% assign requested mandatory channels to subset
for i=1:length(mandatoryChannelsForSet)
setId{i} = mandatoryChannelsForSet{i};
remainingChannels(mandatoryChannelsForSet{i}) = NaN; % flag with Nan so they can be deleted later, this is to keep indexing simple
end
remainingChannels(isnan(remainingChannels)) = [];
r = remainingChannels(randperm(length(remainingChannels)));
% randomly assign remaining channels to subsets
for i=1:length(numberOfChannelsInSubset)
numberOfChannelsTobeAddedToSubset = numberOfChannelsInSubset(i) - length(setId{i})
setId{i} = [setId{i} r(1:numberOfChannelsTobeAddedToSubset)];
r(1:numberOfChannelsTobeAddedToSubset) = [];
end
if length(r) > 0
setId{length(numberOfChannelsInSubset) + 1} = r; % last set gets remaining channels
end
if plotOptimization|plotSubsets
fprintf(['Creating total of ' num2str(length(setId)) ' channel subsets:\n']);
end
if plotOptimization
figure;
xp = floor(sqrt(length(setId)));
yp = ceil(length(setId)/xp);
end
counter = 1;
exchangeHappened = true;
while exchangeHappened
exchangeHappened = false;
for set1 = 1:length(setId)
for set2 = (set1 + 1):length(setId)
for i = 1:length(setId{set1})
for j = 1:length(setId{set2})
chan(1) = setId{set1}(i);
chan(2) = setId{set2}(j);
if cost_of_exchanging_channels(chan,[set1 set2], setId, dist) < 0 && ~any(ismember(chan, allMandatoryChannels))
setId{set1}(find(setId{set1} == chan(1))) = chan(2);
setId{set2}(find(setId{set2} == chan(2))) = chan(1);
sumDistances(counter) = 0;
for s = 1:length(setId)
sumDistances(counter) = sumDistances(counter) ...
+ (sum(sum(dist(setId{s},setId{s}))) / length(setId{s}));
end
if plotOptimization
plot(1:counter,sumDistances,'-b');
xlabel('number of exchanges');
ylabel('sum mean distances within each channel subset');
drawnow;
else
if mod(counter, 20) ==0
fprintf('number of exchanges = %d\nsum of mean distances = %g\n',...
counter, sumDistances(counter));
end
end
counter = counter + 1;
exchangeHappened = true;
end
end
end
end
end
end
for set = 1:length(setId)
idx(setId{set}) = set;
% legendTitle{set} = ['subset ' num2str(set)];
end
subset = setId;
if plotSubsets
FIG_OFFSET = 40;
fpos = get(gcf,'position');
figure('position',[fpos(1)+FIG_OFFSET,fpos(2)-FIG_OFFSET,fpos(3),fpos(4)]);
scatter3(pos(1,:), pos(2,:), pos(3,:),100,idx,'fill');
axis equal;
%legend(legendTitle); it does not work propr
th=title('Channel Subsets');
%set(th,'fontsize',14)
end
if length(r)>0 && (plotOptimization|plotSubsets)
fprintf('The last subset returned contains the %d unused channels.\n',...
length(r));
end
function cost = cost_of_exchanging_channels(chan, betweenSets, setId, dist);
mirrorChan(1) = 2;
mirrorChan(2) = 1;
for i=1:2
newSetId{betweenSets(i)} = setId{betweenSets(i)};
newSetId{betweenSets(i)}(find(newSetId{betweenSets(i)} == chan(i))) = chan(mirrorChan(i));
end
cost = 0;
for i=betweenSets
distSumBefore{i} = sum(sum(dist(setId{i},setId{i})));
distSumAfter{i} = sum(sum(dist(newSetId{i},newSetId{i})));
cost = cost + (distSumBefore{i} - distSumAfter{i}) / length(setId{i});
end
%cost = (distSumAfter{1} > distSumBefore{1}) && (distSumAfter{2} > distSumBefore{2});