-
Notifications
You must be signed in to change notification settings - Fork 13
/
RobustSpectralRegistration.m
254 lines (211 loc) · 9.5 KB
/
RobustSpectralRegistration.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
function [AllFramesFTrealign, MRS_struct] = RobustSpectralRegistration(MRS_struct)
% Align using robust spectral registration (Mikkelsen et al. NMR Biomed.
% 2020;33(10):e4368. doi:10.1002/nbm.4368)
ii = MRS_struct.ii;
% Looping parameters
if MRS_struct.p.HERMES
SpecRegLoop = 3;
SubspecToAlign = repmat([3 2 1 0], [1 size(MRS_struct.fids.data,2)/4]);
else
SpecRegLoop = 1;
SubspecToAlign = MRS_struct.fids.ON_OFF;
end
% Pre-allocate memory
if MRS_struct.p.HERMES
params = zeros(size(MRS_struct.fids.data,2)/4,2);
MSE = zeros(1,size(MRS_struct.fids.data,2)/4);
w = zeros(1,size(MRS_struct.fids.data,2)/4);
else
params = zeros(size(MRS_struct.fids.data,2)/2,2);
MSE = zeros(1,size(MRS_struct.fids.data,2)/2);
w = zeros(1,size(MRS_struct.fids.data,2)/2);
end
MRS_struct.out.SpecReg.freq{ii} = zeros(1,size(MRS_struct.fids.data,2));
MRS_struct.out.SpecReg.phase{ii} = zeros(1,size(MRS_struct.fids.data,2));
MRS_struct.out.SpecReg.MSE{ii} = zeros(1,size(MRS_struct.fids.data,2));
zMSE = zeros(1,size(MRS_struct.fids.data,2));
MRS_struct.fids.data_align = complex(zeros(size(MRS_struct.fids.data)));
DataToAlign = complex(zeros(size(MRS_struct.fids.data)));
MSEfun = @(a,b) sum((a - b).^2) / length(a);
% Optimization options
lsqnonlinopts = optimoptions(@lsqnonlin);
lsqnonlinopts = optimoptions(lsqnonlinopts,'Algorithm','levenberg-marquardt','Display','off');
% Automatic unstable lipid/residual water removal
freqRange = MRS_struct.p.sw(ii) / MRS_struct.p.LarmorFreq(ii);
freq = (MRS_struct.p.npoints(ii) + 1 - (1:MRS_struct.p.npoints(ii))) / MRS_struct.p.npoints(ii) * freqRange + 4.68 - freqRange/2;
waterLim = freq <= 4.68+0.25 & freq >= 4.68-0.25;
lipidLim = freq <= 1.85 & freq >= 0;
noiseLim = freq <= 9 & freq >= 8;
spec = fftshift(fft(MRS_struct.fids.data,[],1),1);
S = mean(real(spec),2);
noise = S(noiseLim);
fit = polyval(polyfit(freq(noiseLim), noise.', 2), freq(noiseLim)); % detrend noise
noise = noise - fit.';
r = std(S(lipidLim)) / std(noise);
if MRS_struct.p.HERMES
ind = all(MRS_struct.fids.ON_OFF' == 0,2);
ind2 = ind;
else
switch MRS_struct.p.target{1}
case {'GABAGlx','GABA','Glx','Lac','EtOH'}
ind = 1:size(MRS_struct.fids.data,2);
ind2 = MRS_struct.fids.ON_OFF == 0;
case 'GSH'
ind = MRS_struct.fids.ON_OFF == 0;
ind2 = ind;
end
end
q = sum(abs(spec(waterLim,ind))) * abs(freq(1) - freq(2));
q = (q - median(q)) / median(q) * 100;
q = sum(abs(q) > 40) / length(q);
% Force-run water filter if very strong water suppression was used
S = mean(abs(spec(:,ind2)),2);
maxNAA = max(S(freq <= 4.25 & freq >= 1.8));
maxWater = max(S(freq <= 4.68+0.22 & freq >= 4.68-0.22));
if maxWater / maxNAA < 1.5
q = 1;
end
r_threshold = 40;
q_threshold = 0.1;
lipid_flag = 0;
water_flag = 0;
if r > r_threshold || q > q_threshold
if r > r_threshold
lipid_flag = 1;
end
if q > q_threshold
water_flag = 1;
end
% Turn off warnings about the legacy random number generator if they are currently on
w1 = warning('query','MATLAB:RandStream:ActivatingLegacyGenerators');
w2 = warning('query','MATLAB:RandStream:ReadingInactiveLegacyGeneratorState');
if strcmp(w1.state,'on')
warning('off','MATLAB:RandStream:ActivatingLegacyGenerators');
end
if strcmp(w2.state,'on')
warning('off','MATLAB:RandStream:ReadingInactiveLegacyGeneratorState');
end
reverseStr = '';
for jj = 1:size(MRS_struct.fids.data,2)
if lipid_flag && ~water_flag
msg = sprintf('\nUnstable lipid contamination detected. Temporarily removing lipids from transient: %d', jj);
elseif ~lipid_flag && water_flag
msg = sprintf('\nUnstable residual water detected. Temporarily removing residual water from transient: %d', jj);
elseif lipid_flag && water_flag
msg = sprintf('\nUnstable lipid contamination and residual water detected.\nTemporarily removing lipids and residual water from transient: %d', jj);
end
fprintf([reverseStr, msg]);
reverseStr = repmat(sprintf('\b'), 1, length(msg));
DataToAlign(:,jj) = SignalFilter(spec(:,jj), lipid_flag, water_flag, jj, MRS_struct);
end
% Turn warnings back on if they were previously on
if strcmp(w1.state,'on')
warning('on','MATLAB:RandStream:ActivatingLegacyGenerators');
end
if strcmp(w2.state,'on')
warning('on','MATLAB:RandStream:ReadingInactiveLegacyGeneratorState');
end
else
DataToAlign = MRS_struct.fids.data;
end
time = (0:(MRS_struct.p.npoints(ii)-1))'/MRS_struct.p.sw(ii);
% Spectral registration
count = 1;
reverseStr = '';
while SpecRegLoop > -1
% Use first n points of time-domain data, where n is the last point where abs(diff(mean(SNR))) > 0.5
signal = abs(DataToAlign(:,SubspecToAlign == SpecRegLoop));
noise = 2*std(signal(ceil(0.75*size(signal,1)):end,:));
SNR = signal ./ repmat(noise, [size(DataToAlign,1) 1]);
SNR = abs(diff(mean(SNR,2)));
SNR = SNR(time <= 0.2); % use no more than 200 ms of data
tMax = find(SNR > 0.5,1,'last');
if isempty(tMax) || tMax < find(time <= 0.05,1,'last') % use at least 50 ms of data
% (shortened this from 100 ms because it seems
% like this helps when there are spurious echoes
% or strong water suppression was used)
tMax = find(time <= 0.05,1,'last');
end
% Flatten complex data for use in spectral registration
clear flatdata
flatdata(:,1,:) = real(DataToAlign(1:tMax,SubspecToAlign == SpecRegLoop));
flatdata(:,2,:) = imag(DataToAlign(1:tMax,SubspecToAlign == SpecRegLoop));
% Determine optimal alignment order by calculating a similarity metric (mean squared error)
if size(MRS_struct.fids.data,2) <= 4
alignOrder = 1;
else
D = zeros(size(flatdata,3));
ind = find(SubspecToAlign == SpecRegLoop);
for jj = 1:size(flatdata,3)
for kk = 1:size(flatdata,3)
D(jj,kk) = MSEfun(real(DataToAlign(1:tMax,ind(jj))), real(DataToAlign(1:tMax,ind(kk))));
end
end
D(~D) = NaN;
d = median(D,'omitnan');
[~,alignOrder] = sort(d);
end
% Set initial reference transient based on similarity index
target = squeeze(flatdata(:,:,alignOrder(1)));
target = target(:);
% Scalar to normalize transients (reduces optimization time)
a = max(abs(target));
% Pre-allocate memory
m = zeros(length(target), size(flatdata,3));
% Starting values for optimization
f0 = MRS_struct.spec.F0freq2{ii}(ind) * MRS_struct.p.LarmorFreq(ii);
f0 = f0(alignOrder);
f0 = f0 - f0(1);
phi0 = zeros(size(f0));
x0 = [f0(:) phi0(:)];
% Determine frequency and phase offsets by spectral registration
t = 0:(1/MRS_struct.p.sw(ii)):(length(target)/2-1)*(1/MRS_struct.p.sw(ii));
iter = 1;
for jj = alignOrder
msg = sprintf('\nRunning robust spectral registration on transient: %d', count);
fprintf([reverseStr, msg]);
reverseStr = repmat(sprintf('\b'), 1, length(msg));
count = count + 1;
transient = squeeze(flatdata(:,:,jj));
fun = @(x) SpecReg(transient(:)/a, target/a, t, x);
params(jj,:) = lsqnonlin(fun, x0(iter,:), [], [], lsqnonlinopts);
f = params(jj,1);
phi = params(jj,2);
m_c = complex(flatdata(:,1,jj), flatdata(:,2,jj));
m_c = m_c .* exp(1i*pi*(t'*f*2+phi/180));
m(:,jj) = [real(m_c); imag(m_c)];
resid = target - m(:,jj);
MSE(jj) = sum(resid.^2) / (length(resid) - 2);
% Update reference
w(jj) = 0.5*corr(target, m(:,jj)).^2;
target = (1 - w(jj))*target + w(jj)*m(:,jj);
iter = iter + 1;
end
ind = find(SubspecToAlign == SpecRegLoop);
MRS_struct.out.SpecReg.freq{ii}(ind) = params(:,1);
MRS_struct.out.SpecReg.phase{ii}(ind) = params(:,2);
MRS_struct.out.SpecReg.MSE{ii}(ind) = MSE;
zMSE(ind(MSE > 0)) = zscore(MSE(MSE > 0)); % standardized MSEs (exclude the first reference)
% Apply frequency and phase corrections to raw data
MRS_struct.fids.data_align(:,ind) = MRS_struct.fids.data(:,ind) ...
.* exp(repmat(1i * params(:,1)' * 2 * pi, [length(time) 1]) .* repmat(time, [1 size(flatdata,3)])) ...
.* repmat(exp(1i * pi/180 * params(:,2)'), [length(time) 1]);
if SpecRegLoop == 0
% Align subspectra
if MRS_struct.p.use_prealign_ref
MRS_struct = AlignSubSpectra_PreAlignRef(MRS_struct, water_flag);
else
MRS_struct = AlignSubSpectra(MRS_struct, water_flag, r);
end
% Line-broadening, zero-filling and FFT
AllFramesFTrealign = MRS_struct.fids.data_align .* repmat((exp(-time * MRS_struct.p.LB * pi)), [1 size(MRS_struct.fids.data,2)]);
AllFramesFTrealign = fftshift(fft(AllFramesFTrealign, MRS_struct.p.ZeroFillTo(ii), 1),1);
% Reject transients that are greater than 3 st. devs. of zMSEs
% (only applies if not using weighted averaging)
MRS_struct.out.SpecReg.zMSE{ii} = zMSE;
MRS_struct.out.reject{ii} = zMSE > 3;
end
SpecRegLoop = SpecRegLoop - 1;
end
fprintf('\n');
end