-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathart.m
2144 lines (1960 loc) · 88.1 KB
/
art.m
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
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
function varargout = art(varargin)
% art - module for automatic and manual detection and removal of outliers.
%
% art
% Lauches the art GUI. This will prompt the user to enter the functional
% volumes and motion pararameters (SPM, FSL, or Siemens supported) for one
% subject (one or multiple sessions). It then displays four graphs:
%
% The top graph is the global brain activation mean as a function of
% time, with all of the identified outlier scans marked. Optionally it also
% overlais the task-related regressors (SPM designs only) and break down of
% number of identified outliers per condition.
%
% The second graph shows the timeseries derived from the global BOLD signal
% and used to identify potential outlier scans (either absolute or
% scan-to-scan differences in the global BOLD signal normalized to
% z-scores)
%
% The third graph shows the timeseries derived from the subject motion
% parameters and used to identify potential outlier scans (either the
% absolute or scan-to-scan differences in the individual motion parameters
% -x,y,z translation parameters and roll,pitch,yaw rotation parameters-, or
% the absolute or scan-to-scan differences in a single 'composite motion'
% measure -maximal movement of any voxel within the brain bounding box-)
%
% Using default threshold values for each of the bottom three graphs
% we define outliers as points that exceed the threshold in at least
% one of the global signal or motion parameter graphs. The thresholds are
% shown as horizontal black lines in each of the graphs.
%
% Points which are identified as outliers, are indicated by a vertical
% black line in the graph that corresponds to the outlying
% parameter(s). For example, the if the absolute value of the Y motion
% parameter for time t=17 is above the motion threshold, it is
% identified as an outlier and indicated by a black vertical line at
% t=17 in the third graph. The union of all outliers is indicated by
% red vertical lines on the top graph. The list of outliers is also
% displayed in the editable text box below the graphs. The current
% values of the thresholds are displayed by the side of the
% corresponding graphs. These values may can be changed by the user
% either by pressing the up/down buttons, which increment/decrement
% the current value by 10%, or by specifying a new value in the text
% box.
%
% In Addition, the user can manually add or remove points from the
% list of outliers by editting the list. Note that the list is only
% updated once the curser points outside the text box (i.e. click the
% mouse somewhere outside the text box). Since any changes made by the
% user are overridden once the thresholds are updated, it is
% recommended to do any manual changes as the last step before saving.
%
% In the 'options' section of the GUI, the user can select to display the
% task-related design, the spectra of the global BOLD signal motion or
% task-related parameters, the level of task-correlated motion and
% task-correlated global BOLD signal, as well as the corrected analysis
% mask (without the influence of the identified outliers)
%
% By default art generates the following output files:
% Regressor files (one per session, stored in the same folders as the
% functional volumes, and named art_regression_outliers_*.mat and
% art_regression_outliers_and_movement_*.mat). This files can be entered
% as covariates in the first-level analyses in order to effectively
% remove the identified outlier scans from further analyses
% Analysis mask (one file named art_mask.img defining the analysis mask
% after disregarding outlier scans). This file can be entered as an
% explicit analysis mask in the first-level analyses in order to
% avoid any influences of the outlier scans on the implicit analysis mask
% computation (on SPM you will also need to modify the defaults in order
% to skip the implicit masking operation, e.g. set defaults.mask.thresh =
% -inf)
%
% Pressing the save button lets the user choose wheter to save the
% list of identified outliers, motion statistics, graphs, outlier
% regressors, or Analysis mask.
%
%
% art filename.cfg
% Uses the configuration file filename.cfg to define the art analysis
% options (see example.cfg file for more information)
%
% See also art_batch
%
% ----------------------------------------------------------------------
% - Added voxel-wise SNR and variability displays; save composite motion
% measure; clean-up code -Alfonso 06/11
% - Added analysis mask option, modified pow spec display and other GUI display
% options, last update 9/25/09 - Sue, Alfonso and Darren
% - if multiple sessions are specified, standard deviations are calculated
% within sessions
% oliver hinds 2008-04-23
% - added support for reading siemens motion paramter file format
% oliver hinds 2008-04-23
% - added ability to read file names and sesions from config file
% oliver hinds 2008-04-23
% - tiny fix to make Matlab 6.5 compatible - Shay 5/14/07
% - added "signal-task correlation" - 5/11/07
% - added "motion-task correlation" and "show spectrum"
% - minor GUI changes to support Windows and open a large graph
% in a separate window. Also fixed starnge motion params filename
% bug. Shay Mozes, 5/2/2007
% - fixed bug in display of motion outlier on the graph, Shay Mozes 4/30/2007
% - superimpose task conditions on the z-graph, Shay Mozes, 4/24/2007
% - added support for SPM5, Shay Mozes, 4/2007
% - now supporting FSL .par format, 4/9/2007
% - new GUI and features Shay Mozes, 2006
% + Mar. 2007 from art_global.m, by Paul Mazaika, April 2004.
% from artdetect4.m, by Jeff Cooper, Nov. 2002
% from artdetect3.m, by Sue Whitfield artdetect.m
% Sue Whitfield 2000
%% --------------- GUI initialization -----------------------------
% GUIDE GUI default creation and callback support
% --------------------------------------------------------------------
gui_Singleton = 0;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @art_OpeningFcn, ...
'gui_OutputFcn', @art_OutputFcn, ...
'gui_LayoutFcn', [] , ...
'gui_Callback', []);
if nargin && ischar(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end
if nargout
[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
end
%% --------------- Main process -----------------------------
% Main functionality of art (this process runs once for each dataset)
% --------------------------------------------------------------------
% -----------------------------------------------------------------------
% ART_OPENINGFCN
% This function is called before opening the gui.
% It executes the main functionality of art, extracting the global signal
% and movement parameters, defining the measures used to identify outliers,
% printing movement statistics, and calling the Update* functions to
% identify outliers and generate the corresponding plots
% -----------------------------------------------------------------------
function art_OpeningFcn(hObject, eventdata, handles, varargin)
% -----------------------
%Initialize
% -----------------------
warning('OFF','all')
%find spm version
if isdeployed, spm_ver='SPM12';
else spm_ver = spm('ver');
end
switch(spm_ver),
case 'SPM99', spm_ver=1;
case 'SPM2', spm_ver=2;
case 'SPM5', spm_ver=5;
case {'SPM8b','SPM8'}, spm_ver=8;
case {'SPM12b','SPM12'}, spm_ver=12;
otherwise, disp(['Warning! unrecognized SPM version ',spm_ver]); spm_ver=12;
end
%clear data from previous sessions
try
setappdata(handles.showDesign,'SPM',[]);
setappdata(handles.showDesign,'sessions',[]);
setappdata(handles.showDesign,'SPMfile',[]);
setappdata(handles.zthresh,'g',[]);
setappdata(handles.mvthresh,'mv_data',[]);
setappdata(handles.mvthresh,'drop_flag',[]);
setappdata(handles.mvthresh,'altval',[]);
setappdata(handles.rtthresh,'altval',[]);
setappdata(handles.savefile,'path',[]);
setappdata(handles.savefile,'datafiles',[]);
setappdata(handles.savefile,'stats_file',[]);
setappdata(handles.savefile,'analyses',[]);
setappdata(handles.mvthresh,'mv_stats',[]);
setappdata(handles.zthresh,'zoutliers',[]);
setappdata(handles.mvthresh,'mv_norm_outliers',[]);
setappdata(handles.mvthresh,'mv_x_outliers',[]);
setappdata(handles.mvthresh,'mv_y_outliers',[]);
setappdata(handles.mvthresh,'mv_z_outliers',[]);
setappdata(handles.rtthresh,'rt_norm_outliers',[]);
setappdata(handles.rtthresh,'rt_p_outliers',[]);
setappdata(handles.rtthresh,'rt_r_outliers',[]);
setappdata(handles.rtthresh,'rt_y_outliers',[]);
catch %#ok<*CTCH>
end
% ------------------------
% Default values for outliers
% ------------------------
z_thresh = 1.5; %global signal threshold
mvmt_thresh = 2.0; %absolute subject motion threshold
rotat_thresh = .05; %absolute subject rotation threshold
mvmt_diff_thresh = 2.0; %scan-to-scan subject motion threshold
rotat_diff_thresh = .02; %scan-to-scan subject rotation threshold
sess_file='';
stats_file='';
output_dir='';
%Passo i parametri del for loop (i e session file)
i = evalin('base','i');
Session_name= evalin('base','Session_name');
out_ppt=evalin('base','folder_out');
output_path=[out_ppt,'/',Session_name{i,1}];
output_dir=output_path;
global path_general;
% look for args in varargin,
for i=1:2:numel(varargin)
if strcmp(varargin{i}, 'sess_file')
sess_file = varargin{i+1};
elseif strcmp(varargin{i}, 'stats_file')
stats_file = varargin{i+1};
elseif strcmp(varargin{i}, 'output_path')
output_path = varargin{i+1};
elseif i==numel(varargin),
if exist(varargin{i},'file'),
[filepath,filename,fileext]=fileparts(varargin{i});
if strcmp(fileext,'.cfg'), sess_file=varargin{i}; end
end
end
end
setappdata(handles.savefile,'stats_file',stats_file);
% ------------------------
% Collect files
% ------------------------
if ~isempty(sess_file) % read config
[num_sess,global_type_flag,drop_flag,motionFileType,motion_threshold,global_threshold,use_diff_motion,use_diff_global,use_norms,SPMfile,mask_file,output_dir,P,M] = ...
read_art_sess_file(sess_file);
% else
% motion_threshold=[];global_threshold=[];SPMfile=[];mask_file=[];output_dir='';
% use_diff_motion=0;use_diff_global=0;use_norms=1;
% num_sess = spm_input('How many sessions?',1,'n',1,1);
%
% global_type_flag = spm_input('Which global mean to use?', 1, 'm', ...
% 'Regular | User Mask',...
% [1 2], 1);
% motionFileType = spm_input('Select type of motion params file.',1,'m',...
% ' txt(SPM) | par(FSL) | txt(Siemens)', ...
% [0 1 2], 0);
%
% P=cell(1,num_sess);
% M=cell(1,num_sess);
% for i = 1:num_sess
% switch spm_ver
% case {1,2}
% P{i} = spm_get(Inf,'.img',['Select functional volumes for session' num2str(i) ':']);
% case {5,8,12}
% P{i} = spm_select(Inf,'image',['Select functional volumes for session' num2str(i) ':']);
% %P{i} = spm_select(Inf,'.*\.nii|.*\.img',['Select functional volumes for session' num2str(i) ':']);
% end
% if motionFileType == 0 %SPM format
% switch spm_ver
% case {1,2}
% mvmt_file = spm_get(1,'.txt',['Select movement params file for session' num2str(i) ':']);
% case {5,8,12}
% mvmt_file = spm_select(1,'^.*\.txt$',['Select movement params file for session' num2str(i) ':']);
% end
% M{i} =load(mvmt_file);
% elseif motionFileType == 1 %FSL format
% switch spm_ver
% case {1,2}
% mvmt_file = spm_get(1,'.par',['Select movement params file for session' num2str(i) ':']);
% case {5,8,12}
% mvmt_file = spm_select(1,'^.*\.par$',['Select movement params file for session' num2str(i) ':']);
% end
% M{i} =load(mvmt_file);
% elseif motionFileType == 2 % Siemens MotionDetectionParameter.txt
% switch spm_ver
% case {1,2}
% mvmt_file = spm_get(1,'.txt',['Select movement params file for session' num2str(i) ':']);
% case {5,8,12}
% mvmt_file = spm_select(1,'^.*\.txt$',['Select movement params file for session' num2str(i) ':']);
% end
% M{i} = read_siemens_motion_parm_file(mvmt_file);
% end
% output_path = fileparts(mvmt_file);
% end
% drop_flag = 0;
end
if global_type_flag==2,
if isempty(mask_file),
mask_file = spm_select(1, '.*\.nii|.*\.img', 'Select mask image in functional space');
end
mask=spm_vol(mask_file);
end
setappdata(handles.showDesign,'sessions',-(1:num_sess));
setappdata(handles.showDesign,'SPMfile',SPMfile);
if ~isempty(SPMfile), temp=load(SPMfile); setappdata(handles.showDesign,'SPM',temp.SPM); clear temp; end
datafiles=cell(1,length(P));
for i=1:length(P),
if ~isempty(P{i}),
datafiles{i}=strtrim(P{i}(1,:));
else
datafiles{i}=[];
end;
end % <alfnie>: keep filenames of functional data (first scan per session only)
if ~isempty(motion_threshold),
mvmt_thresh=motion_threshold(1); mvmt_diff_thresh=motion_threshold(1);
if numel(motion_threshold)>1, rotat_thresh=motion_threshold(2); rotat_diff_thresh=motion_threshold(2); end
end
if ~isempty(global_threshold), z_thresh=global_threshold; end
%if drop_flag, disp('warning: explicitly dropping 1st scan no longer supported. Edit the ''all outliers'' box to specify any number of additional outlier scans'); end
% dropflag extended to add first N scans as outliers in savefile_Callback_SaveRegressor (alfnie 09/15)
mv_data = [];
for i = 1:length(M)
mv_data = vertcat(mv_data,M{i});
end
%translate to SPM format: x y z (in mm) pitch roll yaw (in radians)
if motionFileType == 1, %FSL order of fields is three Euler angles (x,y,z in radians) then three translation params (x,y,z in mm).
tmp = mv_data(:,1:3);
mv_data(:,1:3) = mv_data(:,4:6);
mv_data(:,4:6) = tmp;
elseif motionFileType == 3, %x y z (in mm) pitch roll jaw (in degrees)
mv_data(:,4:6)=mv_data(:,4:6)/180*pi;
end
% NEW VERSION
% ---------------------------------------
% Compute Global signal and analysis mask
% ---------------------------------------
g = cell(1,num_sess); % g is a cell array of the global mean for each scan in each session
gsigma=g;gmean=g;dgsigma=g;dgmean=g;
maskscan={};
maskscan_files=0;
art_mask_temporalfile=['art_mask_temporalfile',datestr(now,'yyyymmddHHMMSSFFF'),'.mat'];
Data_Sum=0;
Data_SumSquared=0;
VY=cell(1,num_sess);
cumdisp;
notcoregistered=false;
for sess=1:num_sess
fprintf('%-4s: ',['Mapping files for session ' num2str(sess) '...']);
VY{sess} = spm_vol(P{sess});
fprintf('%3s\n','...done')
switch spm_ver
case {1,2}
if any(any(diff(cat(1,VY{sess}.dim),1,1),1)&[1,1,1,0])
error('images do not all have the same dimensions')
end
case {5,8,12}
if any(any(diff(cat(1,VY{sess}.dim),1,1),1))
error('images do not all have the same dimensions')
end
end
if ~notcoregistered
for n=1:numel(VY{sess}),
if any(any(VY{sess}(n).mat~=VY{1}(1).mat)), notcoregistered=true; break; end
end
end
% --------------------------------------------------
% Extracts global-signal mask and (optionally) compute analysis mask
% note: scan-specific global-signal mask = voxels above mean(all voxels)/8
% scan-specific analysis mask = voxels above 0.8*mean(global-signal mask)
% global-signal mask = intersection of {scan-specific global-signal mask}
% analysis mask = intersection of {scan-specific analysis mask}
% --------------------------------------------------
nscans = numel(VY{sess});
g{sess} = zeros(nscans,4);
fprintf('%-4s: %3s','Calculating globals...',' ')
if sess==1
VY1inv=pinv(VY{sess}(1).mat);
[tempx,tempy,tempz]=ind2sub(VY{sess}(1).dim(1:3),1:prod(VY{sess}(1).dim(1:3)));
xyz_voxel=[tempx(:),tempy(:),tempz(:),ones(numel(tempx),1)]';
xyz=VY{sess}(1).mat*xyz_voxel;
end
if global_type_flag==1 % regular mean : Global-conjunction (uses conjunction of individual scan masks; individual scan mask are defined as voxels above mean/8 for each scan; see art_maskglobal_scan)
Mask=ones(VY{sess}(1).dim(1:3));
for i = 1:nscans,
if notcoregistered, temp=spm_get_data(VY{sess}(i),pinv(VY{sess}(i).mat)*xyz);
else temp=reshape(spm_get_data(VY{sess}(i),xyz_voxel),VY{sess}(i).dim); %temp=spm_read_vols(VY{sess}(i));
end
[maskscan{end+1},masktemp]=art_maskglobal_scan(temp,VY{sess}(i),VY{sess}(1),VY1inv); %#ok<AGROW>
if numel(maskscan)>=100, % stop maskcan from growing too much
maskscan_files=maskscan_files+1;
try
save(fullfile(output_dir,[art_mask_temporalfile(1:end-4),num2str(maskscan_files),'.mat']),'maskscan');
catch
fprintf('warning: unable to write to ',output_dir,' folder. Writing output files to ',pwd,' instead.');
output_dir=pwd;
save(fullfile(output_dir,[art_mask_temporalfile(1:end-4),num2str(maskscan_files),'.mat']),'maskscan');
end
maskscan={};
end
Mask(masktemp)=0;
cumdisp([num2str(i),'/',num2str(nscans)]);
end
elseif global_type_flag==2 % user-defined mask
Mask=spm_get_data(mask,pinv(mask.mat)*xyz);
end
idxMask=find(Mask);
% --------------------------------------------------
% computes global signal
% --------------------------------------------------
badmask=(numel(idxMask)<numel(Mask)/10) && global_type_flag~=2;
for i = 1:nscans,
if notcoregistered, temp=reshape(spm_get_data(VY{sess}(i),pinv(VY{sess}(i).mat)*xyz),VY{1}(1).dim);
else temp=reshape(spm_get_data(VY{sess}(i),xyz_voxel),VY{sess}(i).dim); %temp=spm_read_vols(VY{sess}(i));
end
if badmask, g{sess}(i) = spm_global(VY{sess}(i));
else g{sess}(i)=mean(temp(idxMask));
end
Data_Sum=Data_Sum+temp;
Data_SumSquared=Data_SumSquared+temp.^2;
%imagesc(Data_Sum(:,:,30));colorbar;drawnow;
end
fprintf('...done');cumdisp;
% OLD VERSION
% % ---------------------------------------
% % Compute Global signal and analysis mask
% % ---------------------------------------
% g = cell(1,num_sess); % g is a cell array of the global mean for each scan in each session
% gsigma=g;gmean=g;dgsigma=g;dgmean=g;
% maskscan={};
% Data_Sum=0;
% Data_SumSquared=0;
% VY=cell(1,num_sess);
% cumdisp;
% for sess=1:num_sess
% fprintf('%-4s: ',['Mapping files for session ' num2str(sess) '...']);
% VY{sess} = spm_vol(P{sess});
% fprintf('%3s\n','...done')
%
% switch spm_ver
% case {1,2}
% if any(any(diff(cat(1,VY{sess}.dim),1,1),1)&[1,1,1,0])
% error('images do not all have the same dimensions')
% end
% case {5,8,12}
% if any(any(diff(cat(1,VY{sess}.dim),1,1),1))
% error('images do not all have the same dimensions')
% end
% end
%
% % --------------------------------------------------
% % Extracts global-signal mask and (optionally) compute analysis mask
% % note: scan-specific global-signal mask = voxels above mean(all voxels)/8
% % scan-specific analysis mask = voxels above 0.8*mean(global-signal mask)
% % global-signal mask = intersection of {scan-specific global-signal mask}
% % analysis mask = intersection of {scan-specific analysis mask}
% % --------------------------------------------------
% nscans = numel(VY{sess});
% g{sess} = zeros(nscans,4);
% fprintf('%-4s: %3s','Calculating globals...',' ')
% VY1inv=pinv(VY{sess}(1).mat);
% [tempx,tempy,tempz]=ind2sub(VY{sess}(1).dim(1:3),1:prod(VY{sess}(1).dim(1:3)));
% xyz_voxel=[tempx(:),tempy(:),tempz(:),ones(numel(tempx),1)]';
% xyz=VY{sess}(1).mat*xyz_voxel;
% if global_type_flag==1 % regular mean : Global-conjunction (uses conjunction of individual scan masks; individual scan mask are defined as voxels above mean/8 for each scan; see art_maskglobal_scan)
% Mask=ones(VY{sess}(1).dim(1:3));
% for i = 1:nscans,
% temp=reshape(spm_get_data(VY{sess}(i),xyz_voxel),VY{sess}(i).dim);
% [maskscan{end+1},masktemp]=art_maskglobal_scan(temp,VY{sess}(i),VY{sess}(1),VY1inv); %#ok<AGROW>
% Mask(masktemp)=0;
% cumdisp([num2str(i),'/',num2str(nscans)]);
% end
% elseif global_type_flag==2 % user-defined mask
% Mask=spm_get_data(mask,pinv(mask.mat)*xyz);
% end
% idxMask=find(Mask);
%
% % --------------------------------------------------
% % computes global signal
% % --------------------------------------------------
% notcoregistered=length(VY{sess})>1&&any(any(VY{sess}(2).mat~=VY{sess}(1).mat));
% badmask=(numel(idxMask)<numel(Mask)/10) && global_type_flag~=2;
% for i = 1:nscans,
% if notcoregistered, temp=spm_get_data(VY{sess}(i),pinv(VY{sess}(i).mat)*xyz);
% else temp=reshape(spm_get_data(VY{sess}(i),xyz_voxel),VY{sess}(i).dim); %temp=spm_read_vols(VY{sess}(i));
% end
% if badmask, g{sess}(i) = spm_global(VY{sess}(i));
% else g{sess}(i)=mean(temp(idxMask));
% end
% Data_Sum=Data_Sum+temp;
% Data_SumSquared=Data_SumSquared+temp.^2;
% %imagesc(Data_Sum(:,:,30));colorbar;drawnow;
% end
% fprintf('...done');cumdisp;
% --------------------------------------------------
% Compute derived signals for outlier identification
% --------------------------------------------------
% g{sess} columns:
% 1: global signal
% 2: standardized global signal
% 3: scan-to-scan differences in global signal
% 4: standardized scan-to-scan differences in global signal
gsigma{sess} = .7413*diff(prctile(g{sess}(:,1),[25,75]));gsigma{sess}(gsigma{sess}==0)=1; % robus standard-deviation
gmean{sess} = median(g{sess}(:,1)); % robust mean
g{sess}(:,2)=(g{sess}(:,1)-gmean{sess})/max(eps,gsigma{sess}); % z-score
g{sess}(2:end,3)=diff(g{sess}(:,1),1,1);
dgsigma{sess} = .7413*diff(prctile(g{sess}(:,3),[25,75]));dgsigma{sess}(dgsigma{sess}==0)=1;
dgmean{sess} = median(g{sess}(:,3));
g{sess}(2:end,4)=(g{sess}(2:end,3)-dgmean{sess})/max(eps,dgsigma{sess});
z_thresh = 0.1*round(z_thresh*10);
end
art_mask_temporalfile=['art_mask_temporalfile',char('0'+floor(10*rand(1,8))),'.mat'];
VY=cat(1,VY{:}); VY1=VY(1); %#ok<NASGU>
try
save(fullfile(output_dir,art_mask_temporalfile),'maskscan','VY1','VY','notcoregistered','xyz','xyz_voxel','Data_Sum','Data_SumSquared');
catch
disp('warning: unable to write to ',output_dir,' folder. Writing output files to ',pwd,' instead.');
output_dir=pwd;
save(fullfile(output_dir,art_mask_temporalfile),'maskscan','VY1','VY','notcoregistered','xyz','xyz_voxel','Data_Sum','Data_SumSquared');
end
set(handles.figure1,'closerequestfcn',['try,if ispc,[nill,ok]=system(''del "',fullfile(output_dir,art_mask_temporalfile),'"'');else [nill,ok]=system(''rm ''''',fullfile(output_dir,art_mask_temporalfile),'''''''); end; end; delete(gcbf);']);
%update text fields
set(handles.data_stdv,'String',num2str(cat(2,gsigma{:}),'%0.1f '));
set(handles.zthresh,'String',num2str(z_thresh));
set(handles.mvthresh,'String',num2str(mvmt_thresh));
set(handles.rtthresh,'String',num2str(rotat_thresh));
% ------------------------------------------------------------------------
% Compute Movement parameters and derived signals for outlier identification
% ------------------------------------------------------------------------
% mv_data columns:
% 1-6 : raw movement parameters
% 7 : euclidean norm of raw movement parameters
% 8-13 : scan-to-scan differences in raw movement parameters
% 14-31 : 6 control points trajectories (placed on center of faces of bounding box; x,y,z coordinates for each control point)
% 32 : composite measure: euclidean norm of control point trajectories
% 33-50 : scan-to-scan differences in 6 control points trajectories
% 51 : composite measure: max scan-to-scan movement across the 6 control points
mv_data=[mv_data,zeros([size(mv_data,1),51-size(mv_data,2)])];
respos=diag([70,70,75]);resneg=diag([-70,-110,-45]);
res=[respos,zeros(3,1),zeros(3,4),zeros(3,4),eye(3),zeros(3,1); % 6 control points: [+x,+y,+z,-x,-y,-z];
zeros(3,4),respos,zeros(3,1),zeros(3,4),eye(3),zeros(3,1);
zeros(3,4),zeros(3,4),respos,zeros(3,1),eye(3),zeros(3,1);
resneg,zeros(3,1),zeros(3,4),zeros(3,4),eye(3),zeros(3,1);
zeros(3,4),resneg,zeros(3,1),zeros(3,4),eye(3),zeros(3,1);
zeros(3,4),zeros(3,4),resneg,zeros(3,1),eye(3),zeros(3,1);];
for i=1:size(mv_data,1)
temp=spm_matrix([1*mv_data(i,1:3),mv_data(i,4:6)]); temp=temp(:)';
mv_data(i,14:31)=temp*res';
end
cur_sess_start=0;
for sess=1:num_sess
n=length(g{sess}(:,1));
mv_data(cur_sess_start+(1:n),7) = sqrt(sum(abs(mv_data(cur_sess_start+(1:n),1:3)).^2,2));
mv_data(cur_sess_start+(2:n),8:13) = diff(mv_data(cur_sess_start+(1:n),1:6),1,1);
mv_data(cur_sess_start+(1:n),32)=sqrt(mean(abs(detrend(mv_data(cur_sess_start+(1:n),14:31),'constant')).^2,2));
mv_data(cur_sess_start+(2:n),33:50)=diff(mv_data(cur_sess_start+(1:n),14:31),1,1);
mv_data(cur_sess_start+(2:n),51)=max(sqrt(sum(reshape(abs(mv_data(cur_sess_start+(2:n),33:50)).^2,[n-1,3,6]),2)),[],3);
cur_sess_start = cur_sess_start + n;
end
%save application data for use in callbacks
setappdata(handles.zthresh,'g',g);
setappdata(handles.mvthresh,'mv_data',mv_data);
setappdata(handles.mvthresh,'altval',num2str(mvmt_diff_thresh));
setappdata(handles.rtthresh,'altval',num2str(rotat_diff_thresh));
setappdata(handles.mvthresh,'drop_flag',drop_flag);
setappdata(handles.savefile,'path',output_path);
setappdata(handles.savefile,'dir',output_dir);
setappdata(handles.savefile,'datafiles',datafiles);
setappdata(handles.savefile,'spm_ver',spm_ver);
setappdata(handles.savefile,'mv_data_raw',M);
setappdata(handles.savefile,'art_mask_temporalfile',art_mask_temporalfile);
if ~isempty(SPMfile), set(handles.showDesign,'Value',get(handles.showDesign,'Max')); end
set(handles.norms,'Value',use_norms);
set(handles.diff1,'value',use_diff_global);
set(handles.diff2,'value',use_diff_motion);
if use_diff_global||use_diff_motion
diffsglobalandmotion_Callback(hObject, [], handles, 2*use_diff_global-1, 2*use_diff_motion-1);
else
UpdateGlobal(hObject, eventdata, handles,1.0);
UpdateMovement(hObject, eventdata, handles,1.0);
UpdateRotation(hObject, eventdata, handles,1.0);
UnionOutliers(hObject, eventdata, handles);
UpdateSummaryplot(hObject, eventdata, handles);
end
idx=str2num(get(handles.all_outliers, 'String')); %#ok<*ST2NM>
for sess=1:num_sess
fprintf('\nSession %d global statistics - mean: %7.4f stdv: %7.4f',sess,gmean{sess},gsigma{sess});
end
fprintf('\n');
fprintf('Outlier detection: %d identified outliers\n',length(idx));
% ------------------------
% Compute and print statistics of movement
%--------------------------
mv_data = getappdata(handles.mvthresh,'mv_data');
mv_stats = [mean(abs(mv_data)); std(abs(mv_data)); max(abs(mv_data)) ];
setappdata(handles.mvthresh,'mv_stats',mv_stats);
fprintf('\n\nStatistics of movement data:\n\n');
fprintf('%5s%10s%10s%10s%11s%10s%9s%10s\n',' ','x','y','z',' pitch','roll','yaw','norm');
fprintf('%7s%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f\n','mean ',mv_stats(1,1:3),mv_stats(1,4:6),mv_stats(1,32));
fprintf('%7s%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f\n','stdv ',mv_stats(2,1:3),mv_stats(2,4:6),mv_stats(2,32));
fprintf('%7s%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f\n\n\n','max ',mv_stats(3,1:3),mv_stats(3,4:6),mv_stats(3,32));
% BEGIN ohinds 2008-04-23: save stats to file
if ~isempty(stats_file)
if isempty(fileparts(stats_file)),stats_file=fullfile(output_dir,stats_file);end
fp = fopen(stats_file,'w');
if fp ~= -1
fprintf('saving global motion stats to %s\n',stats_file);
fprintf(fp,'%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f\n',mv_stats(1,1:3),mv_stats(1,4:6),mv_stats(1,32));
fprintf(fp,'%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f\n',mv_stats(2,1:3),mv_stats(2,4:6),mv_stats(2,32));
fprintf(fp,'%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f\n\n\n',mv_stats(3,1:3),mv_stats(3,4:6),mv_stats(3,32));
fclose(fp);
end
end
% END ohinds 2008-04-23: save stats to file
% Choose default command line output for art
handles.output = hObject;
% Update handles structure
guidata(hObject, handles);
set(handles.figure1,'resizefcn','art(''showOptions_Callback'',gcbo,[],guidata(gcbo))');
% --------------------------------------------------------
% Saves regressor matrix with outliers & new analysis mask
%---------------------------------------------------------
savefile_Callback_SaveRegressor(handles);
try
savefile_Callback_SaveMask(handles);
catch
disp('Warning. Error encountered during ART implicit-mask file creation. Skipping this step');
end
try
nome_path_ptovv=[path_general,'M4_BS_02_Scrub/art_single_sessions/'];
%nome_path_ptovv='/data/00_rsFetMRI_PreProc_1.5T/Dati_fmri/M4_BS_02_Scrub/art_single_sessions/';
%save outliers
out_idx = round(str2num(get(handles.all_outliers, 'String')));
filename='art_outliers.mat';
filename = strcat(char(nome_path_ptovv), char(filename));
save(filename,'out_idx', '-mat');
%save motion_statistics
mv_stats = getappdata(handles.mvthresh,'mv_stats');
filename='art_motion_statistics.mat';
filename = strcat(char(nome_path_ptovv), char(filename));
save(filename ,'mv_stats','-mat');
% save figure
filename='art_figure.jpg';
filename = strcat(char(nome_path_ptovv), char(filename));
saveas(gcf,filename);
catch
disp(':Error encountered during ART file creation. Skipping this step');
end
try % fix Matlab14 new graphic objects issue (uipanel shown on top of uicontrol objects)
h1=findobj(handles.figure1,'type','uipanel');
h2=findobj(handles.figure1,'type','uicontrol');
x1=cell2mat(get(h1,'position'));
x2=cell2mat(get(h2,'position'));
for n=1:size(x2,1)
xpos=repmat(x2(n,1:2),[size(x1,1),1]);
[md,n0]=max(sign(min(min(xpos-x1(:,1:2),x1(:,1:2)+x1(:,3:4)-xpos),[],2)),[],1);
if md>0, set(h2(n),'parent',h1(n0)); end
end
set(h2,'parent',handles.figure1);
end
end
%% --------------- GUI callback functions --------------------------
% Processing steps for each of the art GUI options
% --------------------------------------------------------------------
% -----------------------------------------------------------------------
% UPDATEGLOBAL
% This function identifies outliers based on the BOLD globalsignal and
% generates the corresponding plot (gui 2nd plot from the top)
% -----------------------------------------------------------------------
function UpdateGlobal(hObject, eventdata, handles, incr) %#ok<INUSL>
%get data
z_thresh = str2num(get(handles.zthresh,'String'));
g = getappdata(handles.zthresh,'g');
num_sess = length(g);
%calc new outliers
% BEGIN ohinds 2008-04-23: plot zscores
axes(handles.zvalue);
cla;
hold on;
cur_sess_start=1;
z_thresh = z_thresh*incr;
idxind=2;
out_idx = cell(1,num_sess);
for sess=1:num_sess
if get(handles.diff1,'value'),
out_idx{sess} = cur_sess_start+(find(abs(g{sess}(:,idxind)) > z_thresh|abs([g{sess}(2:end,idxind);0]) > z_thresh))'-1;
else
out_idx{sess} = cur_sess_start+(find(abs(g{sess}(:,idxind)) > z_thresh))'-1;
end
%update plot
plot(cur_sess_start:cur_sess_start+size(g{sess},1)-1, g{sess}(:,idxind));
cur_sess_start = cur_sess_start + size(g{sess},1);
end
out_idx=cat(2,out_idx{:});
set(gca,'xlim',[0,cur_sess_start]);
l=ylabel('global mean\newline [std]');
set(l,'VerticalAlignment','bottom','horizontalalignment','center');
set(gca,'XTickLabel',[]);
thresh_x = 1:cur_sess_start-1;
thresh_y = z_thresh*ones(1,length(thresh_x));
line(thresh_x, thresh_y, 'Color', 'black');
line(thresh_x, -1*thresh_y, 'Color', 'black');
%update text
set(handles.zthresh,'String',num2str(z_thresh));
setappdata(handles.zthresh,'zoutliers',out_idx);
hold off;
% END ohinds 2008-04-23: plot zscores
%plot outliers
axes_lim = get(gca, 'YLim');
axes_height = axes_lim;
for i = 1:length(out_idx)
line((out_idx(i)*ones(1, length(axes_height))), axes_height, 'Color', 'black');
end
end
% -----------------------------------------------------------------------
% UPDATEMOVEMENT
% This function identifies outliers based on the subject movement
% parameters (either translation parameters or composite motion) and
% generates the corresponding plot (gui 3rd plot from the top)
% -----------------------------------------------------------------------
function UpdateMovement(hObject, eventdata, handles, incr) %#ok<INUSL>
%get data
mvmt_thresh = str2num(get(handles.mvthresh,'String'));
mv_data = getappdata(handles.mvthresh,'mv_data');
%calc new outliers
mvmt_thresh = mvmt_thresh*incr;
if get(handles.diff2,'value'),
out_mvmt_idx = (find(abs(mv_data(:,1:3)) > mvmt_thresh | [abs(mv_data(2:end,1:3));[0,0,0]] > mvmt_thresh ))';
else
out_mvmt_idx = (find(abs(mv_data(:,1:3)) > mvmt_thresh))';
end
out_mvmt_idx_X=out_mvmt_idx(out_mvmt_idx<=size(mv_data,1));
out_mvmt_idx_Y=out_mvmt_idx(out_mvmt_idx>size(mv_data,1)&out_mvmt_idx<=2*size(mv_data,1))-size(mv_data,1);
out_mvmt_idx_Z=out_mvmt_idx(out_mvmt_idx>2*size(mv_data,1))-2*size(mv_data,1);
%find norm outliers
if (get(handles.norms,'Value') == get(handles.norms,'Max'))
normv=mv_data(:,32);
if get(handles.diff2,'value')
out_mvmt_idx_norm = find(normv>mvmt_thresh|[normv(2:end);0]>mvmt_thresh);
else
out_mvmt_idx_norm = find(normv>mvmt_thresh);
end
out_mvmt_idx_norm = out_mvmt_idx_norm';
setappdata(handles.mvthresh,'mv_norm_outliers',out_mvmt_idx_norm);
end
%update text
set(handles.mvthresh,'String',num2str(mvmt_thresh));
setappdata(handles.mvthresh,'mv_x_outliers',out_mvmt_idx_X);
setappdata(handles.mvthresh,'mv_y_outliers',out_mvmt_idx_Y);
setappdata(handles.mvthresh,'mv_z_outliers',out_mvmt_idx_Z);
axes(handles.mvmtGraph);
cla;
if (get(handles.norms,'Value') == get(handles.norms,'Max'))
plot(normv);
set(gca,'xlim',[0,length(normv)+1]);
else
plot(mv_data(:,1:3));
set(gca,'xlim',[0,size(mv_data,1)+1]);
end
l=ylabel('movement \newline [mm]');
set(l,'VerticalAlignment','bottom','horizontalalignment','center');
set(gca,'XTickLabel',[]);
if (get(handles.norms,'Value') ~= get(handles.norms,'Max'))
legend('x', 'y', 'z','Location','East');
end
h = gca;
set(h,'Ygrid','on');
thresh_mv_x = 1:size(mv_data,1);
thresh_mv_y = mvmt_thresh*ones(1,size(mv_data,1));
line(thresh_mv_x, thresh_mv_y, 'Color', 'black');
if ~(get(handles.norms,'Value') == get(handles.norms,'Max'))
line(thresh_mv_x, -1*thresh_mv_y, 'Color', 'black');
end
axes_lim = get(gca, 'YLim');
axes_height = axes_lim;
if ~(get(handles.norms,'Value') == get(handles.norms,'Max'))
for i = 1:length(out_mvmt_idx_X)
line((out_mvmt_idx_X(i)*ones(1, length(axes_height))), axes_height, 'Color', 'black');
end
for i = 1:length(out_mvmt_idx_Y)
line((out_mvmt_idx_Y(i)*ones(1, length(axes_height))), axes_height, 'Color', 'black');
end
for i = 1:length(out_mvmt_idx_Z)
line((out_mvmt_idx_Z(i)*ones(1, length(axes_height))), axes_height, 'Color', 'black');
end
else
for i = 1:length(out_mvmt_idx_norm)
line((out_mvmt_idx_norm(i)*ones(1, length(axes_height))), axes_height, 'Color', 'black');
end
end
end
% -----------------------------------------------------------------------
% UPDATEROTATION
% This function identifies outliers based on the subject rotation
% parameters and generates the corresponding plot (gui 4th plot from the top)
% (note: only available when not using composite movement measures)
% -----------------------------------------------------------------------
function UpdateRotation(hObject, eventdata, handles, incr) %#ok<INUSL>
%get data
rotat_thresh = str2num(get(handles.rtthresh,'String'));
mv_data = getappdata(handles.mvthresh,'mv_data');
%calc new outliers
rotat_thresh = rotat_thresh*incr;
out_rotat_idx = (find(abs(mv_data(:,4:6)) > rotat_thresh))';
out_rotat_idx_p=out_rotat_idx(out_rotat_idx<=size(mv_data,1));
out_rotat_idx_r=out_rotat_idx(out_rotat_idx>size(mv_data,1)&out_rotat_idx<=2*size(mv_data,1))-size(mv_data,1);
out_rotat_idx_y=out_rotat_idx(out_rotat_idx>2*size(mv_data,1))-2*size(mv_data,1);
%update text
set(handles.rtthresh,'String',num2str(rotat_thresh));
setappdata(handles.rtthresh,'rt_p_outliers',out_rotat_idx_p);
setappdata(handles.rtthresh,'rt_r_outliers',out_rotat_idx_r);
setappdata(handles.rtthresh,'rt_y_outliers',out_rotat_idx_y);
%if composite measure
if (get(handles.norms,'Value') == get(handles.norms,'Max'))
setappdata(handles.rtthresh,'rt_norm_outliers',[]);
legend(handles.rotatGraph,'off');
cla(handles.rotatGraph);
set([handles.rotatGraph,handles.rt_up,handles.rt_down,handles.rtthresh,handles.text13],'visible','off')
set(handles.axes_mask,'position',get(handles.axes_mask,'position').*[1,1,1,0]+[0,0,0,.40]);
set(handles.all_outliers,'position',get(handles.all_outliers,'position').*[1,1,1,0]+[0,0,0,.29]);
return
end
set(handles.axes_mask,'position',get(handles.axes_mask,'position').*[1,1,1,0]+[0,0,0,.25]);
set(handles.all_outliers,'position',get(handles.all_outliers,'position').*[1,1,1,0]+[0,0,0,.14]);
axes(handles.rotatGraph);
cla;
plot(mv_data(:,4:6));
set(gca,'xlim',[0,size(mv_data,1)+1]);
l=ylabel('rotation \newline [rad]');
set(l,'VerticalAlignment','Bottom');
set(gca,'XTickLabel',[]);
legend('pitch', 'roll', 'yaw', 'Location', 'East');
h = gca;
set(h,'Ygrid','on');
thresh_rt_x = 1:length(mv_data);
thresh_rt_y = rotat_thresh*ones(1,length(mv_data));
y_lim = get(gca, 'YLim');
line(thresh_rt_x, thresh_rt_y, 'Color', 'black');
line(thresh_rt_x, -1*thresh_rt_y, 'Color', 'black');
for i = 1:length(out_rotat_idx_p)
line((out_rotat_idx_p(i)*ones(1, 2)), y_lim, 'Color', 'black');
end
for i = 1:length(out_rotat_idx_r)
line((out_rotat_idx_r(i)*ones(1, 2)), y_lim, 'Color', 'black');
end
for i = 1:length(out_rotat_idx_y)
line((out_rotat_idx_y(i)*ones(1, 2)), y_lim, 'Color', 'black');
end
set([h,handles.rt_up,handles.rt_down,handles.rtthresh,handles.text13],'visible','on')
end
% -----------------------------------------------------------------------
% UPDATESUMMARYPLOT
% This function generates the summary plot (gui 1st plot from the top)
% displaying the global signal and all of the identified outliers, and
% optionally the design matrix information and corresponding breakdown of
% outliers by condition
% -----------------------------------------------------------------------
function UpdateSummaryplot(hObject, eventdata, handles)
g = getappdata(handles.zthresh,'g');
num_sess = length(g);
tmps = get(handles.all_outliers,'String');
if ~isempty(tmps)
nstrings = size(tmps,1);
idx=cell(1,nstrings);
for i=1:nstrings
idx{i}=round(str2num(tmps(i,:)));
end
idx=cat(2,idx{:});
set(handles.all_outliers, 'String', int2str(idx));
else
idx = [];
end
%plot global mean
axes(handles.globalMean);
cla;
% BEGIN ohinds 2008-04-23: plot and print global mean
hold on;
cur_sess_start=1;
rng_mean=0;rng_minmax=[-inf,-inf];
for sess=1:num_sess
%rng{sess} = range(g{sess}(:,1));
rng_mean=rng_mean+mean(g{sess}(:,1));
rng_minmax=max(rng_minmax,[-min(g{sess}(:,1)),max(g{sess}(:,1))]);
plot(cur_sess_start:cur_sess_start+length(g{sess}(:,1))-1, g{sess}(:,1));
%ylabstr = sprintf('%s %f (%d)', ylabstr, rng{sess}, sess); % ohinds: can't put the range for all sessions on the ylabel, not enough room
cur_sess_start = cur_sess_start + length(g{sess}(:,1));
end
set(gca,'xlim',[0,cur_sess_start]);%,'ylim',rng_minmax.*[-1,1]);
rng_mean=rng_mean/num_sess;
ylabel('mean image\newlineintensity');
xlabel('scans');
% END ohinds 2008-04-23: plot global mean
y_lim = get(gca, 'YLim');
cur_sess_start=1;
for sess=1:num_sess
patch(cur_sess_start+[0,0,(length(g{sess}(:,1))-1)*[1,1]],[ylim,fliplr(ylim)],-ones(1,4),.9+.05*rem(sess,2)*[1,1,1],'edgecolor','none');
if num_sess>1
text(cur_sess_start+(length(g{sess}(:,1))-1)/2,ylim*[-.1;1.1],['Session ',num2str(sess)],'horizontalalignment','center');
end
cur_sess_start = cur_sess_start + length(g{sess}(:,1));
end
for i = 1:length(idx)
line((idx(i)*ones(1, 2)), y_lim, 'Color', 'red');
end
hold off;
analyses=getappdata(handles.savefile,'analyses');
analyses.outliers.scans=idx;
setappdata(handles.savefile,'analyses',analyses);
%show design (moved to global plot <alfnie> 2009-01)
if (get(handles.showDesign,'Value') == get(handles.showDesign,'Max'))
[SPM,design,names] = get_design(handles);
stats_file=getappdata(handles.showDesign,'SPMfile');
hold on
colors = {'k:','b:','r:','g:','c:','m:','y:'};
h=plot(1,nan,'.','markersize',1);
for i=1:size(design,2)
h(i+1)=plot(1:size(design,1) , rng_mean+sum(rng_minmax)/2*design(:,i),colors{mod(i,5)+1},'MarkerSize',4);
end
% computes number of outliers per condition <alfnie> 2009-01
out_idx=round(idx(idx>0));
if cur_sess_start-1~=size(design,1),
disp(['warning: incorrect number of scans (design matrix: ',num2str(size(design,1)),' ; functional data: ',num2str(cur_sess_start-1),')']);
outliers_per_condition=length(out_idx);
else
outliers_per_condition=[length(out_idx),sum(abs(design(out_idx,:)>0),1); size(design,1),sum(abs(design(:,:)>0),1)];
end
if size(outliers_per_condition,2)==length(names)+1,
legendnames={['Total :',num2str(outliers_per_condition(1,1)),' outlier scans (',num2str(100*outliers_per_condition(1,1)/max(eps,outliers_per_condition(2,1)),'%0.0f'),'%)']};
for i=1:length(names), legendnames{i+1}=[names{i},' :',num2str(outliers_per_condition(1,i+1),'%0.0f'),' outlier scans (',num2str(100*outliers_per_condition(1,i+1)/max(eps,outliers_per_condition(2,i+1)),'%0.0f'),'%)']; end
legend(h,legendnames{:});
end
analyses=getappdata(handles.savefile,'analyses');
analyses.outliers.condition_effects=outliers_per_condition(1,2:end);
analyses.outliers.condition_names=names;
setappdata(handles.savefile,'analyses',analyses);
[statsfile_path,statsfile_name] = fileparts(stats_file); if isempty(statsfile_path),statsfile_path=pwd;end;
stats_file_outliers=fullfile(statsfile_path,[statsfile_name,'_outliers.txt']);
if ~isempty(stats_file_outliers)
for fpidx = 1:2,
if fpidx==1,
fp=1;
fprintf('Number of outliers\n');
else
fp=fopen(stats_file_outliers,'w');
fprintf('saving outlier statistics to %s\n',stats_file_outliers);
end
if fp ~= -1
fprintf(fp,'%10s ','Total');
fprintf(fp,'%10s ',names{:});
fprintf(fp,'\n');
fprintf(fp,'%10.0f ',outliers_per_condition(1,:));
fprintf(fp,'\n');