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