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