MNE-CPP  0.1.9
A Framework for Electrophysiology
mne_bem.cpp
Go to the documentation of this file.
1 //=============================================================================================================
38 //=============================================================================================================
39 // INCLUDES
40 //=============================================================================================================
41 
42 #include "mne_bem.h"
43 
44 #include <utils/mnemath.h>
45 #include <utils/warp.h>
46 #include <fs/label.h>
47 
48 #include <iostream>
49 
50 //=============================================================================================================
51 // QT INCLUDES
52 //=============================================================================================================
53 
54 #include <QFile>
55 
56 //=============================================================================================================
57 // USED NAMESPACES
58 //=============================================================================================================
59 
60 using namespace UTILSLIB;
61 using namespace FSLIB;
62 using namespace MNELIB;
63 using namespace FIFFLIB;
64 using namespace Eigen;
65 
66 //=============================================================================================================
67 // DEFINE MEMBER METHODS
68 //=============================================================================================================
69 
70 MNEBem::MNEBem()
71 {
72 }
73 
74 //=============================================================================================================
75 
76 MNEBem::MNEBem(const MNEBem &p_MNEBem)
77 : m_qListBemSurface(p_MNEBem.m_qListBemSurface)
78 {
79 }
80 
81 //=============================================================================================================
82 
83 MNEBem::MNEBem(QIODevice &p_IODevice) //const MNEBem &p_MNEBem
84 //: m_qListBemSurface()
85 {
86  FiffStream::SPtr t_pStream(new FiffStream(&p_IODevice));
87 
88  if(!MNEBem::readFromStream(t_pStream, true, *this))
89  {
90  t_pStream->close();
91  std::cout << "Could not read the bem surfaces\n"; // ToDo throw error
92  //ToDo error(me,'Could not read the bem surfaces (%s)',mne_omit_first_line(lasterr));
93 // return false;
94  }
95 
96 // bool testStream =t_pStream->device()->isOpen();
97 }
98 
99 //=============================================================================================================
100 
102 {
103 }
104 
105 //=============================================================================================================
106 
108 {
109  m_qListBemSurface.clear();
110 }
111 
112 //=============================================================================================================
113 
114 bool MNEBem::readFromStream(FiffStream::SPtr& p_pStream, bool add_geom, MNEBem& p_Bem)
115 {
116  //
117  // Open the file, create directory
118  //
119  bool open_here = false;
120  QFile t_file;//ToDo TCPSocket;
121 
122  if (!p_pStream->device()->isOpen())
123  {
124  QString t_sFileName = p_pStream->streamName();
125 
126  t_file.setFileName(t_sFileName);
127  p_pStream = FiffStream::SPtr(new FiffStream(&t_file));
128  if(!p_pStream->open())
129  {
130  return false;
131  }
132  open_here = true;
133 // if(t_pDir)
134 // delete t_pDir;
135  }
136 
137  //
138  // Find all BEM surfaces
139  //
140 
141  QList<FiffDirNode::SPtr> bem = p_pStream->dirtree()->dir_tree_find(FIFFB_BEM);
142  if(bem.isEmpty())
143  {
144  qCritical() << "No BEM block found!";
145  if(open_here)
146  {
147  p_pStream->close();
148  }
149  return false;
150  }
151 
152  QList<FiffDirNode::SPtr> bemsurf = p_pStream->dirtree()->dir_tree_find(FIFFB_BEM_SURF);
153  if(bemsurf.isEmpty())
154  {
155  qCritical() << "No BEM surfaces found!";
156  if(open_here)
157  {
158  p_pStream->close();
159  }
160  return false;
161  }
162 
163  for(int k = 0; k < bemsurf.size(); ++k)
164  {
165  MNEBemSurface p_BemSurface;
166  printf("\tReading a BEM surface...");
167  MNEBem::readBemSurface(p_pStream, bemsurf[k], p_BemSurface);
168  p_BemSurface.addTriangleData();
169  if (add_geom)
170  {
171  p_BemSurface.addVertexNormals();
172  }
173  printf("\t[done]\n" );
174 
175  p_Bem.m_qListBemSurface.append(p_BemSurface);
176 // src(k) = this;
177  }
178 
179  printf("\t%d bem surfaces read\n", bemsurf.size());
180 
181  if(open_here)
182  {
183  p_pStream->close();
184  }
185  return true;
186 }
187 
188 //=============================================================================================================
189 
190 bool MNEBem::readBemSurface(FiffStream::SPtr& p_pStream, const FiffDirNode::SPtr &p_Tree, MNEBemSurface &p_BemSurface)
191 {
192  p_BemSurface.clear();
193 
194  FiffTag::SPtr t_pTag;
195 
196  //=====================================================================
197  if(!p_Tree->find_tag(p_pStream, FIFF_BEM_SURF_ID, t_pTag))
198  {
199  p_BemSurface.id = FIFFV_BEM_SURF_ID_UNKNOWN;
200  }
201  else
202  {
203  p_BemSurface.id = *t_pTag->toInt();
204  }
205 
206 // qDebug() << "Read BemSurface ID; type:" << t_pTag->getType() << "value:" << *t_pTag->toInt();
207 
208  //=====================================================================
209  if(!p_Tree->find_tag(p_pStream, FIFF_BEM_SIGMA, t_pTag))
210  {
211  p_BemSurface.sigma = 1.0;
212  }
213  else
214  {
215  p_BemSurface.sigma = *t_pTag->toFloat();
216  }
217 
218 // qDebug() <<
219 
220  //=====================================================================
221  if(!p_Tree->find_tag(p_pStream, FIFF_BEM_SURF_NNODE, t_pTag))
222  {
223  p_pStream->close();
224  std::cout << "np not found!";
225  return false;
226  }
227  else
228  {
229  p_BemSurface.np = *t_pTag->toInt();
230  }
231 
232 // qDebug() <<
233 
234  //=====================================================================
235  if(!p_Tree->find_tag(p_pStream, FIFF_BEM_SURF_NTRI, t_pTag))
236  {
237  p_pStream->close();
238  std::cout << "ntri not found!";
239  return false;
240  }
241  else
242  {
243  p_BemSurface.ntri = *t_pTag->toInt();
244  }
245 
246 // qDebug() <<
247 
248  //=====================================================================
249  if(!p_Tree->find_tag(p_pStream, FIFF_MNE_COORD_FRAME, t_pTag))
250  {
251  qWarning() << "FIFF_MNE_COORD_FRAME not found, trying FIFF_BEM_COORD_FRAME.";
252  if(!p_Tree->find_tag(p_pStream, FIFF_BEM_COORD_FRAME, t_pTag))
253  {
254  p_pStream->close();
255  std::cout << "Coordinate frame information not found."; //ToDo: throw error.
256  return false;
257  }
258  else
259  {
260  p_BemSurface.coord_frame = *t_pTag->toInt();
261  }
262  }
263  else
264  {
265  p_BemSurface.coord_frame = *t_pTag->toInt();
266  }
267 
268 // qDebug() <<
269 
270  //=====================================================================
271  //
272  // Vertices, normals, and triangles
273  //
274  //=====================================================================
275  if(!p_Tree->find_tag(p_pStream, FIFF_BEM_SURF_NODES, t_pTag))
276  {
277  p_pStream->close();
278  std::cout << "Vertex data not found."; //ToDo: throw error.
279  return false;
280  }
281 
282  p_BemSurface.rr = t_pTag->toFloatMatrix().transpose();
283  qint32 rows_rr = p_BemSurface.rr.rows();
284 
285  if (rows_rr != p_BemSurface.np)
286  {
287  p_pStream->close();
288  std::cout << "Vertex information is incorrect."; //ToDo: throw error.
289  return false;
290  }
291 
292 // qDebug() << "Surf Nodes; type:" << t_pTag->getType();
293 
294  //=====================================================================
295  if(!p_Tree->find_tag(p_pStream, FIFF_BEM_SURF_NORMALS, t_pTag))
296  {
297  if(!p_Tree->find_tag(p_pStream, FIFF_MNE_SOURCE_SPACE_NORMALS, t_pTag))
298  {
299  p_pStream->close();
300  std::cout << "Vertex normals not found."; //ToDo: throw error.
301  return false;
302  }
303 
304  p_BemSurface.nn = t_pTag->toFloatMatrix().transpose();
305  }
306  else
307  {
308  p_BemSurface.nn = t_pTag->toFloatMatrix().transpose();
309  }
310 
311  if (p_BemSurface.nn.rows() != p_BemSurface.np)
312  {
313  p_pStream->close();
314  std::cout << "Vertex normal information is incorrect."; //ToDo: throw error.
315  return false;
316  }
317 
318 // qDebug() << "Bem Vertex Normals; type:" << t_pTag->getType();
319 
320  //=====================================================================
321  if (p_BemSurface.ntri > 0)
322  {
323  if(!p_Tree->find_tag(p_pStream, FIFF_BEM_SURF_TRIANGLES, t_pTag))
324  {
325  if(!p_Tree->find_tag(p_pStream, FIFF_MNE_SOURCE_SPACE_TRIANGLES, t_pTag))
326  {
327  p_pStream->close();
328  std::cout << "Triangulation not found."; //ToDo: throw error.
329  return false;
330  }
331  else
332  {
333  p_BemSurface.tris = t_pTag->toIntMatrix().transpose();
334  p_BemSurface.tris -= MatrixXi::Constant(p_BemSurface.tris.rows(),3,1);//0 based indizes
335  }
336  }
337  else
338  {
339  p_BemSurface.tris = t_pTag->toIntMatrix().transpose();
340  p_BemSurface.tris -= MatrixXi::Constant(p_BemSurface.tris.rows(),3,1);//0 based indizes
341  }
342 
343  if (p_BemSurface.tris.rows() != p_BemSurface.ntri)
344  {
345  p_pStream->close();
346  std::cout << "Triangulation information is incorrect."; //ToDo: throw error.
347  return false;
348  }
349  }
350  else
351  {
352  MatrixXi p_defaultMatrix(0, 0);
353  p_BemSurface.tris = p_defaultMatrix;
354  }
355 
356  return true;
357 }
358 
359 //=============================================================================================================
360 
361 void MNEBem::write(QIODevice &p_IODevice)
362 {
363  //
364  // Open the file, create directory
365  //
366 
367  // Create the file and save the essentials
368  FiffStream::SPtr t_pStream = FiffStream::start_file(p_IODevice);
369  printf("Write BEM surface in %s...\n", t_pStream->streamName().toUtf8().constData());
370  this->writeToStream(t_pStream.data());
371  t_pStream->end_file();
372 }
373 
374 //=============================================================================================================
375 
377 {
378  p_pStream->start_block(FIFFB_BEM);
379  for(qint32 h = 0; h < m_qListBemSurface.size(); ++h)
380  {
381  printf("\tWrite a bem surface... ");
382  p_pStream->start_block(FIFFB_BEM_SURF);
383  m_qListBemSurface[h].writeToStream(p_pStream);
384  p_pStream->end_block(FIFFB_BEM_SURF);
385  printf("[done]\n");
386  }
387  printf("\t%d bem surfaces written\n", m_qListBemSurface.size());
388  p_pStream->end_block(FIFFB_BEM);
389 }
390 
391 //=============================================================================================================
392 
393 const MNEBemSurface& MNEBem::operator[] (qint32 idx) const
394 {
395  if (idx>=m_qListBemSurface.length())
396  {
397  qWarning("Warning: Required surface doesn't exist! Returning surface '0'.");
398  idx=0;
399  }
400  return m_qListBemSurface[idx];
401 }
402 
403 //=============================================================================================================
404 
406 {
407  if (idx >= m_qListBemSurface.length())
408  {
409  qWarning("Warning: Required surface doesn't exist! Returning surface '0'.");
410  idx = 0;
411  }
412  return m_qListBemSurface[idx];
413 }
414 
415 //=============================================================================================================
416 
418 {
419  this->m_qListBemSurface.append(surf);
420  return *this;
421 }
422 
423 //=============================================================================================================
424 
426 {
427  this->m_qListBemSurface.append(*surf);
428  return *this;
429 }
430 
431 //=============================================================================================================
432 
433 void MNEBem::warp(const MatrixXf & sLm, const MatrixXf &dLm)
434 {
435  Warp help;
436  QList<MatrixXf> vertList;
437  for (int i=0; i<this->m_qListBemSurface.size(); i++)
438  {
439  vertList.append(this->m_qListBemSurface[i].rr);
440  }
441 
442  help.calculate(sLm, dLm, vertList);
443 
444  for (int i=0; i<this->m_qListBemSurface.size(); i++)
445  {
446  this->m_qListBemSurface[i].rr = vertList.at(i);
447  }
448  return;
449 }
450 
451 //=============================================================================================================
452 
454 {
455  MatrixX3f vert;
456  for (int i=0; i<this->m_qListBemSurface.size(); i++)
457  {
458  vert = this->m_qListBemSurface[i].rr;
459  vert = trans.apply_trans(vert);
460  this->m_qListBemSurface[i].rr = vert;
461  }
462  return;
463 }
464 
465 //=============================================================================================================
466 
468 {
469  MatrixX3f vert;
470  for (int i=0; i<this->m_qListBemSurface.size(); i++)
471  {
472  vert = this->m_qListBemSurface[i].rr;
473  vert = trans.apply_inverse_trans(vert);
474  this->m_qListBemSurface[i].rr = vert;
475  }
476  return;
477 }
FIFFLIB::fiff_int_t id
Eigen::MatrixX3f nn
#define FIFF_BEM_SURF_TRIANGLES
Definition: fiff_file.h:739
FIFFLIB::fiff_float_t sigma
#define FIFF_BEM_SURF_NNODE
Definition: fiff_file.h:736
BEM surface provides geometry information.
QSharedPointer< FiffDirNode > SPtr
Definition: fiff_dir_node.h:77
MNEBem & operator<<(const MNEBemSurface &surf)
Definition: mne_bem.cpp:417
Eigen::MatrixX3f rr
#define FIFF_BEM_SURF_ID
Definition: fiff_file.h:734
Coordinate transformation description.
static bool readFromStream(FIFFLIB::FiffStream::SPtr &p_pStream, bool add_geom, MNEBem &p_Bem)
Definition: mne_bem.cpp:114
BEM descritpion.
Definition: mne_bem.h:89
void warp(const Eigen::MatrixXf &sLm, const Eigen::MatrixXf &dLm)
Definition: mne_bem.cpp:433
Eigen::MatrixXf calculate(const Eigen::MatrixXf &sLm, const Eigen::MatrixXf &dLm, const Eigen::MatrixXf &sVert)
FIFFLIB::fiff_int_t np
#define FIFF_MNE_SOURCE_SPACE_TRIANGLES
#define FIFF_MNE_SOURCE_SPACE_NORMALS
fiff_long_t start_block(fiff_int_t kind)
void clear()
Definition: mne_bem.cpp:107
MNEBem class declaration.
QSharedPointer< FiffTag > SPtr
Definition: fiff_tag.h:152
#define FIFFB_BEM
Definition: fiff_file.h:403
Thin Plate Spline Warp.
Definition: warp.h:72
FIFF File I/O routines.
Definition: fiff_stream.h:104
const MNEBemSurface & operator[](qint32 idx) const
Definition: mne_bem.cpp:393
MNEMath class declaration.
#define FIFFB_BEM_SURF
Definition: fiff_file.h:404
#define FIFF_BEM_SURF_NTRI
Definition: fiff_file.h:737
void write(QIODevice &p_IODevice)
Definition: mne_bem.cpp:361
void invtransform(const FIFFLIB::FiffCoordTrans &trans)
Definition: mne_bem.cpp:467
Eigen::MatrixX3f apply_trans(const Eigen::MatrixX3f &rr, bool do_move=true) const
void transform(const FIFFLIB::FiffCoordTrans &trans)
Definition: mne_bem.cpp:453
Eigen::MatrixX3f apply_inverse_trans(const Eigen::MatrixX3f &rr, bool do_move=true) const
#define FIFF_BEM_SURF_NORMALS
Definition: fiff_file.h:740
Eigen::MatrixX3i tris
QSharedPointer< FiffStream > SPtr
Definition: fiff_stream.h:107
FIFFLIB::fiff_int_t ntri
Warp class declaration.
Label class declaration.
void writeToStream(FIFFLIB::FiffStream *p_pStream)
Definition: mne_bem.cpp:376
fiff_long_t end_block(fiff_int_t kind, fiff_int_t next=FIFFV_NEXT_SEQ)
#define FIFF_BEM_COORD_FRAME
Definition: fiff_file.h:746
#define FIFF_BEM_SURF_NODES
Definition: fiff_file.h:738
static bool readBemSurface(FIFFLIB::FiffStream::SPtr &p_pStream, const FIFFLIB::FiffDirNode::SPtr &p_Tree, MNEBemSurface &p_BemSurface)
Definition: mne_bem.cpp:190
#define FIFF_BEM_SIGMA
Definition: fiff_file.h:747
#define FIFF_MNE_COORD_FRAME
FIFFLIB::fiff_int_t coord_frame