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