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