-
Notifications
You must be signed in to change notification settings - Fork 198
Expand file tree
/
Copy pathphase_encoding.cpp
More file actions
388 lines (322 loc) · 14.7 KB
/
phase_encoding.cpp
File metadata and controls
388 lines (322 loc) · 14.7 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
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
/* Copyright (c) 2008-2025 the MRtrix3 contributors.
*
* This Source Code Form is subject to the terms of the Mozilla Public
* License, v. 2.0. If a copy of the MPL was not distributed with this
* file, You can obtain one at http://mozilla.org/MPL/2.0/.
*
* Covered Software is provided under this License on an "as is"
* basis, without warranty of any kind, either expressed, implied, or
* statutory, including, without limitation, warranties that the
* Covered Software is free of defects, merchantable, fit for a
* particular purpose or non-infringing.
* See the Mozilla Public License v. 2.0 for more details.
*
* For more details, see http://www.mrtrix.org/.
*/
#include "metadata/phase_encoding.h"
#include "exception.h"
namespace MR {
namespace Metadata {
namespace PhaseEncoding {
using namespace App;
const OptionGroup ImportOptions =
OptionGroup("Options for importing phase-encode tables")
+ Option("import_pe_table", "import a phase-encoding table from file")
+ Argument("file").type_file_in()
+ Option("import_pe_eddy", "import phase-encoding information from an EDDY-style config / index file pair")
+ Argument("config").type_file_in()
+ Argument("indices").type_file_in();
const OptionGroup SelectOptions =
OptionGroup("Options for selecting volumes based on phase-encoding")
+ Option("pe",
"select volumes with a particular phase encoding;"
" this can be three comma-separated values"
" (for i,j,k components of vector direction)"
" or four (direction & total readout time)")
+ Argument("desc").type_sequence_float();
const OptionGroup ExportOptions =
OptionGroup("Options for exporting phase-encode tables")
+ Option("export_pe_table", "export phase-encoding table to file")
+ Argument("file").type_file_out()
+ Option("export_pe_eddy", "export phase-encoding information to an EDDY-style config / index file pair")
+ Argument("config").type_file_out()
+ Argument("indices").type_file_out();
void check(const scheme_type& PE) {
if (PE.rows() == 0)
throw Exception("No valid phase encoding table found");
if (PE.cols() < 3)
throw Exception("Phase-encoding matrix must have at least 3 columns");
for (ssize_t row = 0; row != PE.rows(); ++row) {
for (ssize_t axis = 0; axis != 3; ++axis) {
if (std::round(PE(row, axis)) != PE(row, axis))
throw Exception("Phase-encoding matrix contains non-integral axis designation");
}
}
}
void check(const scheme_type& PE, const Header& header) {
check(PE);
const ssize_t num_volumes = (header.ndim() > 3) ? header.size(3) : 1;
if (num_volumes != PE.rows())
throw Exception("Number of volumes in image \"" + header.name() + "\" (" + str(num_volumes) + ")"
+ " does not match that in phase encoding table (" + str(PE.rows()) + ")");
}
namespace {
void erase(KeyValues& keyval, const std::string& s) {
auto it = keyval.find(s);
if (it != keyval.end())
keyval.erase(it);
};
}
void set_scheme(KeyValues& keyval, const scheme_type& PE) {
if (PE.rows() == 0) {
erase(keyval, "pe_scheme");
erase(keyval, "PhaseEncodingDirection");
erase(keyval, "TotalReadoutTime");
return;
}
std::string pe_scheme;
std::string first_line;
bool variation = false;
for (ssize_t row = 0; row < PE.rows(); ++row) {
std::string line = str(PE(row, 0));
for (ssize_t col = 1; col < PE.cols(); ++col)
line += "," + str(PE(row, col), 3);
add_line(pe_scheme, line);
if (first_line.empty())
first_line = line;
else if (line != first_line)
variation = true;
}
if (variation) {
keyval["pe_scheme"] = pe_scheme;
erase(keyval, "PhaseEncodingDirection");
erase(keyval, "TotalReadoutTime");
} else {
erase(keyval, "pe_scheme");
const Metadata::BIDS::axis_vector_type dir{int(PE(0, 0)), int(PE(0, 1)), int(PE(0, 2))};
keyval["PhaseEncodingDirection"] = Metadata::BIDS::vector2axisid(dir);
if (PE.cols() >= 4)
keyval["TotalReadoutTime"] = str(PE(0, 3), 3);
else
erase(keyval, "TotalReadoutTime");
}
}
void clear_scheme(KeyValues& keyval) {
auto erase = [&](const std::string& s) {
auto it = keyval.find(s);
if (it != keyval.end())
keyval.erase(it);
};
erase("pe_scheme");
erase("PhaseEncodingDirection");
erase("TotalReadoutTime");
}
scheme_type parse_scheme(const KeyValues& keyval, const Header& header) {
scheme_type PE;
const auto it = keyval.find("pe_scheme");
if (it != keyval.end()) {
try {
PE = MR::parse_matrix(it->second);
} catch (Exception& e) {
throw Exception(e, "malformed PE scheme associated with image \"" + header.name() + "\"");
}
if (ssize_t(PE.rows()) != ((header.ndim() > 3) ? header.size(3) : 1))
throw Exception("malformed PE scheme associated with image \"" + header.name() + "\": "
+ "number of rows does not equal number of volumes");
} else {
const auto it_dir = keyval.find("PhaseEncodingDirection");
if (it_dir != keyval.end()) {
const auto it_time = keyval.find("TotalReadoutTime");
const size_t cols = it_time == keyval.end() ? 3 : 4;
Eigen::Matrix<default_type, Eigen::Dynamic, 1> row(cols);
try {
row.head(3) = BIDS::axisid2vector(it_dir->second).cast<default_type>();
} catch (Exception& e) {
throw Exception(e, "malformed phase encoding direction associated with image \"" + header.name() + "\"");
}
if (it_time != keyval.end()) {
try {
row[3] = to<default_type>(it_time->second);
} catch (Exception& e) {
throw Exception(e, "Error adding readout time to phase encoding table");
}
}
PE.resize((header.ndim() > 3) ? header.size(3) : 1, cols);
PE.rowwise() = row.transpose();
}
}
return PE;
}
scheme_type get_scheme(const Header& header) {
DEBUG("searching for suitable phase encoding data...");
using namespace App;
scheme_type result;
try {
const auto opt_table = get_options("import_pe_table");
if (!opt_table.empty())
result = load(opt_table[0][0], header);
const auto opt_eddy = get_options("import_pe_eddy");
if (!opt_eddy.empty()) {
if (!opt_table.empty())
throw Exception("Phase encoding table can be provided"
" using either -import_pe_table or -import_pe_eddy option,"
" but NOT both");
result = load_eddy(opt_eddy[0][0], opt_eddy[0][1], header);
}
if (opt_table.empty() && opt_eddy.empty())
result = parse_scheme(header.keyval(), header);
} catch (Exception &e) {
throw Exception(e, "error importing phase encoding table for image \"" + header.name() + "\"");
}
if (result.rows() == 0)
return result;
if (result.cols() < 3)
throw Exception("unexpected phase encoding table matrix dimensions");
INFO("found " + str(result.rows()) + "x" + str(result.cols()) + " phase encoding table");
return result;
}
void transform_for_image_load(KeyValues& keyval, const Header& H) {
scheme_type pe_scheme;
try {
pe_scheme = parse_scheme(keyval, H);
} catch (Exception& e) {
WARN("Unable to conform phase encoding information to image realignment "
" for image \"" + H.name() + "\"; erasing");
clear_scheme(keyval);
return;
}
if (pe_scheme.rows() == 0) {
DEBUG("No phase encoding information found for transformation with load of image \"" + H.name() + "\"");
return;
}
if (H.realignment().is_identity()) {
INFO("No transformation of phase encoding data for load of image \"" + H.name() + "\" required");
return;
}
set_scheme(keyval, transform_for_image_load(pe_scheme, H));
INFO("Phase encoding data transformed to match RAS realignment of image \"" + H.name() + "\"");
}
scheme_type transform_for_image_load(const scheme_type& pe_scheme, const Header& H) {
if (H.realignment().is_identity())
return pe_scheme;
scheme_type result(pe_scheme.rows(), pe_scheme.cols());
for (ssize_t row = 0; row != pe_scheme.rows(); ++row) {
Eigen::VectorXd new_line = pe_scheme.row(row);
new_line.head<3>() = (H.realignment().applied_transform() * new_line.head<3>().cast<int>()).cast<default_type>();
result.row(row) = new_line;
}
return result;
}
void transform_for_nifti_write(KeyValues& keyval, const Header& H) {
const scheme_type pe_scheme = parse_scheme(keyval, H);
if (pe_scheme.rows() == 0) {
DEBUG("No phase encoding information found for transformation with save of NIfTI image \"" + H.name() + "\"");
return;
}
set_scheme(keyval, transform_for_nifti_write(pe_scheme, H));
}
scheme_type transform_for_nifti_write(const scheme_type& pe_scheme, const Header& H) {
if (pe_scheme.rows() == 0)
return pe_scheme;
Axes::Shuffle shuffle = File::NIfTI::axes_on_write(H);
if (shuffle.is_identity()) {
INFO("No transformation of phase encoding data required for export to file:"
" output image will be RAS");
return pe_scheme;
}
scheme_type result(pe_scheme.rows(), pe_scheme.cols());
for (ssize_t row = 0; row != pe_scheme.rows(); ++row) {
Eigen::VectorXd new_line = pe_scheme.row(row);
for (ssize_t axis = 0; axis != 3; ++axis)
new_line[axis] =
pe_scheme(row, shuffle.permutations[axis]) != 0.0 && shuffle.flips[axis] ?
-pe_scheme(row, shuffle.permutations[axis]) :
pe_scheme(row, shuffle.permutations[axis]);
result.row(row) = new_line;
}
INFO("Phase encoding data transformed to match NIfTI / MGH image export prior to writing to file");
return result;
}
void scheme2eddy(const scheme_type& PE, Eigen::MatrixXd& config, Eigen::Array<int, Eigen::Dynamic, 1>& indices) {
try {
check(PE);
} catch (Exception& e) {
throw Exception(e, "Cannot convert phase-encoding scheme to eddy format");
}
if (PE.cols() != 4)
throw Exception("Phase-encoding matrix requires 4 columns to convert to eddy format");
config.resize(0, 0);
indices = Eigen::Array<int, Eigen::Dynamic, 1>::Constant (PE.rows(), PE.rows());
for (ssize_t PE_row = 0; PE_row != PE.rows(); ++PE_row) {
for (ssize_t config_row = 0; config_row != config.rows(); ++config_row) {
const bool dir_match = PE.template block<1, 3>(PE_row, 0).isApprox(config.block<1, 3>(config_row, 0));
const bool time_match = abs(PE(PE_row, 3) - config(config_row, 3)) < 1e-3;
if (dir_match && time_match) {
// FSL-style index file indexes from 1
indices[PE_row] = config_row + 1;
break;
}
}
if (indices[PE_row] == PE.rows()) {
// No corresponding match found in config matrix; create a new entry
config.conservativeResize(config.rows() + 1, 4);
config.row(config.rows() - 1) = PE.row(PE_row);
indices[PE_row] = config.rows();
}
}
}
scheme_type eddy2scheme(const Eigen::MatrixXd& config, const Eigen::Array<int, Eigen::Dynamic, 1>& indices) {
if (config.cols() != 4)
throw Exception("Expected 4 columns in EDDY-format phase-encoding config file");
scheme_type result(indices.size(), 4);
for (ssize_t row = 0; row != indices.size(); ++row) {
if (indices[row] > config.rows())
throw Exception("Malformed EDDY-style phase-encoding information:"
" index exceeds number of config entries");
result.row(row) = config.row(indices[row] - 1);
}
return result;
}
void export_commandline(const Header& header) {
auto check = [&](const scheme_type& m) -> const scheme_type & {
if (m.rows() == 0)
throw Exception("no phase-encoding information found within image \"" + header.name() + "\"");
return m;
};
auto scheme = parse_scheme(header.keyval(), header);
auto opt = get_options("export_pe_table");
if (!opt.empty())
save(check(scheme), header, opt[0][0]);
opt = get_options("export_pe_eddy");
if (!opt.empty())
save_eddy(check(scheme), header, opt[0][0], opt[0][1]);
}
scheme_type load(const std::string& path, const Header& header) {
const scheme_type PE = load_matrix(path);
check(PE, header);
// As with JSON import, need to query the header to discover if the
// strides / transform were modified on image load to make the image
// data appear approximately axial, in which case we need to apply the
// same transforms to the phase encoding data on load
return transform_for_image_load(PE, header);
}
scheme_type load_eddy(const std::string& config_path, const std::string& index_path, const Header& header) {
const Eigen::MatrixXd config = load_matrix(config_path);
const Eigen::Array<int, Eigen::Dynamic, 1> indices = load_vector<int>(index_path);
const scheme_type PE = eddy2scheme(config, indices);
check(PE, header);
return transform_for_image_load(PE, header);
}
void save(const scheme_type& PE, const std::string& path) {
File::OFStream out(path);
for (ssize_t row = 0; row != PE.rows(); ++row) {
// Write phase-encode direction as integers; other information as floating-point
out << PE.template block<1, 3>(row, 0).template cast<int>();
if (PE.cols() > 3)
out << " " << PE.block(row, 3, 1, PE.cols() - 3);
out << "\n";
}
}
}
}
}