v2.0.0
Loading...
Searching...
No Matches
fs_atlas_lookup.cpp
Go to the documentation of this file.
1//=============================================================================================================
12
13//=============================================================================================================
14// INCLUDES
15//=============================================================================================================
16
17#include "fs_atlas_lookup.h"
18
19//=============================================================================================================
20// EIGEN INCLUDES
21//=============================================================================================================
22
23#include <Eigen/LU>
24
25//=============================================================================================================
26// QT INCLUDES
27//=============================================================================================================
28
29#include <QFile>
30#include <QDataStream>
31#include <QDebug>
32#ifndef WASMBUILD
33#include <QProcess>
34#endif
35#include <QtEndian>
36
37//=============================================================================================================
38// USED NAMESPACES
39//=============================================================================================================
40
41using namespace FSLIB;
42using namespace Eigen;
43
44//=============================================================================================================
45// DEFINE MEMBER METHODS
46//=============================================================================================================
47
48FsAtlasLookup::FsAtlasLookup()
49: m_ras2vox(Matrix4f::Identity())
50{
51 initLookupTable();
52}
53
54//=============================================================================================================
55
56bool FsAtlasLookup::load(const QString& sParcellationPath)
57{
58 m_loaded = false;
59
60 // Support .mgz (gzip-compressed MGH) and .mgh (uncompressed)
61 // For .mgz, use QProcess to decompress, or read via gzip.
62 // Here we support both by trying QFile directly (works for .mgh)
63 // and using a pipe for .mgz.
64
65 QByteArray data;
66
67 if (sParcellationPath.endsWith(QStringLiteral(".mgz"), Qt::CaseInsensitive)) {
68#ifdef WASMBUILD
69 qWarning() << "[FsAtlasLookup::load] .mgz files not supported in WebAssembly build (QProcess unavailable):" << sParcellationPath;
70 return false;
71#else
72 // Decompress mgz via gzip
73 QProcess gzipProc;
74 gzipProc.start(QStringLiteral("gzip"), QStringList() << QStringLiteral("-dc") << sParcellationPath);
75 if (!gzipProc.waitForFinished(30000)) {
76 qWarning() << "[FsAtlasLookup::load] gzip decompression timed out for" << sParcellationPath;
77 return false;
78 }
79 if (gzipProc.exitCode() != 0) {
80 qWarning() << "[FsAtlasLookup::load] gzip decompression failed for" << sParcellationPath;
81 return false;
82 }
83 data = gzipProc.readAllStandardOutput();
84#endif
85 } else {
86 QFile file(sParcellationPath);
87 if (!file.open(QIODevice::ReadOnly)) {
88 qWarning() << "[FsAtlasLookup::load] Cannot open" << sParcellationPath;
89 return false;
90 }
91 data = file.readAll();
92 file.close();
93 }
94
95 if (data.size() < 284 + 64) { // Minimum header size
96 qWarning() << "[FsAtlasLookup::load] File too small:" << sParcellationPath;
97 return false;
98 }
99
100 // Parse MGH header (all big-endian)
101 // Offset 0: version (int32) — must be 1
102 // Offset 4: width (int32)
103 // Offset 8: height (int32)
104 // Offset 12: depth (int32)
105 // Offset 16: nframes (int32)
106 // Offset 20: type (int32) — 0=uchar, 1=int, 3=float, 4=short
107 // Offset 284: vox2ras matrix (12 floats, row-major 3x4, then last row implicit [0,0,0,1])
108
109 const char* raw = data.constData();
110
111 auto readInt32 = [&](int offset) -> qint32 {
112 return qFromBigEndian<qint32>(raw + offset);
113 };
114
115 qint32 version = readInt32(0);
116 if (version != 1) {
117 qWarning() << "[FsAtlasLookup::load] Unsupported MGH version:" << version;
118 return false;
119 }
120
121 m_dimX = readInt32(4);
122 m_dimY = readInt32(8);
123 m_dimZ = readInt32(12);
124 qint32 nframes = readInt32(16);
125 qint32 type = readInt32(20);
126
127 if (m_dimX <= 0 || m_dimY <= 0 || m_dimZ <= 0) {
128 qWarning() << "[FsAtlasLookup::load] Invalid dimensions:" << m_dimX << m_dimY << m_dimZ;
129 return false;
130 }
131
132 // Read vox2ras matrix at offset 284 (stored as 4x4 float32 big-endian starting with ras_good_flag at 28)
133 // Actually MGH format: offset 284 has the vox2ras 4x4 as 16 big-endian floats
134 // But first check ras_good_flag at offset 284 - 4*16 = no...
135 // MGH header layout:
136 // 0-3: version (int32)
137 // 4-7: width
138 // 8-11: height
139 // 12-15: depth
140 // 16-19: nframes
141 // 20-23: type
142 // 24-27: dof
143 // 28-29: ras_good_flag (short)
144 // 30-41: x_size, y_size, z_size (3 floats — voxel sizes)
145 // 42-77: x_ras, y_ras, z_ras (3x3 floats — direction cosines)
146 // 78-89: c_ras (3 floats — center RAS)
147 // The actual data starts at offset 284.
148
149 // Read ras_good_flag
150 qint16 rasGoodFlag = qFromBigEndian<qint16>(raw + 28);
151
152 Matrix4f vox2ras = Matrix4f::Identity();
153
154 if (rasGoodFlag > 0) {
155 // Read voxel sizes and direction cosines to construct vox2ras
156 // Properly read floats by memcpy + endian swap
157 auto readBEFloat = [&](int offset) -> float {
158 quint32 bits = qFromBigEndian<quint32>(raw + offset);
159 float val;
160 memcpy(&val, &bits, sizeof(float));
161 return val;
162 };
163
164 float xsize = readBEFloat(30);
165 float ysize = readBEFloat(34);
166 float zsize = readBEFloat(38);
167
168 // Direction cosines (column vectors)
169 Vector3f x_ras(readBEFloat(42), readBEFloat(46), readBEFloat(50));
170 Vector3f y_ras(readBEFloat(54), readBEFloat(58), readBEFloat(62));
171 Vector3f z_ras(readBEFloat(66), readBEFloat(70), readBEFloat(74));
172
173 // Center RAS
174 Vector3f c_ras(readBEFloat(78), readBEFloat(82), readBEFloat(86));
175
176 // Construct vox2ras from direction cosines, voxel sizes, and center
177 // vox2ras maps voxel (i,j,k) to RAS (x,y,z):
178 // M = [Mdc * diag(sizes), Pcrs_c] where Pcrs_c is computed from c_ras
179 Matrix3f Mdc;
180 Mdc.col(0) = x_ras;
181 Mdc.col(1) = y_ras;
182 Mdc.col(2) = z_ras;
183
184 Matrix3f MdcD = Mdc;
185 MdcD.col(0) *= xsize;
186 MdcD.col(1) *= ysize;
187 MdcD.col(2) *= zsize;
188
189 // Center voxel
190 Vector3f Pcrs_center(static_cast<float>(m_dimX) / 2.0f,
191 static_cast<float>(m_dimY) / 2.0f,
192 static_cast<float>(m_dimZ) / 2.0f);
193
194 Vector3f Pxyz_0 = c_ras - MdcD * Pcrs_center;
195
196 vox2ras.block<3, 3>(0, 0) = MdcD;
197 vox2ras.block<3, 1>(0, 3) = Pxyz_0;
198 }
199
200 // Compute ras2vox = inverse(vox2ras)
201 m_ras2vox = vox2ras.inverse();
202
203 // Read voxel data starting at offset 284
204 const int headerSize = 284;
205 qint64 nVoxels = static_cast<qint64>(m_dimX) * m_dimY * m_dimZ * nframes;
206 m_voxelData.resize(static_cast<int>(nVoxels));
207
208 const char* voxPtr = raw + headerSize;
209 qint64 dataSize = data.size() - headerSize;
210
211 switch (type) {
212 case 0: { // MRI_UCHAR
213 if (dataSize < nVoxels) {
214 qWarning() << "[FsAtlasLookup::load] Insufficient data for uchar volume";
215 return false;
216 }
217 for (qint64 i = 0; i < nVoxels; ++i)
218 m_voxelData[static_cast<int>(i)] = static_cast<int>(static_cast<unsigned char>(voxPtr[i]));
219 break;
220 }
221 case 1: { // MRI_INT
222 if (dataSize < nVoxels * 4) {
223 qWarning() << "[FsAtlasLookup::load] Insufficient data for int volume";
224 return false;
225 }
226 for (qint64 i = 0; i < nVoxels; ++i)
227 m_voxelData[static_cast<int>(i)] = qFromBigEndian<qint32>(voxPtr + i * 4);
228 break;
229 }
230 case 3: { // MRI_FLOAT
231 if (dataSize < nVoxels * 4) {
232 qWarning() << "[FsAtlasLookup::load] Insufficient data for float volume";
233 return false;
234 }
235 for (qint64 i = 0; i < nVoxels; ++i) {
236 quint32 bits = qFromBigEndian<quint32>(voxPtr + i * 4);
237 float val;
238 memcpy(&val, &bits, sizeof(float));
239 m_voxelData[static_cast<int>(i)] = static_cast<int>(std::round(val));
240 }
241 break;
242 }
243 case 4: { // MRI_SHORT
244 if (dataSize < nVoxels * 2) {
245 qWarning() << "[FsAtlasLookup::load] Insufficient data for short volume";
246 return false;
247 }
248 for (qint64 i = 0; i < nVoxels; ++i)
249 m_voxelData[static_cast<int>(i)] = qFromBigEndian<qint16>(voxPtr + i * 2);
250 break;
251 }
252 default:
253 qWarning() << "[FsAtlasLookup::load] Unsupported MGH data type:" << type;
254 return false;
255 }
256
257 m_loaded = true;
258 return true;
259}
260
261//=============================================================================================================
262
263QString FsAtlasLookup::labelAtRas(const Vector3f& ras) const
264{
265 if (!m_loaded)
266 return QStringLiteral("Unknown");
267
268 Vector3i vox = rasToVoxel(ras);
269
270 // Bounds check
271 if (vox(0) < 0 || vox(0) >= m_dimX ||
272 vox(1) < 0 || vox(1) >= m_dimY ||
273 vox(2) < 0 || vox(2) >= m_dimZ)
274 return QStringLiteral("Unknown");
275
276 // MGH stores data in column-major Fortran order: index = x + dimX * (y + dimY * z)
277 int idx = vox(0) + m_dimX * (vox(1) + m_dimY * vox(2));
278 int label = m_voxelData.at(idx);
279
280 return m_lookupTable.value(label, QStringLiteral("Unknown"));
281}
282
283//=============================================================================================================
284
285QStringList FsAtlasLookup::labelsForPositions(const QVector<Vector3f>& positions) const
286{
287 QStringList result;
288 result.reserve(positions.size());
289 for (const auto& pos : positions)
290 result.append(labelAtRas(pos));
291 return result;
292}
293
294//=============================================================================================================
295
296bool FsAtlasLookup::isLoaded() const
297{
298 return m_loaded;
299}
300
301//=============================================================================================================
302
303Vector3i FsAtlasLookup::rasToVoxel(const Vector3f& ras) const
304{
305 Vector4f rasH(ras(0), ras(1), ras(2), 1.0f);
306 Vector4f voxH = m_ras2vox * rasH;
307 return Vector3i(static_cast<int>(std::round(voxH(0))),
308 static_cast<int>(std::round(voxH(1))),
309 static_cast<int>(std::round(voxH(2))));
310}
311
312//=============================================================================================================
313
314void FsAtlasLookup::initLookupTable()
315{
316 // Subcortical / standard FreeSurfer labels
317 m_lookupTable[0] = QStringLiteral("Unknown");
318 m_lookupTable[2] = QStringLiteral("Left-Cerebral-White-Matter");
319 m_lookupTable[3] = QStringLiteral("Left-Cerebral-Cortex");
320 m_lookupTable[4] = QStringLiteral("Left-Lateral-Ventricle");
321 m_lookupTable[5] = QStringLiteral("Left-Inf-Lat-Vent");
322 m_lookupTable[7] = QStringLiteral("Left-Cerebellum-White-Matter");
323 m_lookupTable[8] = QStringLiteral("Left-Cerebellum-Cortex");
324 m_lookupTable[10] = QStringLiteral("Left-Thalamus");
325 m_lookupTable[11] = QStringLiteral("Left-Caudate");
326 m_lookupTable[12] = QStringLiteral("Left-Putamen");
327 m_lookupTable[13] = QStringLiteral("Left-Pallidum");
328 m_lookupTable[14] = QStringLiteral("3rd-Ventricle");
329 m_lookupTable[15] = QStringLiteral("4th-Ventricle");
330 m_lookupTable[16] = QStringLiteral("Brain-Stem");
331 m_lookupTable[17] = QStringLiteral("Left-Hippocampus");
332 m_lookupTable[18] = QStringLiteral("Left-Amygdala");
333 m_lookupTable[24] = QStringLiteral("CSF");
334 m_lookupTable[26] = QStringLiteral("Left-Accumbens-area");
335 m_lookupTable[28] = QStringLiteral("Left-VentralDC");
336 m_lookupTable[30] = QStringLiteral("Left-vessel");
337 m_lookupTable[31] = QStringLiteral("Left-choroid-plexus");
338 m_lookupTable[41] = QStringLiteral("Right-Cerebral-White-Matter");
339 m_lookupTable[42] = QStringLiteral("Right-Cerebral-Cortex");
340 m_lookupTable[43] = QStringLiteral("Right-Lateral-Ventricle");
341 m_lookupTable[44] = QStringLiteral("Right-Inf-Lat-Vent");
342 m_lookupTable[46] = QStringLiteral("Right-Cerebellum-White-Matter");
343 m_lookupTable[47] = QStringLiteral("Right-Cerebellum-Cortex");
344 m_lookupTable[49] = QStringLiteral("Right-Thalamus");
345 m_lookupTable[50] = QStringLiteral("Right-Caudate");
346 m_lookupTable[51] = QStringLiteral("Right-Putamen");
347 m_lookupTable[52] = QStringLiteral("Right-Pallidum");
348 m_lookupTable[53] = QStringLiteral("Right-Hippocampus");
349 m_lookupTable[54] = QStringLiteral("Right-Amygdala");
350 m_lookupTable[58] = QStringLiteral("Right-Accumbens-area");
351 m_lookupTable[60] = QStringLiteral("Right-VentralDC");
352 m_lookupTable[62] = QStringLiteral("Right-vessel");
353 m_lookupTable[63] = QStringLiteral("Right-choroid-plexus");
354 m_lookupTable[72] = QStringLiteral("5th-Ventricle");
355 m_lookupTable[77] = QStringLiteral("WM-hypointensities");
356 m_lookupTable[80] = QStringLiteral("non-WM-hypointensities");
357 m_lookupTable[85] = QStringLiteral("Optic-Chiasm");
358 m_lookupTable[251] = QStringLiteral("CC_Posterior");
359 m_lookupTable[252] = QStringLiteral("CC_Mid_Posterior");
360 m_lookupTable[253] = QStringLiteral("CC_Central");
361 m_lookupTable[254] = QStringLiteral("CC_Mid_Anterior");
362 m_lookupTable[255] = QStringLiteral("CC_Anterior");
363
364 // Left hemisphere cortical parcellation (aparc, Desikan-Killiany atlas)
365 m_lookupTable[1001] = QStringLiteral("ctx-lh-bankssts");
366 m_lookupTable[1002] = QStringLiteral("ctx-lh-caudalanteriorcingulate");
367 m_lookupTable[1003] = QStringLiteral("ctx-lh-caudalmiddlefrontal");
368 m_lookupTable[1005] = QStringLiteral("ctx-lh-cuneus");
369 m_lookupTable[1006] = QStringLiteral("ctx-lh-entorhinal");
370 m_lookupTable[1007] = QStringLiteral("ctx-lh-fusiform");
371 m_lookupTable[1008] = QStringLiteral("ctx-lh-inferiorparietal");
372 m_lookupTable[1009] = QStringLiteral("ctx-lh-inferiortemporal");
373 m_lookupTable[1010] = QStringLiteral("ctx-lh-isthmuscingulate");
374 m_lookupTable[1011] = QStringLiteral("ctx-lh-lateraloccipital");
375 m_lookupTable[1012] = QStringLiteral("ctx-lh-lateralorbitofrontal");
376 m_lookupTable[1013] = QStringLiteral("ctx-lh-lingual");
377 m_lookupTable[1014] = QStringLiteral("ctx-lh-medialorbitofrontal");
378 m_lookupTable[1015] = QStringLiteral("ctx-lh-middletemporal");
379 m_lookupTable[1016] = QStringLiteral("ctx-lh-parahippocampal");
380 m_lookupTable[1017] = QStringLiteral("ctx-lh-paracentral");
381 m_lookupTable[1018] = QStringLiteral("ctx-lh-parsopercularis");
382 m_lookupTable[1019] = QStringLiteral("ctx-lh-parsorbitalis");
383 m_lookupTable[1020] = QStringLiteral("ctx-lh-parstriangularis");
384 m_lookupTable[1021] = QStringLiteral("ctx-lh-pericalcarine");
385 m_lookupTable[1022] = QStringLiteral("ctx-lh-postcentral");
386 m_lookupTable[1023] = QStringLiteral("ctx-lh-posteriorcingulate");
387 m_lookupTable[1024] = QStringLiteral("ctx-lh-precentral");
388 m_lookupTable[1025] = QStringLiteral("ctx-lh-precuneus");
389 m_lookupTable[1026] = QStringLiteral("ctx-lh-rostralanteriorcingulate");
390 m_lookupTable[1027] = QStringLiteral("ctx-lh-rostralmiddlefrontal");
391 m_lookupTable[1028] = QStringLiteral("ctx-lh-superiorfrontal");
392 m_lookupTable[1029] = QStringLiteral("ctx-lh-superiorparietal");
393 m_lookupTable[1030] = QStringLiteral("ctx-lh-superiortemporal");
394 m_lookupTable[1031] = QStringLiteral("ctx-lh-supramarginal");
395 m_lookupTable[1032] = QStringLiteral("ctx-lh-frontalpole");
396 m_lookupTable[1033] = QStringLiteral("ctx-lh-temporalpole");
397 m_lookupTable[1034] = QStringLiteral("ctx-lh-transversetemporal");
398 m_lookupTable[1035] = QStringLiteral("ctx-lh-insula");
399
400 // Right hemisphere cortical parcellation (aparc, Desikan-Killiany atlas)
401 m_lookupTable[2001] = QStringLiteral("ctx-rh-bankssts");
402 m_lookupTable[2002] = QStringLiteral("ctx-rh-caudalanteriorcingulate");
403 m_lookupTable[2003] = QStringLiteral("ctx-rh-caudalmiddlefrontal");
404 m_lookupTable[2005] = QStringLiteral("ctx-rh-cuneus");
405 m_lookupTable[2006] = QStringLiteral("ctx-rh-entorhinal");
406 m_lookupTable[2007] = QStringLiteral("ctx-rh-fusiform");
407 m_lookupTable[2008] = QStringLiteral("ctx-rh-inferiorparietal");
408 m_lookupTable[2009] = QStringLiteral("ctx-rh-inferiortemporal");
409 m_lookupTable[2010] = QStringLiteral("ctx-rh-isthmuscingulate");
410 m_lookupTable[2011] = QStringLiteral("ctx-rh-lateraloccipital");
411 m_lookupTable[2012] = QStringLiteral("ctx-rh-lateralorbitofrontal");
412 m_lookupTable[2013] = QStringLiteral("ctx-rh-lingual");
413 m_lookupTable[2014] = QStringLiteral("ctx-rh-medialorbitofrontal");
414 m_lookupTable[2015] = QStringLiteral("ctx-rh-middletemporal");
415 m_lookupTable[2016] = QStringLiteral("ctx-rh-parahippocampal");
416 m_lookupTable[2017] = QStringLiteral("ctx-rh-paracentral");
417 m_lookupTable[2018] = QStringLiteral("ctx-rh-parsopercularis");
418 m_lookupTable[2019] = QStringLiteral("ctx-rh-parsorbitalis");
419 m_lookupTable[2020] = QStringLiteral("ctx-rh-parstriangularis");
420 m_lookupTable[2021] = QStringLiteral("ctx-rh-pericalcarine");
421 m_lookupTable[2022] = QStringLiteral("ctx-rh-postcentral");
422 m_lookupTable[2023] = QStringLiteral("ctx-rh-posteriorcingulate");
423 m_lookupTable[2024] = QStringLiteral("ctx-rh-precentral");
424 m_lookupTable[2025] = QStringLiteral("ctx-rh-precuneus");
425 m_lookupTable[2026] = QStringLiteral("ctx-rh-rostralanteriorcingulate");
426 m_lookupTable[2027] = QStringLiteral("ctx-rh-rostralmiddlefrontal");
427 m_lookupTable[2028] = QStringLiteral("ctx-rh-superiorfrontal");
428 m_lookupTable[2029] = QStringLiteral("ctx-rh-superiorparietal");
429 m_lookupTable[2030] = QStringLiteral("ctx-rh-superiortemporal");
430 m_lookupTable[2031] = QStringLiteral("ctx-rh-supramarginal");
431 m_lookupTable[2032] = QStringLiteral("ctx-rh-frontalpole");
432 m_lookupTable[2033] = QStringLiteral("ctx-rh-temporalpole");
433 m_lookupTable[2034] = QStringLiteral("ctx-rh-transversetemporal");
434 m_lookupTable[2035] = QStringLiteral("ctx-rh-insula");
435}
fs_atlas_lookup.h
RAS-coordinate lookup against a FreeSurfer volumetric parcellation (e.g. aparc+aseg....
FSLIB
FreeSurfer surface, annotation and parcellation I/O for mne-cpp.
Definition connectivitysettings.h:64
FSLIB::FsAtlasLookup::labelAtRas
QString labelAtRas(const Eigen::Vector3f &ras) const
Look up the anatomical label at a RAS coordinate.
Definition fs_atlas_lookup.cpp:263
FSLIB::FsAtlasLookup::isLoaded
bool isLoaded() const
Check whether a parcellation volume has been loaded.
Definition fs_atlas_lookup.cpp:296
FSLIB::FsAtlasLookup::load
bool load(const QString &sParcellationPath)
Load a FreeSurfer MGH/MGZ volume parcellation.
Definition fs_atlas_lookup.cpp:56
FSLIB::FsAtlasLookup::labelsForPositions
QStringList labelsForPositions(const QVector< Eigen::Vector3f > &positions) const
Look up labels for multiple RAS positions.
Definition fs_atlas_lookup.cpp:285
Generated on for MNE-CPP by doxygen 1.16.1