-
Notifications
You must be signed in to change notification settings - Fork 2
/
art_threshdown.m
113 lines (94 loc) · 3.36 KB
/
art_threshdown.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
function art_threshdown
% Function art_threshdown (v.2.2)
%
% Helper function to art_global. Grabs the saved global data
% from the main figure and extracts the relevant current
% threshold info from the displayed current threshold, then
% decreases the threshold by a small increment for both movement
% and global intensity, re-calculates the outlier list, and
% redraws the text list and graphical outlier display.
% ------------------------------
% Jeff Cooper, Nov. 2002
% Updated for motion. Paul Mazaika, Aug. 2006
% fix plot v2.2
% v3 pkm. Adjust global intensity thresh by 0.05%. mar 2009
% -------------------
% Defaults increment
% -------------------
incr = 0.1; % For global intensity in z-score (obsolete)
incrpct = 0.05; % For global intensity in percent
incrmv = 0.05; % For scan to scan motion
% -------------------
% Interrogate GUI
% -------------------
handles = guihandles;
g = guidata(gcf);
curr_thresh = str2num(get(handles.threshnum, 'String'));
curr_threshpct = str2num(get(handles.threshnumpct, 'String'));
curr_threshmv = str2num(get(handles.threshnummv, 'String'));
pctmap = curr_threshpct/curr_thresh;
incr = incrpct/pctmap;
% -------------------
% Main action
% -------------------
%keyboard;
delta = squeeze(g(:,2));
g = squeeze(g(:,1));
new_thresh = curr_thresh - incr;
zscoreA = (g - mean(g))./std(g); % in case Matlab zscore is not available
glout_idx = (find(abs(zscoreA) > new_thresh))';
new_threshmv = curr_threshmv - incrmv;
mvout_idx = (find(delta > new_threshmv))';
% All volumes to be repaired from global intensity and motion
out_idx = unique([glout_idx mvout_idx]);
outdw_idx = out_idx;
% Update text
set(handles.threshnum, 'String', num2str(new_thresh));
set(handles.threshnumpct, 'String', num2str(new_thresh*pctmap));
set(handles.threshnummv, 'String', num2str(new_threshmv));
% Update edit text box
set(handles.indexedit, 'String', int2str(out_idx));
set(handles.deweightlist, 'String', int2str(outdw_idx));
% Update graphics for global intensity
subplot(5,1,2);
cla;
zscoreA = (g - mean(g))./std(g); % in case Matlab zscore is not available
plot(abs(zscoreA));
ylabel('Std away from mean');
thresh_x = 1:length(g);
thresh_y = new_thresh*ones(1,length(g));
line(thresh_x, thresh_y, 'Color', 'r');
axes_lim = get(gca, 'YLim');
axes_height = [axes_lim(1) axes_lim(2)];
for i = 1:length(glout_idx)
line((glout_idx(i)*ones(1, 2)), axes_height, 'Color', 'r');
end
% Update graphics for movement
subplot(5,1,4);
cla;
plot(delta);
ylabel('Motion (mm/TR)');
xlabel('Scan to scan movement (~mm). Rotation assumes 65 mm from origin');
thresh_x = 1:length(g);
thresh_y = new_threshmv*ones(1,length(g));
line(thresh_x, thresh_y, 'Color', 'r');
axes_lim = get(gca, 'YLim');
axes_height = [ axes_lim(1) axes_lim(2)];
for i = 1:length(mvout_idx)
line((mvout_idx(i)*ones(1, 2)), axes_height, 'Color', 'r');
end
% Update top chart for deweighting indices
subplot(5,1,1)
axes_lim = get(gca, 'YLim');
axes_height = [axes_lim(1) axes_lim(2)];
outerase = [1:length(delta)];
for i = 1:length(outerase)
line((outerase(i)*ones(1, 2)), axes_height, 'Color', 'w');
end
plot(g,'b');
ylabel('Global Avg. Signal');
xlabel('Red vertical lines are to repair. Green vertical lines are to deweight.');
title('ArtifactRepair GUI to repair outliers and identify scans to deweight')
for i = 1:length(outdw_idx)
line((outdw_idx(i)*ones(1, 2)), axes_height, 'Color', 'r');
end