MNE-CPP  0.1.9
A Framework for Electrophysiology
mne_sourcespace.cpp
Go to the documentation of this file.
1 //=============================================================================================================
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "mne_sourcespace.h"
42 
43 #include <utils/mnemath.h>
44 #include <fs/label.h>
45 
46 #include <iostream>
47 
48 //=============================================================================================================
49 // QT INCLUDES
50 //=============================================================================================================
51 
52 #include <QFile>
53 
54 //=============================================================================================================
55 // USED NAMESPACES
56 //=============================================================================================================
57 
58 using namespace UTILSLIB;
59 using namespace FSLIB;
60 using namespace MNELIB;
61 using namespace Eigen;
62 using namespace FIFFLIB;
63 
64 //=============================================================================================================
65 // DEFINE MEMBER METHODS
66 //=============================================================================================================
67 
68 MNESourceSpace::MNESourceSpace()
69 {
70 }
71 
72 //=============================================================================================================
73 
74 MNESourceSpace::MNESourceSpace(const MNESourceSpace &p_MNESourceSpace)
75 : m_qListHemispheres(p_MNESourceSpace.m_qListHemispheres)
76 {
77 }
78 
79 //=============================================================================================================
80 
82 {
83 }
84 
85 //=============================================================================================================
86 
88 {
89  m_qListHemispheres.clear();
90 }
91 
92 //=============================================================================================================
93 
94 QList<VectorXi> MNESourceSpace::get_vertno() const
95 {
96  QList<VectorXi> p_vertices;
97  for(qint32 i = 0; i < m_qListHemispheres.size(); ++i)
98  p_vertices.push_back(m_qListHemispheres[i].vertno);
99  return p_vertices;
100 }
101 
102 //=============================================================================================================
103 
104 QList<VectorXi> MNESourceSpace::label_src_vertno_sel(const Label &p_label, VectorXi &src_sel) const
105 {
106 // if(src[0].['type'] != 'surf')
107 // return Exception('Label are only supported with surface source spaces')
108 
109  QList<VectorXi> vertno;
110  vertno << this->m_qListHemispheres[0].vertno << this->m_qListHemispheres[1].vertno;
111 
112  if (p_label.hemi == 0) //lh
113  {
114  VectorXi vertno_sel = MNEMath::intersect(vertno[0], p_label.vertices, src_sel);
115  vertno[0] = vertno_sel;
116  vertno[1] = VectorXi();
117  }
118  else if (p_label.hemi == 1) //rh
119  {
120  VectorXi vertno_sel = MNEMath::intersect(vertno[1], p_label.vertices, src_sel);
121  src_sel.array() += p_label.vertices.size();
122  vertno[0] = VectorXi();
123  vertno[1] = vertno_sel;
124  }
125 
126 // if (p_label.hemi == 0) //lh
127 // {
128 // VectorXi vertno_sel = MNEMath::intersect(vertno[0], p_label.vertices[0], src_sel);
129 // vertno[0] = vertno_sel;
130 // vertno[1] = VectorXi();
131 // }
132 // else if (p_label.hemi == 1) //rh
133 // {
134 // VectorXi vertno_sel = MNEMath::intersect(vertno[1], p_label.vertices[1], src_sel);
135 // src_sel.array() += p_label.vertices[0].size();
136 // vertno[0] = VectorXi();
137 // vertno[1] = vertno_sel;
138 // }
139 // else if (p_label.hemi == 2) //both
140 // {
141 // VectorXi src_sel_lh, src_sel_rh;
142 // VectorXi vertno_sel_lh = MNEMath::intersect(vertno[0], p_label.vertices[0], src_sel_lh);
143 // VectorXi vertno_sel_rh = MNEMath::intersect(vertno[1], p_label.vertices[1], src_sel_rh);
144 // src_sel.resize(src_sel_lh.size() + src_sel_rh.size());
145 // src_sel.block(0,0,src_sel_lh.size(),1) = src_sel_lh;
146 // src_sel.block(src_sel_lh.size(),0,src_sel_rh.size(),1) = src_sel_rh;
147 // vertno[0] = vertno_sel_lh;
148 // vertno[0] = vertno_sel_rh;
149 // }
150  else
151  {
152  qWarning("Unknown hemisphere type\n");
153  vertno[0] = VectorXi::Zero(0);
154  vertno[1] = VectorXi::Zero(0);
155  }
156 
157  return vertno;
158 }
159 
160 //=============================================================================================================
161 
162 MNESourceSpace MNESourceSpace::pick_regions(const QList<Label> &p_qListLabels) const
163 {
164  Q_UNUSED(p_qListLabels);
165 
166  MNESourceSpace selectedSrc(*this);
167 
168  for(qint32 h = 0; h < 2; ++h)
169  {
170  VectorXi selVertices;
171 
172  //get vertices indeces for new selection
173  qint32 iSize = 0;
174  for(qint32 i = 0; i < p_qListLabels.size(); ++i)
175  {
176  if(p_qListLabels[i].hemi == h)
177  {
178  VectorXi currentSelection;
179 
180  MNEMath::intersect(m_qListHemispheres[h].vertno, p_qListLabels[i].vertices, currentSelection);
181 
182  selVertices.conservativeResize(iSize+currentSelection.size());
183  selVertices.block(iSize,0,currentSelection.size(),1) = currentSelection;
184  iSize = selVertices.size();
185  }
186  }
187 
188  MNEMath::sort(selVertices, false);
189 
190  VectorXi newVertno(selVertices.size());
191 
192  selectedSrc.m_qListHemispheres[h].inuse = VectorXi::Zero(selectedSrc.m_qListHemispheres[h].np);
193 
194  for(qint32 i = 0; i < selVertices.size(); ++i)
195  {
196  selectedSrc.m_qListHemispheres[h].inuse[selVertices[i]] = 1;
197  newVertno[i] = this->m_qListHemispheres[h].vertno[selVertices[i]];
198  }
199 
200  selectedSrc.m_qListHemispheres[h].nuse = selVertices.size();
201  selectedSrc.m_qListHemispheres[h].vertno = newVertno;
202 
203  //
204  // Tris
205  //
206  VectorXi idx_select = VectorXi::Zero(this->m_qListHemispheres[h].use_tris.rows());
207  for(qint32 i = 0; i < 3; ++i)
208  {
209  VectorXi tri_dim = this->m_qListHemispheres[h].use_tris.col(i);
210  VectorXi idx_dim;
211  MNEMath::intersect(tri_dim, newVertno, idx_dim);
212 
213  for(qint32 j = 0; j < idx_dim.size(); ++j)
214  idx_select[idx_dim[j]] = 1;
215  }
216 
217  qint32 countSel = 0;
218  for(qint32 i = 0; i < idx_select.size(); ++i)
219  if(idx_select[i] == 1)
220  ++countSel;
221 
222  selectedSrc.m_qListHemispheres[h].nuse_tri = countSel;
223 
224  MatrixX3i use_tris_new(countSel,3);
225  MatrixX3d use_tri_cent_new(countSel,3);
226  MatrixX3d use_tri_nn_new(countSel,3);
227  VectorXd use_tri_area_new(countSel);
228 
229  countSel = 0;
230  for(qint32 i = 0; i < idx_select.size(); ++i)
231  {
232  if(idx_select[i] == 1)
233  {
234  use_tris_new.row(countSel) = this->m_qListHemispheres[h].use_tris.row(i);
235  use_tri_cent_new.row(countSel) = this->m_qListHemispheres[h].use_tri_cent.row(i);
236  use_tri_nn_new.row(countSel) = this->m_qListHemispheres[h].use_tri_nn.row(i);
237  use_tri_area_new[countSel] = this->m_qListHemispheres[h].use_tri_area[i];
238  ++countSel;
239  }
240  }
241 
242  selectedSrc.m_qListHemispheres[h].use_tris = use_tris_new;
243  selectedSrc.m_qListHemispheres[h].use_tri_cent = use_tri_cent_new;
244  selectedSrc.m_qListHemispheres[h].use_tri_nn = use_tri_nn_new;
245  selectedSrc.m_qListHemispheres[h].use_tri_area = use_tri_area_new;
246  }
247 
248  return selectedSrc;
249 }
250 
251 //=============================================================================================================
252 
254  bool add_geom,
255  MNESourceSpace& p_SourceSpace)
256 {
257 // if (p_pSourceSpace != NULL)
258 // delete p_pSourceSpace;
259  p_SourceSpace = MNESourceSpace();
260 
261  //
262  // Open the file, create directory
263  //
264  bool open_here = false;
265  QFile t_file;//ToDo TCPSocket;
266 
267  if (!p_pStream->device()->isOpen())
268  {
269  QString t_sFileName = p_pStream->streamName();
270 
271  t_file.setFileName(t_sFileName);
272  p_pStream = FiffStream::SPtr(new FiffStream(&t_file));
273  if(!p_pStream->open())
274  return false;
275  open_here = true;
276 // if(t_pDir)
277 // delete t_pDir;
278  }
279  //
280  // Find all source spaces
281  //
282  QList<FiffDirNode::SPtr> spaces = p_pStream->dirtree()->dir_tree_find(FIFFB_MNE_SOURCE_SPACE);
283  if (spaces.size() == 0)
284  {
285  if(open_here)
286  p_pStream->close();
287  std::cout << "No source spaces found";
288  return false;
289  }
290 
291  for(int k = 0; k < spaces.size(); ++k)
292  {
293  MNEHemisphere p_Hemisphere;
294  printf("\tReading a source space...");
295  MNESourceSpace::read_source_space(p_pStream, spaces[k], p_Hemisphere);
296  printf("\t[done]\n" );
297  if (add_geom)
298  complete_source_space_info(p_Hemisphere);
299 
300  p_SourceSpace.m_qListHemispheres.append(p_Hemisphere);
301 
302 // src(k) = this;
303  }
304 
305  printf("\t%d source spaces read\n", spaces.size());
306 
307  if(open_here)
308  p_pStream->close();
309 
310  return true;
311 }
312 
313 //=============================================================================================================
314 
316 {
317  double xave = p_Hemisphere.rr.col(0).sum();
318 
319  qint32 hemi;
320  if (xave < 0)
321  hemi = FIFFV_MNE_SURF_LEFT_HEMI;
322  else
323  hemi = FIFFV_MNE_SURF_RIGHT_HEMI;
324 
325  return hemi;
326 }
327 
328 //=============================================================================================================
329 
331 {
332  for(int k = 0; k < this->m_qListHemispheres.size(); ++k)
333  {
334  if(!this->m_qListHemispheres[k].transform_hemisphere_to(dest,trans))
335  {
336  printf("Could not transform source space.\n");
337  return false;
338  }
339  }
340  return true;
341 }
342 
343 //=============================================================================================================
344 
345 bool MNESourceSpace::read_source_space(FiffStream::SPtr& p_pStream, const FiffDirNode::SPtr& p_Tree, MNEHemisphere& p_Hemisphere)
346 {
347  p_Hemisphere.clear();
348 
349  FiffTag::SPtr t_pTag;
350 
351  //=====================================================================
352  if(!p_Tree->find_tag(p_pStream, FIFF_MNE_SOURCE_SPACE_ID, t_pTag))
353  p_Hemisphere.id = FIFFV_MNE_SURF_UNKNOWN;
354  else
355  p_Hemisphere.id = *t_pTag->toInt();
356 
357 // qDebug() << "Read SourceSpace ID; type:" << t_pTag->getType() << "value:" << *t_pTag->toInt();
358 
359  //=====================================================================
360  if(!p_Tree->find_tag(p_pStream, FIFF_MNE_SOURCE_SPACE_NPOINTS, t_pTag))
361  {
362  p_pStream->close();
363  std::cout << "error: Number of vertices not found."; //ToDo: throw error.
364  return false;
365  }
366 // qDebug() << "Number of vertice; type:" << t_pTag->getType() << "value:" << *t_pTag->toInt();
367  p_Hemisphere.np = *t_pTag->toInt();
368 
369  //=====================================================================
370  if(!p_Tree->find_tag(p_pStream, FIFF_BEM_SURF_NTRI, t_pTag))
371  {
372  if(!p_Tree->find_tag(p_pStream, FIFF_MNE_SOURCE_SPACE_NTRI, t_pTag))
373  p_Hemisphere.ntri = 0;
374  else
375  p_Hemisphere.ntri = *t_pTag->toInt();
376  }
377  else
378  {
379  p_Hemisphere.ntri = *t_pTag->toInt();
380  }
381 // qDebug() << "Number of Tris; type:" << t_pTag->getType() << "value:" << *t_pTag->toInt();
382 
383  //=====================================================================
384  if(!p_Tree->find_tag(p_pStream, FIFF_MNE_COORD_FRAME, t_pTag))
385  {
386  p_pStream->close();
387  std::cout << "Coordinate frame information not found."; //ToDo: throw error.
388  return false;
389  }
390  p_Hemisphere.coord_frame = *t_pTag->toInt();
391 // qDebug() << "Coord Frame; type:" << t_pTag->getType() << "value:" << *t_pTag->toInt();
392 
393  //=====================================================================
394  //
395  // Vertices, normals, and triangles
396  //
397  if(!p_Tree->find_tag(p_pStream, FIFF_MNE_SOURCE_SPACE_POINTS, t_pTag))
398  {
399  p_pStream->close();
400  std::cout << "Vertex data not found."; //ToDo: throw error.
401  return false;
402  }
403 
404  p_Hemisphere.rr = t_pTag->toFloatMatrix().transpose();
405  qint32 rows_rr = p_Hemisphere.rr.rows();
406 // qDebug() << "last element rr: " << p_Hemisphere.rr(rows_rr-1, 0) << p_Hemisphere.rr(rows_rr-1, 1) << p_Hemisphere.rr(rows_rr-1, 2);
407 
408  if (rows_rr != p_Hemisphere.np)
409  {
410  p_pStream->close();
411  std::cout << "Vertex information is incorrect."; //ToDo: throw error.
412  return false;
413  }
414 // qDebug() << "Source Space Points; type:" << t_pTag->getType();
415 
416  //=====================================================================
417  if(!p_Tree->find_tag(p_pStream, FIFF_MNE_SOURCE_SPACE_NORMALS, t_pTag))
418  {
419  p_pStream->close();
420  std::cout << "Vertex normals not found."; //ToDo: throw error.
421  return false;
422  }
423 
424  p_Hemisphere.nn = t_pTag->toFloatMatrix().transpose();
425  qint32 rows_nn = p_Hemisphere.nn.rows();
426 
427  if (rows_nn != p_Hemisphere.np)
428  {
429  p_pStream->close();
430  std::cout << "Vertex normal information is incorrect."; //ToDo: throw error.
431  return false;
432  }
433 // qDebug() << "Source Space Normals; type:" << t_pTag->getType();
434 
435  //=====================================================================
436  if (p_Hemisphere.ntri > 0)
437  {
438  if(!p_Tree->find_tag(p_pStream, FIFF_BEM_SURF_TRIANGLES, t_pTag))
439  {
440  if(!p_Tree->find_tag(p_pStream, FIFF_MNE_SOURCE_SPACE_TRIANGLES, t_pTag))
441  {
442  p_pStream->close();
443  std::cout << "Triangulation not found."; //ToDo: throw error.
444  return false;
445  }
446  else
447  {
448  p_Hemisphere.tris = t_pTag->toIntMatrix().transpose();
449  p_Hemisphere.tris -= MatrixXi::Constant(p_Hemisphere.tris.rows(),3,1);//0 based indizes
450  }
451  }
452  else
453  {
454  p_Hemisphere.tris = t_pTag->toIntMatrix().transpose();
455  p_Hemisphere.tris -= MatrixXi::Constant(p_Hemisphere.tris.rows(),3,1);//0 based indizes
456  }
457  if (p_Hemisphere.tris.rows() != p_Hemisphere.ntri)
458  {
459  p_pStream->close();
460  std::cout << "Triangulation information is incorrect."; //ToDo: throw error.
461  return false;
462  }
463  }
464  else
465  {
466  MatrixXi p_defaultMatrix(0, 0);
467  p_Hemisphere.tris = p_defaultMatrix;
468  }
469 // qDebug() << "Triangles; type:" << t_pTag->getType() << "rows:" << p_Hemisphere.tris.rows() << "cols:" << p_Hemisphere.tris.cols();
470 
471  //
472  // Which vertices are active
473  //
474  if(!p_Tree->find_tag(p_pStream, FIFF_MNE_SOURCE_SPACE_NUSE, t_pTag))
475  {
476  p_Hemisphere.nuse = 0;
477  p_Hemisphere.inuse = VectorXi::Zero(p_Hemisphere.nuse);
478  VectorXi p_defaultVector;
479  p_Hemisphere.vertno = p_defaultVector;
480  }
481  else
482  {
483  p_Hemisphere.nuse = *t_pTag->toInt();
484  if(!p_Tree->find_tag(p_pStream, FIFF_MNE_SOURCE_SPACE_SELECTION, t_pTag))
485  {
486  p_pStream->close();
487  std::cout << "Source selection information missing."; //ToDo: throw error.
488  return false;
489  }
490  p_Hemisphere.inuse = VectorXi(Map<VectorXi>(t_pTag->toInt(), t_pTag->size()/4, 1));//use copy constructor, for the sake of easy memory management
491 
492  p_Hemisphere.vertno = VectorXi::Zero(p_Hemisphere.nuse);
493  if (p_Hemisphere.inuse.rows() != p_Hemisphere.np)
494  {
495  p_pStream->close();
496  std::cout << "Incorrect number of entries in source space selection."; //ToDo: throw error.
497  return false;
498  }
499  int pp = 0;
500  for (int p = 0; p < p_Hemisphere.np; ++p)
501  {
502  if(p_Hemisphere.inuse(p) == 1)
503  {
504  p_Hemisphere.vertno(pp) = p;
505  ++pp;
506  }
507  }
508  }
509 // qDebug() << "Vertices; type:" << t_pTag->getType() << "nuse:" << p_Hemisphere.nuse;
510 
511  //
512  // Use triangulation
513  //
514  FiffTag::SPtr t_pTag1;
515  FiffTag::SPtr t_pTag2;
516  if(!p_Tree->find_tag(p_pStream, FIFF_MNE_SOURCE_SPACE_NUSE_TRI, t_pTag1) || !p_Tree->find_tag(p_pStream, FIFF_MNE_SOURCE_SPACE_USE_TRIANGLES, t_pTag2))
517  {
518  MatrixX3i p_defaultMatrix;
519  p_Hemisphere.nuse_tri = 0;
520  p_Hemisphere.use_tris = p_defaultMatrix;
521  }
522  else
523  {
524  p_Hemisphere.nuse_tri = *t_pTag1->toInt();
525  p_Hemisphere.use_tris = t_pTag2->toIntMatrix().transpose();
526  p_Hemisphere.use_tris -= MatrixXi::Constant(p_Hemisphere.use_tris.rows(),3,1); //0 based indizes
527  }
528 // qDebug() << "triangulation; type:" << t_pTag2->getType() << "use_tris:" << p_Hemisphere.use_tris.rows()<< "x" << p_Hemisphere.use_tris.cols();
529 
530  //
531  // Patch-related information
532  //
533  if(!p_Tree->find_tag(p_pStream, FIFF_MNE_SOURCE_SPACE_NEAREST, t_pTag1) || !p_Tree->find_tag(p_pStream, FIFF_MNE_SOURCE_SPACE_NEAREST_DIST, t_pTag2))
534  {
535  VectorXi p_defaultVector;
536  p_Hemisphere.nearest = p_defaultVector;
537  VectorXd p_defaultFloatVector;
538  p_Hemisphere.nearest_dist = p_defaultFloatVector;
539  }
540  else
541  {
542  //res.nearest = tag1.data + 1;
543  p_Hemisphere.nearest = VectorXi(Map<VectorXi>(t_pTag1->toInt(), t_pTag1->size()/4, 1));//use copy constructor, for the sake of easy memory management
544  p_Hemisphere.nearest_dist = VectorXd((Map<VectorXf>(t_pTag2->toFloat(), t_pTag1->size()/4, 1)).cast<double>());//use copy constructor, for the sake of easy memory management
545  }
546 
547 // patch_info(p_Hemisphere.nearest, p_Hemisphere.pinfo);
548  if (patch_info(p_Hemisphere))
549  printf("\tPatch information added...");
550 
551  //
552  // Distances
553  //
554 // if(p_Hemisphere.dist)
555 // delete p_Hemisphere.dist;
556  if(!p_Tree->find_tag(p_pStream, FIFF_MNE_SOURCE_SPACE_DIST, t_pTag1) || !p_Tree->find_tag(p_pStream, FIFF_MNE_SOURCE_SPACE_DIST_LIMIT, t_pTag2))
557  {
558  p_Hemisphere.dist = SparseMatrix<double>();//NULL;
559  p_Hemisphere.dist_limit = 0;
560  }
561  else
562  {
563  p_Hemisphere.dist = t_pTag1->toSparseFloatMatrix();
564  p_Hemisphere.dist_limit = *t_pTag2->toFloat(); //ToDo Check if this is realy always a float and not a matrix
565  //
566  // Add the upper triangle
567  //
568  SparseMatrix<double> distT = p_Hemisphere.dist.transpose();
569  p_Hemisphere.dist += distT;
570  }
571 
572  return true;
573 }
574 
575 //=============================================================================================================
576 
577 bool MNESourceSpace::patch_info(MNEHemisphere &p_Hemisphere)//VectorXi& nearest, QList<VectorXi>& pinfo)
578 {
579  if (p_Hemisphere.nearest.rows() == 0)
580  {
581  p_Hemisphere.pinfo.clear();
582  p_Hemisphere.patch_inds = VectorXi();
583  return false;
584  }
585 
586  printf("\tComputing patch statistics...");
587 
588  std::vector< std::pair<int,int> > t_vIndn;
589 
590  for(qint32 i = 0; i < p_Hemisphere.nearest.rows(); ++i)
591  {
592  std::pair<int,int> t_pair(i, p_Hemisphere.nearest(i));
593  t_vIndn.push_back(t_pair);
594  }
595  std::sort(t_vIndn.begin(),t_vIndn.end(), MNEMath::compareIdxValuePairSmallerThan<int> );
596 
597  VectorXi nearest_sorted(t_vIndn.size());
598 
599  qint32 current = 0;
600  std::vector<qint32> t_vfirsti;
601  t_vfirsti.push_back(current);
602  std::vector<qint32> t_vlasti;
603 
604  for(quint32 i = 0; i < t_vIndn.size(); ++i)
605  {
606  nearest_sorted[i] = t_vIndn[i].second;
607  if (t_vIndn[current].second != t_vIndn[i].second)
608  {
609  current = i;
610  t_vlasti.push_back(i-1);
611  t_vfirsti.push_back(current);
612  }
613  }
614  t_vlasti.push_back(static_cast<int>(t_vIndn.size()-1));
615 
616  for(quint32 k = 0; k < t_vfirsti.size(); ++k)
617  {
618  std::vector<int> t_vIndex;
619 
620  for(int l = t_vfirsti[k]; l <= t_vlasti[k]; ++l)
621  t_vIndex.push_back(t_vIndn[l].first);
622 
623  std::sort(t_vIndex.begin(),t_vIndex.end());
624 
625  int* t_pV = &t_vIndex[0];
626  Eigen::Map<Eigen::VectorXi> t_vPInfo(t_pV, t_vIndex.size());
627 
628  p_Hemisphere.pinfo.append(t_vPInfo);
629  }
630 
631  // compute patch indices of the in-use source space vertices
632  std::vector<qint32> patch_verts;
633  patch_verts.reserve(t_vlasti.size());
634  for(quint32 i = 0; i < t_vlasti.size(); ++i)
635  patch_verts.push_back(nearest_sorted[t_vlasti[i]]);
636 
637  p_Hemisphere.patch_inds.resize(p_Hemisphere.vertno.size());
638  std::vector<qint32>::iterator it;
639  for(qint32 i = 0; i < p_Hemisphere.vertno.size(); ++i)
640  {
641  it = std::find(patch_verts.begin(), patch_verts.end(), p_Hemisphere.vertno[i]);
642  p_Hemisphere.patch_inds[i] = it-patch_verts.begin();
643  }
644 
645  return true;
646 }
647 
648 //=============================================================================================================
649 
650 bool MNESourceSpace::complete_source_space_info(MNEHemisphere& p_Hemisphere)
651 {
652  //
653  // Main triangulation
654  //
655  printf("\tCompleting triangulation info...");
656  p_Hemisphere.tri_cent = MatrixX3d::Zero(p_Hemisphere.ntri,3);
657  p_Hemisphere.tri_nn = MatrixX3d::Zero(p_Hemisphere.ntri,3);
658  p_Hemisphere.tri_area = VectorXd::Zero(p_Hemisphere.ntri);
659 
660  Matrix3d r;
661  Vector3d a, b;
662  int k = 0;
663  float size = 0;
664  for (qint32 i = 0; i < p_Hemisphere.ntri; ++i)
665  {
666  for ( qint32 j = 0; j < 3; ++j)
667  {
668  k = p_Hemisphere.tris(i, j);
669 
670  r(j,0) = p_Hemisphere.rr(k, 0);
671  r(j,1) = p_Hemisphere.rr(k, 1);
672  r(j,2) = p_Hemisphere.rr(k, 2);
673 
674  p_Hemisphere.tri_cent(i, 0) += p_Hemisphere.rr(k, 0);
675  p_Hemisphere.tri_cent(i, 1) += p_Hemisphere.rr(k, 1);
676  p_Hemisphere.tri_cent(i, 2) += p_Hemisphere.rr(k, 2);
677  }
678  p_Hemisphere.tri_cent.row(i) /= 3.0f;
679 
680  //cross product {cross((r2-r1),(r3-r1))}
681  a = r.row(1) - r.row(0 );
682  b = r.row(2) - r.row(0);
683  p_Hemisphere.tri_nn(i,0) = a(1)*b(2)-a(2)*b(1);
684  p_Hemisphere.tri_nn(i,1) = a(2)*b(0)-a(0)*b(2);
685  p_Hemisphere.tri_nn(i,2) = a(0)*b(1)-a(1)*b(0);
686 
687  //area
688  size = p_Hemisphere.tri_nn.row(i)*p_Hemisphere.tri_nn.row(i).transpose();
689  size = std::pow(size, 0.5f );
690 
691  p_Hemisphere.tri_area(i) = size/2.0f;
692  p_Hemisphere.tri_nn.row(i) /= size;
693 
694  }
695  printf("[done]\n");
696 
697 // qDebug() << "p_Hemisphere.tri_cent:" << p_Hemisphere.tri_cent(0,0) << p_Hemisphere.tri_cent(0,1) << p_Hemisphere.tri_cent(0,2);
698 // qDebug() << "p_Hemisphere.tri_cent:" << p_Hemisphere.tri_cent(2,0) << p_Hemisphere.tri_cent(2,1) << p_Hemisphere.tri_cent(2,2);
699 
700 // qDebug() << "p_Hemisphere.tri_nn:" << p_Hemisphere.tri_nn(0,0) << p_Hemisphere.tri_nn(0,1) << p_Hemisphere.tri_nn(0,2);
701 // qDebug() << "p_Hemisphere.tri_nn:" << p_Hemisphere.tri_nn(2,0) << p_Hemisphere.tri_nn(2,1) << p_Hemisphere.tri_nn(2,2);
702 
703  //
704  // Selected triangles
705  //
706  printf("\tCompleting selection triangulation info...");
707  if (p_Hemisphere.nuse_tri > 0)
708  {
709  p_Hemisphere.use_tri_cent = MatrixX3d::Zero(p_Hemisphere.nuse_tri,3);
710  p_Hemisphere.use_tri_nn = MatrixX3d::Zero(p_Hemisphere.nuse_tri,3);
711  p_Hemisphere.use_tri_area = VectorXd::Zero(p_Hemisphere.nuse_tri);
712 
713  for (qint32 i = 0; i < p_Hemisphere.nuse_tri; ++i)
714  {
715  for ( qint32 j = 0; j < 3; ++j)
716  {
717  k = p_Hemisphere.use_tris(i, j);
718 
719  r(j,0) = p_Hemisphere.rr(k, 0);
720  r(j,1) = p_Hemisphere.rr(k, 1);
721  r(j,2) = p_Hemisphere.rr(k, 2);
722 
723  p_Hemisphere.use_tri_cent(i, 0) += p_Hemisphere.rr(k, 0);
724  p_Hemisphere.use_tri_cent(i, 1) += p_Hemisphere.rr(k, 1);
725  p_Hemisphere.use_tri_cent(i, 2) += p_Hemisphere.rr(k, 2);
726  }
727  p_Hemisphere.use_tri_cent.row(i) /= 3.0f;
728 
729  //cross product {cross((r2-r1),(r3-r1))}
730  a = r.row(1) - r.row(0 );
731  b = r.row(2) - r.row(0);
732  p_Hemisphere.use_tri_nn(i,0) = a(1)*b(2)-a(2)*b(1);
733  p_Hemisphere.use_tri_nn(i,1) = a(2)*b(0)-a(0)*b(2);
734  p_Hemisphere.use_tri_nn(i,2) = a(0)*b(1)-a(1)*b(0);
735 
736  //area
737  size = p_Hemisphere.use_tri_nn.row(i)*p_Hemisphere.use_tri_nn.row(i).transpose();
738  size = std::pow(size, 0.5f );
739 
740  p_Hemisphere.use_tri_area(i) = size/2.0f;
741  }
742 
743  }
744  printf("[done]\n");
745 
746 // qDebug() << "p_Hemisphere.use_tri_cent:" << p_Hemisphere.use_tri_cent(0,0) << p_Hemisphere.use_tri_cent(0,1) << p_Hemisphere.use_tri_cent(0,2);
747 // qDebug() << "p_Hemisphere.use_tri_cent:" << p_Hemisphere.use_tri_cent(2,0) << p_Hemisphere.use_tri_cent(2,1) << p_Hemisphere.use_tri_cent(2,2);
748 
749 // qDebug() << "p_Hemisphere.use_tri_nn:" << p_Hemisphere.use_tri_nn(0,0) << p_Hemisphere.use_tri_nn(0,1) << p_Hemisphere.use_tri_nn(0,2);
750 // qDebug() << "p_Hemisphere.use_tri_nn:" << p_Hemisphere.use_tri_nn(2,0) << p_Hemisphere.use_tri_nn(2,1) << p_Hemisphere.use_tri_nn(2,2);
751 
752  printf("\tCompleting triangle and vertex neighboring info...");
753  p_Hemisphere.add_geometry_info();
754  printf("[done]\n");
755 
756  return true;
757 }
758 
759 //=============================================================================================================
760 
762 {
763  for(qint32 h = 0; h < m_qListHemispheres.size(); ++h)
764  {
765  printf("\tWrite a source space... ");
766  p_pStream->start_block(FIFFB_MNE_SOURCE_SPACE);
767  m_qListHemispheres[h].writeToStream(p_pStream);
768  p_pStream->end_block(FIFFB_MNE_SOURCE_SPACE);
769  printf("[done]\n");
770  }
771  printf("\t%d source spaces written\n", m_qListHemispheres.size());
772 }
773 
774 //=============================================================================================================
775 
777 {
778  if(m_qListHemispheres.size() > idx)
779  return m_qListHemispheres[idx];
780  else
781  {
782  qWarning("Warning: Index out of bound! Returning last element.");
783  return m_qListHemispheres[m_qListHemispheres.size()-1];
784  }
785 }
786 
787 //=============================================================================================================
788 
790 {
791  if(m_qListHemispheres.size() > idx)
792  return m_qListHemispheres[idx];
793  else
794  {
795  qWarning("Warning: Index out of bound! Returning last element.");
796  return m_qListHemispheres[m_qListHemispheres.size()-1];
797  }
798 }
799 
800 //=============================================================================================================
801 
803 {
804  if(idt.compare("lh") == 0)
805  return m_qListHemispheres[0];
806  else if(idt.compare("rh") == 0)
807  return m_qListHemispheres[1];
808  else
809  {
810  qWarning("Warning: Identifier is not 'lh' or 'rh'! Returning 'lh'.");
811  return m_qListHemispheres[0];
812  }
813 }
814 
815 //=============================================================================================================
816 
818 {
819  if(idt.compare("lh") == 0)
820  return m_qListHemispheres[0];
821  else if(idt.compare("rh") == 0)
822  return m_qListHemispheres[1];
823  else
824  {
825  qWarning("Warning: Identifier is not 'lh' or 'rh'! Returning 'lh'.");
826  return m_qListHemispheres[0];
827  }
828 }
Eigen::MatrixX3f rr
Eigen::MatrixX3d use_tri_cent
Eigen::VectorXi patch_inds
FIFFLIB::fiff_int_t nuse
FIFFLIB::fiff_int_t coord_frame
Eigen::VectorXi inuse
Eigen::VectorXi nearest
MNEHemisphere & operator[](qint32 idx)
#define FIFF_BEM_SURF_TRIANGLES
Definition: fiff_file.h:739
QList< Eigen::VectorXi > pinfo
FIFFLIB::fiff_int_t np
#define FIFF_MNE_SOURCE_SPACE_USE_TRIANGLES
Eigen::MatrixX3i tris
QSharedPointer< FiffDirNode > SPtr
Definition: fiff_dir_node.h:77
#define FIFF_MNE_SOURCE_SPACE_NUSE
Coordinate transformation description.
Eigen::MatrixX3d use_tri_nn
Eigen::SparseMatrix< double > dist
#define FIFF_MNE_SOURCE_SPACE_SELECTION
bool transform_source_space_to(FIFFLIB::fiff_int_t dest, FIFFLIB::FiffCoordTrans &trans)
FIFFLIB::fiff_int_t ntri
Freesurfer/MNE label.
Definition: label.h:80
#define FIFF_MNE_SOURCE_SPACE_TRIANGLES
Eigen::MatrixX3d tri_cent
MNESourceSpace class declaration.
#define FIFF_MNE_SOURCE_SPACE_NORMALS
static qint32 find_source_space_hemi(MNEHemisphere &p_Hemisphere)
#define FIFF_MNE_SOURCE_SPACE_NEAREST_DIST
Eigen::VectorXi vertno
fiff_long_t start_block(fiff_int_t kind)
#define FIFF_MNE_SOURCE_SPACE_NTRI
#define FIFF_MNE_SOURCE_SPACE_POINTS
QSharedPointer< FiffTag > SPtr
Definition: fiff_tag.h:152
Source Space descritpion.
Hemisphere provides geometry information.
FIFF File I/O routines.
Definition: fiff_stream.h:104
Eigen::VectorXd tri_area
MNEMath class declaration.
static bool readFromStream(FIFFLIB::FiffStream::SPtr &p_pStream, bool add_geom, MNESourceSpace &p_SourceSpace)
#define FIFF_BEM_SURF_NTRI
Definition: fiff_file.h:737
Eigen::VectorXd nearest_dist
Eigen::VectorXd use_tri_area
Eigen::VectorXi vertices
Definition: label.h:166
Eigen::MatrixX3d tri_nn
void writeToStream(FIFFLIB::FiffStream *p_pStream)
qint32 hemi
Definition: label.h:169
FIFFLIB::fiff_int_t id
#define FIFF_MNE_SOURCE_SPACE_DIST
QSharedPointer< FiffStream > SPtr
Definition: fiff_stream.h:107
QList< Eigen::VectorXi > label_src_vertno_sel(const FSLIB::Label &p_label, Eigen::VectorXi &src_sel) const
#define FIFF_MNE_SOURCE_SPACE_DIST_LIMIT
Label class declaration.
QList< Eigen::VectorXi > get_vertno() const
#define FIFF_MNE_SOURCE_SPACE_ID
#define FIFF_MNE_SOURCE_SPACE_NEAREST
MNESourceSpace pick_regions(const QList< FSLIB::Label > &p_qListLabels) const
Eigen::MatrixX3i use_tris
fiff_long_t end_block(fiff_int_t kind, fiff_int_t next=FIFFV_NEXT_SEQ)
static bool patch_info(MNEHemisphere &p_Hemisphere)
#define FIFF_MNE_SOURCE_SPACE_NPOINTS
#define FIFF_MNE_SOURCE_SPACE_NUSE_TRI
Eigen::MatrixX3f nn
#define FIFF_MNE_COORD_FRAME