Attachment 'process_beamformer_mcb.m'
Download 1 function varargout = process_beamformer_mcb( varargin )
2 % PROCESS_BEAMFORMER_TEST:
3
4 % @=============================================================================
5 % This software is part of the Brainstorm software:
6 % http://neuroimage.usc.edu/brainstorm
7 %
8 % Copyright (c)2000-2015 University of Southern California & McGill University
9 % This software is distributed under the terms of the GNU General Public License
10 % as published by the Free Software Foundation. Further details on the GPL
11 % license can be found at http://www.gnu.org/copyleft/gpl.html.
12 %
13 % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
14 % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
15 % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
16 % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
17 % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
18 %
19 % For more information type "brainstorm license" at command prompt.
20 % =============================================================================@
21 %
22 % Authors:
23
24 macro_methodcall;
25 end
26
27
28 %% ===== GET DESCRIPTION =====
29 function sProcess = GetDescription() %#ok<DEFNU>
30 % Description the process
31 sProcess.Comment = 'Beamformer: Maximum Contrast [Experimental]';
32 sProcess.FileTag = '';
33 sProcess.Category = 'Custom';
34 sProcess.SubGroup = 'Sources';
35 sProcess.Index = 330;
36 % Definition of the input accepted by this process
37 sProcess.InputTypes = {'data'};
38 sProcess.OutputTypes = {'results'};
39 sProcess.nInputs = 1;
40 sProcess.nMinFiles = 1;
41 sProcess.isSeparator = 1;
42 % Definition of the options
43 sProcess.options.result_comm.Comment = 'Comment: ';
44 sProcess.options.result_comm.Type = 'text';
45 sProcess.options.result_comm.Value = '';
46 % === CORTICAL CONSTRAINED ORIENTATION
47 sProcess.options.label_orient.Comment = '<HTML><BR> Sources orientation:';
48 sProcess.options.label_orient.Type = 'label';
49 sProcess.options.oriconstraint.Comment = {'Unconstrained', 'Constrained (normal to cortex)'};
50 sProcess.options.oriconstraint.Type = 'radio';
51 sProcess.options.oriconstraint.Value = 1;
52 % === F-MAP TIME RANGE
53 sProcess.options.label_range.Comment = '<HTML><BR> ';
54 sProcess.options.label_range.Type = 'label';
55 sProcess.options.fmaps_range.Comment = 'Time range of active state: ';
56 sProcess.options.fmaps_range.Type = 'poststim';
57 sProcess.options.fmaps_range.Value = [];
58 % === BASELINE TIME RANGE FOR Z STATISTICS CALCULATION
59 sProcess.options.baseline_time.Comment = 'Baseline for F-statistics: ';
60 sProcess.options.baseline_time.Type = 'baseline';
61 sProcess.options.baseline_time.Value = [];
62 % === F-MAP TIME WINDOW SIZE
63 sProcess.options.fmap_size.Comment = 'F-statistic window size: ';
64 sProcess.options.fmap_size.Type = 'value';
65 sProcess.options.fmap_size.Value = {0.2958, 'ms', 1};
66 % === F-MAP TEMPORAL RESOLUTION
67 sProcess.options.fmaps_tresolution.Comment = 'Temporal resolution: ';
68 sProcess.options.fmaps_tresolution.Type = 'value';
69 sProcess.options.fmaps_tresolution.Value = {0.001, 'ms', 1};
70 % === REGULARIZATION
71 sProcess.options.reg.Comment = 'Regularization parameter: ';
72 sProcess.options.reg.Type = 'value';
73 sProcess.options.reg.Value = {0.03, '%', 4};
74 % === SENSOR TYPES
75 sProcess.options.label_sensor.Comment = '<HTML><BR> ';
76 sProcess.options.label_sensor.Type = 'label';
77 sProcess.options.sensortypes.Comment = 'Sensor types or names (empty=all): ';
78 sProcess.options.sensortypes.Type = 'text';
79 sProcess.options.sensortypes.Value = 'MEG';
80 % === Output types
81 sProcess.options.outputValue.Comment = {'Power / Variance', 'Variance / Variance', 'Variance (user defined baseline) / Variance'};
82 sProcess.options.outputValue.Type = 'radio';
83 sProcess.options.outputValue.Value = 2;
84
85 end
86
87
88 %% ===== FORMAT COMMENT =====
89 function Comment = FormatComment(sProcess) %#ok<DEFNU>
90 Comment = sProcess.Comment;
91 end
92
93
94 %% ===== RUN =====
95 function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
96 % Initialize returned list of files
97 OutputFiles = {};
98
99 % Get option values
100 BaselineTime = sProcess.options.baseline_time.Value{1};
101 FmapRange = sProcess.options.fmaps_range.Value{1};
102 Reg = sProcess.options.reg.Value{1};
103 SensorTypes = sProcess.options.sensortypes.Value;
104
105 FmapSize = sProcess.options.fmap_size.Value{1};
106 FmapTResolu = sProcess.options.fmaps_tresolution.Value{1};
107
108 result_comment = sProcess.options.result_comm.Value;
109 if ~isempty(result_comment)
110 result_comment = [result_comment ':'];
111 end
112 % ===== LOAD THE DATA =====
113 % Read the first file in the list, to initialize the loop
114 DataMat = in_bst(sInputs(1).FileName, [], 0);
115 nChannels = size(DataMat.F,1);
116 nTime = size(DataMat.F,2);
117 Time = DataMat.Time;
118
119 % ===== PROCESS THE TIME WINDOWS =====
120 if FmapRange(1) > FmapRange(2)
121 bst_report('Error', sProcess, sInputs, 'The setting of time range of active state is incorrect.');
122 end
123
124 if FmapRange(1) < Time(1)
125 bst_report('Warning', sProcess, sInputs, 'The start for time range of active state is reset to the first time point of data');
126 FmapRange(1) = Time(1);
127 end
128
129 if FmapRange(2) > Time(end)
130 bst_report('Warning', sProcess, sInputs, 'The end for time range of active state is reset to the end point of data');
131 FmapRange(2) = Time(end);
132 end
133
134 if (FmapRange(1)+FmapSize) > FmapRange(2)
135 bst_report('Warning', sProcess, sInputs, 'The F statistic window size is too large and reset to the same as the time range of active state.');
136 FmapSize = FmapRange(2) - FmapRange(1);
137 end
138
139 ActiveTime = FmapRange;
140
141 if FmapTResolu == 0
142 nFmaps = 1;
143 FmapTResolu = 1;
144 FmapSize = FmapRange(2) - FmapRange(1);
145 if FmapRange(1)+FmapSize ~= FmapRange(2)
146 bst_report('Warning', sProcess, sInputs, 'Temporal resolution should not be 0 ms. The F statistic window size is reset to the same as the time range of interest.');
147 end
148 end
149
150 FmapPoints = panel_time('GetTimeIndices', Time, [FmapRange(1)+FmapSize FmapRange(2)]);
151 if length(FmapPoints)<=1
152 nFmaps = 1;
153 else
154 nFmaps = length((FmapRange(1)+FmapSize):FmapTResolu:FmapRange(2));
155 end
156
157
158 HalfFmapSize= FmapSize/2;
159 FmapTimeList= zeros(nFmaps,2);
160 for i=1:nFmaps
161 FmapTimeList(i,:) = FmapRange(1) + FmapTResolu*(i-1) + [0 FmapSize] ;
162 end
163
164
165
166 % ===== LOAD CHANNEL FILE =====
167 % Load channel file
168 ChannelMat = in_bst_channel(sInputs(1).ChannelFile);
169 % Find the MEG channels
170 % iMEG = good_channel(ChannelMat.Channel, [], 'MEG');
171 % iEEG = good_channel(ChannelMat.Channel, [], 'EEG');
172 % iSEEG = good_channel(ChannelMat.Channel, [], 'SEEG');
173 % iECOG = good_channel(ChannelMat.Channel, [], 'ECOG');
174 iChannels = channel_find(ChannelMat.Channel, SensorTypes);
175
176 % ===== LOAD HEAD MODEL =====
177 % Get channel study
178 [sChannelStudy, iChannelStudy] = bst_get('ChannelFile', sInputs(1).ChannelFile);
179 % Load the default head model
180 HeadModelFile = sChannelStudy.HeadModel(sChannelStudy.iHeadModel).FileName;
181 sHeadModel = load(file_fullpath(HeadModelFile));
182 % Get number of sources
183 nSources = length(sHeadModel.GridLoc);
184
185 if (sProcess.options.oriconstraint.Value == 2) && (isequal(sHeadModel.HeadModelType,'volume') || isempty(sHeadModel.GridOrient))
186 bst_report('Error', sProcess, sInputs, 'No dipole orientation for cortical constrained beamformer estimation.');
187 % Stop the process
188 return;
189 end
190
191
192 % ===== LOAD THE DATA =====
193 % Find the indices for covariance calculation
194 iActiveTime = panel_time('GetTimeIndices', Time, ActiveTime);
195 iBaselineTime = panel_time('GetTimeIndices', Time, BaselineTime);
196 % Initialize the covariance matrices
197 ActiveCov = zeros(nChannels, nChannels);
198 nTotalActive = zeros(nChannels, nChannels);
199
200
201 nFmapPoints = length(panel_time('GetTimeIndices', Time, FmapTimeList(1,:)));
202 iFmapTime = zeros(nFmaps, nFmapPoints);
203 for i = 1:nFmaps
204 if FmapTimeList(i,1) < Time(1) || FmapTimeList(i,2) > Time(end)
205 % Add an error message to the report
206 bst_report('Error', sProcess, sInputs, 'One fmap time window is not within the data time range.');
207 % Stop the process
208 return;
209 end
210 single_iFmap = panel_time('GetTimeIndices', Time, FmapTimeList(i,:));
211 iFmapTime(i,:) = single_iFmap(1:nFmapPoints);
212 end
213 % Initialize the covariance matrices
214 FmapActiveCov = zeros(nChannels, nChannels, nFmaps);
215 nTotalFmapActive = zeros(nChannels, nChannels, nFmaps);
216
217 bst_progress('start', 'Applying process: MCB', 'Calulating covariance matrix...', 0, length(sInputs)*nFmaps*4+2);
218 % Reading all the input files in a big matrix
219 for i = 1:length(sInputs)
220 % Read the file #i
221 DataMat = in_bst(sInputs(i).FileName, [], 0);
222 % Check the dimensions of the recordings matrix in this file
223 if (size(DataMat.F,1) ~= nChannels) || (size(DataMat.F,2) ~= nTime)
224 % Add an error message to the report
225 bst_report('Error', sProcess, sInputs, 'One file has a different number of channels or a different number of time samples.');
226 % Stop the process
227 return;
228 end
229
230 % Get good channels
231 iGoodChan = find(DataMat.ChannelFlag == 1);
232
233
234 % Average baseline values
235 FavgActive = mean(DataMat.F(iGoodChan,iActiveTime), 2);
236 % Remove average
237 DataActive = bst_bsxfun(@minus, DataMat.F(iGoodChan,iActiveTime), FavgActive);
238 % Compute covariance for this file
239 fileActiveCov = DataMat.nAvg .* (DataActive * DataActive');
240 % Add file covariance to accumulator
241 ActiveCov(iGoodChan,iGoodChan) = ActiveCov(iGoodChan,iGoodChan) + fileActiveCov;
242 nTotalActive(iGoodChan,iGoodChan) = nTotalActive(iGoodChan,iGoodChan) + length(iActiveTime);
243
244
245 for j = 1:nFmaps
246 % Average baseline values
247 if sProcess.options.outputValue.Value == 1
248 FavgFmapActive = 0;
249 elseif sProcess.options.outputValue.Value == 2
250 FavgFmapActive = mean(DataMat.F(iGoodChan,iFmapTime(j,:)), 2);
251 elseif sProcess.options.outputValue.Value == 3
252 FavgFmapActive = mean(DataMat.F(iGoodChan,iBaselineTime), 2);
253 end
254
255 % Remove average
256 DataFmapActive = bst_bsxfun(@minus, DataMat.F(iGoodChan,iFmapTime(j,:)), FavgFmapActive);
257 bst_progress('inc',1);
258 % Compute covariance for this file
259 fileFmapActiveCov = DataMat.nAvg .* (DataFmapActive * DataFmapActive');
260 bst_progress('inc',1);
261 % Add file covariance to accumulator
262 FmapActiveCov(iGoodChan,iGoodChan,j) = FmapActiveCov(iGoodChan,iGoodChan,j) + fileFmapActiveCov;
263 bst_progress('inc',1);
264 nTotalFmapActive(iGoodChan,iGoodChan,j) = nTotalFmapActive(iGoodChan,iGoodChan,j) + nFmapPoints;
265 bst_progress('inc',1);
266 end
267
268 end
269 % Remove zeros from N matrix
270 if sProcess.options.oriconstraint.Value == 1,
271 nTotalActive(nTotalActive <= 1) = 2;
272 end
273 nTotalFmapActive(nTotalFmapActive <= 1) = 2;
274 bst_progress('inc',1);
275 % Divide final matrix by number of samples
276 if sProcess.options.oriconstraint.Value == 1,
277 ActiveCov = ActiveCov ./ (nTotalActive - 1);
278 end
279 FmapActiveCov = FmapActiveCov ./ (nTotalFmapActive - 1);
280 bst_progress('inc',1);
281 bst_progress('stop');
282 % ===== PROCESS =====
283 % Number of channels used to compute sources
284 nUsedChannels = length(iChannels);
285
286 % Extract the covariance matrix of the used channels
287 MinVarCov = ActiveCov;
288 ActiveCov = ActiveCov(iChannels,iChannels);
289 NoiseCovMat = load(file_fullpath(sChannelStudy.NoiseCov(1).FileName));
290 NoiseCov = NoiseCovMat.NoiseCov(iChannels,iChannels);
291 MinVarCov = MinVarCov(iChannels,iChannels);
292 FmapActiveCov = FmapActiveCov(iChannels,iChannels,:);
293
294 %% Calculate the inverse of (C+alpha*I)
295 % eigValues = eig(MinVarCov);
296 % Reg_alpha = Reg / 100 * max(eigValues);
297 % invMinVarCovI = inv(MinVarCov+Reg_alpha*eye(nUsedChannels));
298
299 [U,S,V] = svd(MinVarCov);
300 S = diag(S); % Covariance = Cm = U*S*U'
301 Si = diag(1 ./ (S + S(1) * Reg / 100)); % 1/(S + lambda I)
302 invMinVarCovI = V*Si*U'; % V * 1/(S + lambda I) * U' = Cm^(-1)
303
304 % Initilize the ImagingKernel (spatial filter) and ImageGridAmp (f value)
305 ImagingKernel = zeros(nSources, nUsedChannels);
306 ImageGridAmp = zeros(nSources, nFmaps);
307
308 % Set the dipole positions for the computation of sources
309 GridLoc = sHeadModel.GridLoc;
310
311
312 % Get forward field
313 if sProcess.options.oriconstraint.Value == 2; % constrained LCMV
314 Kernel = bst_gain_orient(sHeadModel.Gain(iChannels,:), sHeadModel.GridOrient);
315 else
316 Kernel = sHeadModel.Gain(iChannels,:);
317 end
318 %Kernel(abs(Kernel(:)) < eps) = eps; % Set zero elements to strictly non-zero
319 bst_progress('start', 'Applying process: MCB', 'Calculating source significance...', 0, nSources*(nFmaps+1));
320
321 %% Compute the spatial filter and f value for each position
322 if sProcess.options.oriconstraint.Value == 1;
323 % Initial the index for gain matrix
324 iGain = [1 2 3];
325
326 for i = 1:nSources
327 % Compute A = inv(C+alpha*I)*Lr
328 A = invMinVarCovI * Kernel(:,iGain);
329
330 % Compute B = Lr'*A
331 B = Kernel(:,iGain)' * A;
332
333
334 % == Compute orientation using maxium constrast criterion ==
335 % Compute P = A'*Ca*A
336 P = A' * ActiveCov * A;
337 % Compute Q = A'*Cc*A
338 Q = A' * NoiseCov * A;
339
340 % Regularize the matrix Q to avoid singular problem
341 [U,S,V] = svd(Q);
342 S = diag(S); % Covariance = Cm = V*S*U'
343 Si = diag(1 ./ (S + S(1) * 0.0000000001)); % 1/(S + lambda I)
344 invQ = V*Si*U'; % V * 1/(S + lambda I) * U' = Cm^(-1)
345
346 % Compute the dipole orientation
347 % (the eigenvector corresponding to maximum eigenvalue of inv(Q)*P)
348 [eigVectors,eigValues] = eig(invQ*P);
349 % check whether eigValues are saved as matrix or vector
350 if(min(size(eigValues))==1)
351 [tmp, imax] = max(eigValues);
352 else
353 [tmp, imax] = max(diag(eigValues));
354 end
355 DipoleOri = eigVectors(:,imax);
356
357 % Compute the spatial filter
358 SpatialFilter = (A * DipoleOri) / (DipoleOri' * B * DipoleOri);
359
360 varControl = SpatialFilter'* NoiseCov * SpatialFilter;
361 bst_progress('inc',1);
362 for j = 1:nFmaps
363 % Compute the contrast of source power during active state and control state
364 varActive = SpatialFilter'*FmapActiveCov(:,:,j)*SpatialFilter;
365
366 ImageGridAmp(i,j) = varActive / varControl;
367
368 bst_progress('inc',1);
369 end
370 % Save teh result
371 ImagingKernel(i,:) = SpatialFilter / sqrt(varControl);
372
373
374 % The index of gain for the next dipole positions
375 iGain = iGain + 3;
376 end
377
378 else
379 for i = 1:nSources
380 % Compute the spatial filter with cortical constrained dipole orientation
381
382 % Compute A = inv(C+alpha*I)*Lr
383 A = invMinVarCovI * Kernel(:,i);
384
385 % Compute B = Lr'*A
386 B = Kernel(:,i)'*A;
387
388 % Compute the spatial filter
389 SpatialFilter = A / B;
390
391
392 % compute the variance of control state
393 varControl = SpatialFilter'* NoiseCov * SpatialFilter;
394 bst_progress('inc',1);
395 for j = 1:nFmaps
396 % Compute the contrast of source power during active state and control state
397 varActive = SpatialFilter'*FmapActiveCov(:,:,j)*SpatialFilter;
398
399 ImageGridAmp(i,j) = varActive / varControl;
400 %ImageGridAmp(i,j) = Fvalue;
401 bst_progress('inc',1);
402 end
403 % Save teh result
404 ImagingKernel(i,:) = SpatialFilter / sqrt(varControl);
405
406 end
407 end
408
409 bst_progress('stop');
410
411 bst_progress('start', 'Applying process: MCB', 'Interpolating results...', 0, 1);
412 %bst_progress('text', ['Applying process: ' sProcess.Comment ' [Interpolating results]']);
413 if nFmaps == 1
414 FmapRangePoints = panel_time('GetTimeIndices', Time, FmapRange);
415 ImageGridAmpOriginal = ImageGridAmp;
416 ImageGridAmp = zeros(nSources,nTime);
417 for i = FmapRangePoints(1):FmapRangePoints(end)
418 ImageGridAmp(:,i) = ImageGridAmpOriginal;
419 end
420 TimeIndex = Time;
421 else % Interpolate the f maps to have the same temporal resolution as the data
422
423 ImageGridAmpOriginal = ImageGridAmp;
424 ImageGridAmp = zeros(nSources,nTime);
425 for i=1:(nFmaps-1)
426 InterpolateTimeWindow = FmapRange(1)+HalfFmapSize+[(i-1)*FmapTResolu i*FmapTResolu];
427 iInterpolateTime = panel_time('GetTimeIndices', Time, InterpolateTimeWindow);
428 nInterpolateTime = length(iInterpolateTime);
429
430 ImageGridAmp(:,iInterpolateTime(1)) = ImageGridAmpOriginal(:,i);
431 ImageGridAmp(:,iInterpolateTime(end)) = ImageGridAmpOriginal(:,i+1);
432
433 for j=2:(nInterpolateTime-1)
434 InterpolatePercentage = (j-1)/(nInterpolateTime-1);
435 ImageGridAmp(:,iInterpolateTime(j)) = ImageGridAmpOriginal(:,i+1)*InterpolatePercentage + ImageGridAmpOriginal(:,i)*(1-InterpolatePercentage);
436 end
437 end
438 %%%%%%%%%%%%
439
440 % InterpolateTimeWindow = FmapRange(1)+ [0 HalfFmapSize];
441 % iInterpolateTime = panel_time('GetTimeIndices', Time, InterpolateTimeWindow);
442 % nInterpolateTime = length(iInterpolateTime);
443 %
444 % for j=1:(nInterpolateTime-1)
445 % ImageGridAmp(:,iInterpolateTime(j)) = ImageGridAmpOriginal(:,1);
446 % end
447 %
448 % InterpolateTimeWindow = FmapRange(1)+HalfFmapSize+(nFmaps-1)*FmapTResolu+[0 HalfFmapSize];
449 % iInterpolateTime = panel_time('GetTimeIndices', Time, InterpolateTimeWindow);
450 % nInterpolateTime = length(iInterpolateTime);
451 %
452 % for j=2:nInterpolateTime
453 % ImageGridAmp(:,iInterpolateTime(j)) = ImageGridAmpOriginal(:,end);
454 % end
455 %
456
457 TimeIndex = Time;
458 end
459 bst_progress('inc',1);
460 bst_progress('stop');
461
462 %bst_progress('text', ['Applying process: ' sProcess.Comment ' [Saving results]']);
463 % ===== SAVE THE RESULTS =====
464 % Create a new data file structure
465 ResultsMat = db_template('resultsmat');
466 ResultsMat.ImagingKernel = [];
467 ResultsMat.ImageGridAmp = ImageGridAmp;
468 ResultsMat.nComponents = 1; % 1 or 3
469
470 timestring = sprintf('%.01f_%.01f%s',(FmapRange(1)+HalfFmapSize)*1000,(FmapRange(1)+HalfFmapSize+nFmaps*FmapTResolu)*1000,'ms');
471
472 if nFmaps == 1
473 if sProcess.options.outputValue.Value == 3
474 timestring2 = sprintf('(bs:%.1f_%.1f%s)',BaselineTime(1)*1000,BaselineTime(2)*1000,'ms');
475 else
476 timestring2 = [];
477 end
478 else
479 if sProcess.options.outputValue.Value == 3
480 timestring2 = sprintf('(ws:%.1f%s,tr:%.1f%s,bs:%.1f_%.1f%s)',FmapSize*1000,'ms',FmapTResolu*1000,'ms',BaselineTime(1)*1000,BaselineTime(2)*1000,'ms');
481 else
482 timestring2 = sprintf('(ws:%.1f%s,tr:%.1f%s)',FmapSize*1000,'ms',FmapTResolu*1000,'ms');
483 end
484 end
485
486 if sProcess.options.oriconstraint.Value == 1;
487 ostr = 'Unconstr';
488 else
489 ostr = 'Constr';
490 end
491
492
493 if sProcess.options.outputValue.Value == 1
494 ResultsMat.Comment = [result_comment 'MCB: F-statistics pow/var(' ostr ')|' timestring timestring2] ;
495 elseif sProcess.options.outputValue.Value == 2
496 ResultsMat.Comment = [result_comment 'MCB: F-statistics var/var(' ostr ')|' timestring timestring2];
497 elseif sProcess.options.outputValue.Value == 3
498 ResultsMat.Comment = [result_comment 'MCB: F-statistics var(b)/var(' ostr ')|' timestring timestring2];
499 end
500
501
502 ResultsMat.Function = 'MaximumContrastBeamformerResult';
503 ResultsMat.Time = TimeIndex; % Leave it empty if using ImagingKernel
504 ResultsMat.DataFile = [];
505 ResultsMat.HeadModelFile = HeadModelFile;
506 ResultsMat.HeadModelType = sHeadModel.HeadModelType;
507 ResultsMat.ChannelFlag = [];
508 ResultsMat.GoodChannel = iChannels;
509 ResultsMat.SurfaceFile = sHeadModel.SurfaceFile;
510 ResultsMat.GridLoc = GridLoc;
511
512 % === NOT SHARED ===
513 % Get the output study (pick the one from the first file)
514 iStudy = sInputs(1).iStudy;
515 % Create a default output filename
516 OutputFiles{1} = bst_process('GetNewFilename', fileparts(sInputs(1).FileName), 'results_MCB_amp');
517
518 % Save on disk
519 save(OutputFiles{1}, '-struct', 'ResultsMat');
520 % Register in database
521 db_add_data(iStudy, OutputFiles{1}, ResultsMat);
522
523
524 % == Save the spatial filter as ImagingKernel ==
525 % Create a new data file structure
526 ResultsMat2 = db_template('resultsmat');
527 ResultsMat2.ImagingKernel = ImagingKernel;
528 ResultsMat2.ImageGridAmp = [];
529 ResultsMat2.nComponents = 1; % 1 or 3
530 timestring = sprintf('%.01f_%.01f%s',ActiveTime(1)*1000,ActiveTime(2)*1000,'ms');
531 if sProcess.options.oriconstraint.Value == 1;
532 ResultsMat2.Comment = [result_comment 'MCB: SpatilFilter(Unconstr)|' timestring];
533 else
534 ResultsMat2.Comment = [result_comment 'MCB: SpatilFilter(Constr)|' timestring];
535 end
536 ResultsMat2.Function = 'MaximumContrastBeamformerFilter';
537 ResultsMat2.Time = []; % Leave it empty if using ImagingKernel
538 ResultsMat2.DataFile = [];
539 ResultsMat2.HeadModelFile = HeadModelFile;
540 ResultsMat2.HeadModelType = sHeadModel.HeadModelType;
541 ResultsMat2.ChannelFlag = [];
542 ResultsMat2.GoodChannel = iChannels;
543 ResultsMat2.SurfaceFile = sHeadModel.SurfaceFile;
544 ResultsMat2.GridLoc = GridLoc;
545
546 % === SHARED ==
547 % Get the output study (pick the one from the first file)
548 iStudy = iChannelStudy;
549 % Create a default output filename
550 OutputFiles{2} = bst_process('GetNewFilename', fileparts(sInputs(1).ChannelFile), 'results_MCB_KERNEL');
551
552 % Save on disk
553 save(OutputFiles{2}, '-struct', 'ResultsMat2');
554 % Register in database
555 db_add_data(iStudy, OutputFiles{2}, ResultsMat2);
556 %%===========
557 %
558 % if (sProcess.options.oriconstraint.Value == 1) && (isempty(sHeadModel.GridOrient)==0)
559 % % ===== SAVE THE RESULTS =====
560 % % Create a new data file structure
561 % ResultsMat3 = db_template('resultsmat');
562 % ResultsMat3.ImagingKernel = [];
563 % ResultsMat3.ImageGridAmp = [OriDifference OriDifference];
564 % ResultsMat3.nComponents = 1; % 1 or 3
565 % ResultsMat3.Comment = 'MCB: Orientation Difference(Unconstr)';
566 % ResultsMat3.Function = 'MaximumContrastBeamformerOriDiff';
567 % ResultsMat3.Time = [1 1]; % Leave it empty if using ImagingKernel
568 % ResultsMat3.DataFile = [];
569 % ResultsMat3.HeadModelFile = HeadModelFile;
570 % ResultsMat3.HeadModelType = sHeadModel.HeadModelType;
571 % ResultsMat3.ChannelFlag = [];
572 % ResultsMat3.GoodChannel = iChannels;
573 % ResultsMat3.SurfaceFile = sHeadModel.SurfaceFile;
574 % ResultsMat3.GridLoc = GridLoc;
575 %
576 % % === NOT SHARED ===
577 % % Get the output study (pick the one from the first file)
578 % iStudy = sInputs(1).iStudy;
579 % % Create a default output filename
580 % OutputFiles{3} = bst_process('GetNewFilename', fileparts(sInputs(1).FileName), 'results_MCB_oriDiff');
581 % % Save on disk
582 % save(OutputFiles{3}, '-struct', 'ResultsMat3');
583 % % Register in database
584 % db_add_data(iStudy, OutputFiles{3}, ResultsMat3);
585 % end
586
587 end
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.