MNE-CPP  0.1.9
A Framework for Electrophysiology
surface.cpp
Go to the documentation of this file.
1 //=============================================================================================================
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "surface.h"
42 #include <utils/ioutils.h>
43 
44 #include <iostream>
45 
46 //=============================================================================================================
47 // QT INCLUDES
48 //=============================================================================================================
49 
50 #include <QFile>
51 #include <QDataStream>
52 #include <QTextStream>
53 
54 //=============================================================================================================
55 // USED NAMESPACES
56 //=============================================================================================================
57 
58 using namespace UTILSLIB;
59 using namespace FSLIB;
60 using namespace Eigen;
61 
62 //=============================================================================================================
63 // DEFINE MEMBER METHODS
64 //=============================================================================================================
65 
66 Surface::Surface()
67 : m_sFilePath("")
68 , m_sFileName("")
69 , m_iHemi(-1)
70 , m_sSurf("")
71 , m_vecOffset(Vector3f::Zero(3))
72 {
73 }
74 
75 //=============================================================================================================
76 
77 Surface::Surface(const QString& p_sFile)
78 : m_sFilePath("")
79 , m_sFileName("")
80 , m_iHemi(-1)
81 , m_sSurf("")
82 , m_vecOffset(Vector3f::Zero(3))
83 {
84  Surface::read(p_sFile, *this);
85 }
86 
87 //=============================================================================================================
88 
89 Surface::Surface(const QString &subject_id, qint32 hemi, const QString &surf, const QString &subjects_dir)
90 : m_sFilePath("")
91 , m_sFileName("")
92 , m_iHemi(-1)
93 , m_sSurf("")
94 , m_vecOffset(Vector3f::Zero(3))
95 {
96  Surface::read(subject_id, hemi, surf, subjects_dir, *this);
97 }
98 
99 //=============================================================================================================
100 
101 Surface::Surface(const QString &path, qint32 hemi, const QString &surf)
102 : m_sFilePath("")
103 , m_sFileName("")
104 , m_iHemi(-1)
105 , m_sSurf("")
106 , m_vecOffset(Vector3f::Zero(3))
107 {
108  Surface::read(path, hemi, surf, *this);
109 }
110 
111 //=============================================================================================================
112 
113 Surface::~Surface()
114 {
115 }
116 
117 //=============================================================================================================
118 
119 void Surface::clear()
120 {
121  m_sFilePath.clear();
122  m_sFileName.clear();
123  m_iHemi = -1;
124  m_sSurf.clear();
125  m_matRR.resize(0,3);
126  m_matTris.resize(0,3);
127  m_matNN.resize(0,3);
128  m_vecCurv.resize(0);
129 }
130 
131 //=============================================================================================================
132 
133 MatrixX3f Surface::compute_normals(const MatrixX3f& rr, const MatrixX3i& tris)
134 {
135  printf("\tcomputing normals\n");
136  // first, compute triangle normals
137  MatrixX3f r1(tris.rows(),3); MatrixX3f r2(tris.rows(),3); MatrixX3f r3(tris.rows(),3);
138 
139  for(qint32 i = 0; i < tris.rows(); ++i)
140  {
141  r1.row(i) = rr.row(tris(i, 0));
142  r2.row(i) = rr.row(tris(i, 1));
143  r3.row(i) = rr.row(tris(i, 2));
144  }
145 
146  MatrixX3f x = r2 - r1;
147  MatrixX3f y = r3 - r1;
148  MatrixX3f tri_nn(x.rows(),y.cols());
149  tri_nn.col(0) = x.col(1).cwiseProduct(y.col(2)) - x.col(2).cwiseProduct(y.col(1));
150  tri_nn.col(1) = x.col(2).cwiseProduct(y.col(0)) - x.col(0).cwiseProduct(y.col(2));
151  tri_nn.col(2) = x.col(0).cwiseProduct(y.col(1)) - x.col(1).cwiseProduct(y.col(0));
152 
153  // Triangle normals and areas
154  MatrixX3f tmp = tri_nn.cwiseProduct(tri_nn);
155  VectorXf normSize = tmp.rowwise().sum();
156  normSize = normSize.cwiseSqrt();
157 
158  for(qint32 i = 0; i < normSize.size(); ++i)
159  if(normSize(i) != 0)
160  tri_nn.row(i) /= normSize(i);
161 
162  MatrixX3f nn = MatrixX3f::Zero(rr.rows(), 3);
163 
164  for(qint32 p = 0; p < tris.rows(); ++p)
165  {
166  Vector3i verts = tris.row(p);
167  for(qint32 j = 0; j < verts.size(); ++j)
168  nn.row(verts(j)) = tri_nn.row(p);
169  }
170 
171  tmp = nn.cwiseProduct(nn);
172  normSize = tmp.rowwise().sum();
173  normSize = normSize.cwiseSqrt();
174 
175  for(qint32 i = 0; i < normSize.size(); ++i)
176  if(normSize(i) != 0)
177  nn.row(i) /= normSize(i);
178 
179  return nn;
180 }
181 
182 //=============================================================================================================
183 
184 bool Surface::read(const QString &subject_id, qint32 hemi, const QString &surf, const QString &subjects_dir, Surface &p_Surface, bool p_bLoadCurvature)
185 {
186  if(hemi != 0 && hemi != 1)
187  return false;
188 
189  QString p_sFile = QString("%1/%2/surf/%3.%4").arg(subjects_dir).arg(subject_id).arg(hemi == 0 ? "lh" : "rh").arg(surf);
190 
191  return read(p_sFile, p_Surface, p_bLoadCurvature);
192 }
193 
194 //=============================================================================================================
195 
196 bool Surface::read(const QString &path, qint32 hemi, const QString &surf, Surface &p_Surface, bool p_bLoadCurvature)
197 {
198  if(hemi != 0 && hemi != 1)
199  return false;
200 
201  QString p_sFile = QString("%1/%2.%3").arg(path).arg(hemi == 0 ? "lh" : "rh").arg(surf);
202 
203  return read(p_sFile, p_Surface, p_bLoadCurvature);
204 }
205 
206 //=============================================================================================================
207 
208 bool Surface::read(const QString &p_sFile, Surface &p_Surface, bool p_bLoadCurvature)
209 {
210  p_Surface.clear();
211 
212  QFile t_File(p_sFile);
213 
214  if (!t_File.open(QIODevice::ReadOnly))
215  {
216  printf("\tError: Couldn't open the surface file\n");
217  return false;
218  }
219 
220  printf("Reading surface...\n");
221 
222  //Strip file name and path
223  qint32 t_NameIdx = 0;
224  if(p_sFile.contains("lh."))
225  t_NameIdx = p_sFile.indexOf("lh.");
226  else if(p_sFile.contains("rh."))
227  t_NameIdx = p_sFile.indexOf("rh.");
228  else
229  return false;
230 
231  p_Surface.m_sFilePath = p_sFile.mid(0,t_NameIdx);
232  p_Surface.m_sFileName = p_sFile.mid(t_NameIdx,p_sFile.size()-t_NameIdx);
233 
234  QDataStream t_DataStream(&t_File);
235  t_DataStream.setByteOrder(QDataStream::BigEndian);
236 
237  //
238  // Magic numbers to identify QUAD and TRIANGLE files
239  //
240  // QUAD_FILE_MAGIC_NUMBER = (-1 & 0x00ffffff) ;
241  // NEW_QUAD_FILE_MAGIC_NUMBER = (-3 & 0x00ffffff) ;
242  //
243  qint32 NEW_QUAD_FILE_MAGIC_NUMBER = 16777213;
244  qint32 TRIANGLE_FILE_MAGIC_NUMBER = 16777214;
245  qint32 QUAD_FILE_MAGIC_NUMBER = 16777215;
246 
247  qint32 magic = IOUtils::fread3(t_DataStream);
248 
249  qint32 nvert = 0;
250  qint32 nquad = 0;
251  qint32 nface = 0;
252  MatrixXf verts;
253  MatrixXi faces;
254 
255  if(magic == QUAD_FILE_MAGIC_NUMBER || magic == NEW_QUAD_FILE_MAGIC_NUMBER)
256  {
257  nvert = IOUtils::fread3(t_DataStream);
258  nquad = IOUtils::fread3(t_DataStream);
259  if(magic == QUAD_FILE_MAGIC_NUMBER)
260  printf("\t%s is a quad file (nvert = %d nquad = %d)\n", p_sFile.toUtf8().constData(),nvert,nquad);
261  else
262  printf("\t%s is a new quad file (nvert = %d nquad = %d)\n", p_sFile.toUtf8().constData(),nvert,nquad);
263 
264  //vertices
265  verts.resize(nvert, 3);
266  if(magic == QUAD_FILE_MAGIC_NUMBER)
267  {
268  qint16 iVal;
269  for(qint32 i = 0; i < nvert; ++i)
270  {
271  for(qint32 j = 0; j < 3; ++j)
272  {
273  t_DataStream >> iVal;
274  IOUtils::swap_short(iVal);
275  verts(i,j) = ((float)iVal) / 100;
276  }
277  }
278  }
279  else
280  {
281  t_DataStream.readRawData((char *)verts.data(), nvert*3*sizeof(float));
282  for(qint32 i = 0; i < nvert; ++i)
283  for(qint32 j = 0; j < 3; ++j)
284  IOUtils::swap_floatp(&verts(i,j));
285  }
286 
287  MatrixXi quads = IOUtils::fread3_many(t_DataStream, nquad*4);
288  MatrixXi quads_new(4, nquad);
289  qint32 count = 0;
290  for(qint32 j = 0; j < quads.cols(); ++j)
291  {
292  for(qint32 i = 0; i < quads.rows(); ++i)
293  {
294  quads_new(i,j) = quads(count, 0);
295  ++count;
296  }
297  }
298  quads = quads_new.transpose();
299  //
300  // Face splitting follows
301  //
302  faces = MatrixXi::Zero(2*nquad,3);
303  for(qint32 k = 0; k < nquad; ++k)
304  {
305  RowVectorXi quad = quads.row(k);
306  if ((quad[0] % 2) == 0)
307  {
308  faces(nface,0) = quad[0];
309  faces(nface,1) = quad[1];
310  faces(nface,2) = quad[3];
311  ++nface;
312 
313  faces(nface,0) = quad[2];
314  faces(nface,1) = quad[3];
315  faces(nface,2) = quad[1];
316  ++nface;
317  }
318  else
319  {
320  faces(nface,0) = quad(0);
321  faces(nface,1) = quad(1);
322  faces(nface,2) = quad(2);
323  ++nface;
324 
325  faces(nface,0) = quad(0);
326  faces(nface,1) = quad(2);
327  faces(nface,2) = quad(3);
328  ++nface;
329  }
330  }
331  }
332  else if(magic == TRIANGLE_FILE_MAGIC_NUMBER)
333  {
334  QString s = t_File.readLine();
335  t_File.readLine();
336 
337  t_DataStream >> nvert;
338  t_DataStream >> nface;
339  IOUtils::swap_int(nvert);
340  IOUtils::swap_int(nface);
341 
342  printf("\t%s is a triangle file (nvert = %d ntri = %d)\n", p_sFile.toUtf8().constData(), nvert, nface);
343  printf("\t%s", s.toUtf8().constData());
344 
345  //vertices
346  verts.resize(3, nvert);
347  t_DataStream.readRawData((char *)verts.data(), nvert*3*sizeof(float));
348  for(qint32 i = 0; i < 3; ++i)
349  for(qint32 j = 0; j < nvert; ++j)
350  IOUtils::swap_floatp(&verts(i,j));
351 
352  //faces
353  faces.resize(nface, 3);
354  qint32 iVal;
355  for(qint32 i = 0; i < nface; ++i)
356  {
357  for(qint32 j = 0; j < 3; ++j)
358  {
359  t_DataStream >> iVal;
360  IOUtils::swap_int(iVal);
361  faces(i,j) = iVal;
362  }
363  }
364  }
365  else
366  {
367  qWarning("Bad magic number (%d) in surface file %s",magic,p_sFile.toUtf8().constData());
368  return false;
369  }
370 
371  verts.transposeInPlace();
372  verts.array() *= 0.001f;
373 
374  p_Surface.m_matRR = verts.block(0,0,verts.rows(),3);
375  p_Surface.m_matTris = faces.block(0,0,faces.rows(),3);
376 
377  //-> not needed since qglbuilder is doing that for us
378  p_Surface.m_matNN = compute_normals(p_Surface.m_matRR, p_Surface.m_matTris);
379 
380  // hemi info
381  if(t_File.fileName().contains("lh."))
382  p_Surface.m_iHemi = 0;
383  else if(t_File.fileName().contains("rh."))
384  p_Surface.m_iHemi = 1;
385  else
386  {
387  p_Surface.m_iHemi = -1;
388  return false;
389  }
390 
391  //Loaded surface
392  p_Surface.m_sSurf = t_File.fileName().mid((t_NameIdx+3),t_File.fileName().size() - (t_NameIdx+3));
393 
394  //Load curvature
395  if(p_bLoadCurvature)
396  {
397  QString t_sCurvatureFile = QString("%1%2.curv").arg(p_Surface.m_sFilePath).arg(p_Surface.m_iHemi == 0 ? "lh" : "rh");
398  printf("\t");
399  p_Surface.m_vecCurv = Surface::read_curv(t_sCurvatureFile);
400  }
401 
402  t_File.close();
403  printf("\tRead a surface with %d vertices from %s\n[done]\n",nvert,p_sFile.toUtf8().constData());
404 
405  return true;
406 }
407 
408 //=============================================================================================================
409 
410 VectorXf Surface::read_curv(const QString &p_sFileName)
411 {
412  VectorXf curv;
413 
414  printf("Reading curvature...");
415  QFile t_File(p_sFileName);
416 
417  if (!t_File.open(QIODevice::ReadOnly))
418  {
419  printf("\tError: Couldn't open the curvature file\n");
420  return curv;
421  }
422 
423  QDataStream t_DataStream(&t_File);
424  t_DataStream.setByteOrder(QDataStream::BigEndian);
425 
426  qint32 vnum = IOUtils::fread3(t_DataStream);
427  qint32 NEW_VERSION_MAGIC_NUMBER = 16777215;
428 
429  if(vnum == NEW_VERSION_MAGIC_NUMBER)
430  {
431  qint32 fnum, vals_per_vertex;
432  t_DataStream >> vnum;
433 
434  t_DataStream >> fnum;
435  t_DataStream >> vals_per_vertex;
436 
437  curv.resize(vnum, 1);
438  t_DataStream.readRawData((char *)curv.data(), vnum*sizeof(float));
439  for(qint32 i = 0; i < vnum; ++i)
440  IOUtils::swap_floatp(&curv(i));
441  }
442  else
443  {
444  qint32 fnum = IOUtils::fread3(t_DataStream);
445  Q_UNUSED(fnum)
446  qint16 iVal;
447  curv.resize(vnum, 1);
448  for(qint32 i = 0; i < vnum; ++i)
449  {
450  t_DataStream >> iVal;
451  IOUtils::swap_short(iVal);
452  curv(i) = ((float)iVal) / 100;
453  }
454  }
455  t_File.close();
456 
457  printf("[done]\n");
458 
459  return curv;
460 }
FSLIB::Surface::read_curv
static Eigen::VectorXf read_curv(const QString &p_sFileName)
Definition: surface.cpp:410
FSLIB::Surface::Surface
Surface()
Definition: surface.cpp:66
FSLIB::Surface::hemi
qint32 hemi() const
Definition: surface.h:301
k
int k
Definition: fiff_tag.cpp:322
FSLIB::Surface::rr
const Eigen::MatrixX3f & rr() const
Definition: surface.h:322
FSLIB::Surface::nn
const Eigen::MatrixX3f & nn() const
Definition: surface.h:336
FSLIB::Surface::tris
const Eigen::MatrixX3i & tris() const
Definition: surface.h:329
FSLIB::Surface::~Surface
~Surface()
Definition: surface.cpp:113
FSLIB::Surface::surf
QString surf() const
Definition: surface.h:315
FSLIB::Surface::read
static bool read(const QString &subject_id, qint32 hemi, const QString &surf, const QString &subjects_dir, Surface &p_Surface, bool p_bLoadCurvature=true)
Definition: surface.cpp:184
FSLIB::Surface::compute_normals
static Eigen::MatrixX3f compute_normals(const Eigen::MatrixX3f &rr, const Eigen::MatrixX3i &tris)
Definition: surface.cpp:133
surface.h
Surface class declaration.
FSLIB::Surface
FreeSurfer surface mesh.
Definition: surface.h:75
FSLIB::Surface::curv
const Eigen::VectorXf & curv() const
Definition: surface.h:343
ioutils.h
IOUtils class declaration.
FSLIB::Surface::clear
void clear()
Definition: surface.cpp:119
Generated on Wed Aug 9 2023 17:55:41 for MNE-CPP by   doxygen 1.8.17