MNE-CPP  beta 1.0
mne_forwardsolution.cpp
Go to the documentation of this file.
1 //=============================================================================================================
36 //*************************************************************************************************************
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "mne_forwardsolution.h"
42 
43 
44 //*************************************************************************************************************
45 //=============================================================================================================
46 // FIFF INCLUDES
47 //=============================================================================================================
48 
49 #include <fs/colortable.h>
50 #include <fs/label.h>
51 #include <utils/mnemath.h>
52 #include <utils/kmeans.h>
53 
54 
55 //*************************************************************************************************************
56 //=============================================================================================================
57 // STL INCLUDES
58 //=============================================================================================================
59 
60 #include <iostream>
61 #include <QtConcurrent>
62 #include <QFuture>
63 
64 
65 //*************************************************************************************************************
66 //=============================================================================================================
67 // USED NAMESPACES
68 //=============================================================================================================
69 
70 using namespace MNELIB;
71 using namespace UTILSLIB;
72 using namespace FSLIB;
73 
74 
75 //*************************************************************************************************************
76 //=============================================================================================================
77 // DEFINE MEMBER METHODS
78 //=============================================================================================================
79 
81 : source_ori(-1)
82 , surf_ori(false)
83 , coord_frame(-1)
84 , nsource(-1)
85 , nchan(-1)
86 , sol(new FiffNamedMatrix)
87 , sol_grad(new FiffNamedMatrix)
88 //, mri_head_t(NULL)
89 //, src(NULL)
90 , source_rr(MatrixX3f::Zero(0,3))
91 , source_nn(MatrixX3f::Zero(0,3))
92 {
93 
94 }
95 
96 
97 //*************************************************************************************************************
98 
99 MNEForwardSolution::MNEForwardSolution(QIODevice &p_IODevice, bool force_fixed, bool surf_ori, const QStringList& include, const QStringList& exclude, bool bExcludeBads)
100 : source_ori(-1)
101 , surf_ori(surf_ori)
102 , coord_frame(-1)
103 , nsource(-1)
104 , nchan(-1)
105 , sol(new FiffNamedMatrix)
106 , sol_grad(new FiffNamedMatrix)
107 //, mri_head_t(NULL)
108 //, src(NULL)
109 , source_rr(MatrixX3f::Zero(0,3))
110 , source_nn(MatrixX3f::Zero(0,3))
111 {
112  if(!read(p_IODevice, *this, force_fixed, surf_ori, include, exclude, bExcludeBads))
113  {
114  printf("\tForward solution not found.\n");//ToDo Throw here
115  return;
116  }
117 }
118 
119 
120 //*************************************************************************************************************
121 
123 : info(p_MNEForwardSolution.info)
124 , source_ori(p_MNEForwardSolution.source_ori)
125 , surf_ori(p_MNEForwardSolution.surf_ori)
126 , coord_frame(p_MNEForwardSolution.coord_frame)
127 , nsource(p_MNEForwardSolution.nsource)
128 , nchan(p_MNEForwardSolution.nchan)
129 , sol(p_MNEForwardSolution.sol)
130 , sol_grad(p_MNEForwardSolution.sol_grad)
131 , mri_head_t(p_MNEForwardSolution.mri_head_t)
132 , src(p_MNEForwardSolution.src)
133 , source_rr(p_MNEForwardSolution.source_rr)
134 , source_nn(p_MNEForwardSolution.source_nn)
135 {
136 
137 }
138 
139 
140 //*************************************************************************************************************
141 
143 {
144 
145 }
146 
147 
148 //*************************************************************************************************************
149 
151 {
152  info.clear();
153  source_ori = -1;
154  surf_ori = false;
155  coord_frame = -1;
156  nsource = -1;
157  nchan = -1;
160  mri_head_t.clear();
161  src.clear();
162  source_rr = MatrixX3f(0,3);
163  source_nn = MatrixX3f(0,3);
164 }
165 
166 
167 //*************************************************************************************************************
168 
169 MNEForwardSolution MNEForwardSolution::cluster_forward_solution(const AnnotationSet &p_AnnotationSet, qint32 p_iClusterSize, MatrixXd& p_D, const FiffCov &p_pNoise_cov, const FiffInfo &p_pInfo, QString p_sMethod) const
170 {
171  MNEForwardSolution p_fwdOut = MNEForwardSolution(*this);
172 
173  printf("Cluster forward solution using %s.\n", p_sMethod.toUtf8().constData());
174 
175 // qDebug() << "this->sol->data" << this->sol->data.rows() << "x" << this->sol->data.cols();
176 
177  //
178  // Check consisty
179  //
180  if(this->isFixedOrient())
181  {
182  printf("Error: Fixed orientation not implemented jet!\n");
183  return p_fwdOut;
184  }
185 
186 // for(qint32 h = 0; h < this->src.hemispheres.size(); ++h )//obj.sizeForwardSolution)
187 // {
188 // if(this->src[h]->vertno.rows() != t_listAnnotation[h]->getLabel()->rows())
189 // {
190 // printf("Error: Annotation doesn't fit to Forward Solution: Vertice number is different!");
191 // return false;
192 // }
193 // }
194 
195  MatrixXd t_G_Whitened(0,0);
196  bool t_bUseWhitened = false;
197  //
198  //Whiten gain matrix before clustering -> cause diffenerent units Magnetometer, Gradiometer and EEG
199  //
200  if(!p_pNoise_cov.isEmpty() && !p_pInfo.isEmpty())
201  {
202  FiffInfo p_outFwdInfo;
203  FiffCov p_outNoiseCov;
204  MatrixXd p_outWhitener;
205  qint32 p_outNumNonZero;
206  //do whitening with noise cov
207  this->prepare_forward(p_pInfo, p_pNoise_cov, false, p_outFwdInfo, t_G_Whitened, p_outNoiseCov, p_outWhitener, p_outNumNonZero);
208  printf("\tWhitening the forward solution.\n");
209 
210  t_G_Whitened = p_outWhitener*t_G_Whitened;
211  t_bUseWhitened = true;
212  }
213 
214 
215  //
216  // Assemble input data
217  //
218  qint32 count;
219  qint32 offset;
220 
221  MatrixXd t_G_new;
222 
223  for(qint32 h = 0; h < this->src.size(); ++h )
224  {
225 
226  count = 0;
227  offset = 0;
228 
229  // Offset for continuous indexing;
230  if(h > 0)
231  for(qint32 j = 0; j < h; ++j)
232  offset += this->src[j].nuse;
233 
234  if(h == 0)
235  printf("Cluster Left Hemisphere\n");
236  else
237  printf("Cluster Right Hemisphere\n");
238 
239  Colortable t_CurrentColorTable = p_AnnotationSet[h].getColortable();
240  VectorXi label_ids = t_CurrentColorTable.getLabelIds();
241 
242  // Get label ids for every vertex
243  VectorXi vertno_labeled = VectorXi::Zero(this->src[h].vertno.rows());
244 
245  //ToDo make this more universal -> using Label instead of annotations - obsolete when using Labels
246  for(qint32 i = 0; i < vertno_labeled.rows(); ++i)
247  vertno_labeled[i] = p_AnnotationSet[h].getLabelIds()[this->src[h].vertno[i]];
248 
249  //Qt Concurrent List
250  QList<RegionData> m_qListRegionDataIn;
251 
252  //
253  // Generate cluster input data
254  //
255  for (qint32 i = 0; i < label_ids.rows(); ++i)
256  {
257  if (label_ids[i] != 0)
258  {
259  QString curr_name = t_CurrentColorTable.struct_names[i];//obj.label2AtlasName(label(i));
260  printf("\tCluster %d / %li %s...", i+1, label_ids.rows(), curr_name.toUtf8().constData());
261 
262  //
263  // Get source space indeces
264  //
265  VectorXi idcs = VectorXi::Zero(vertno_labeled.rows());
266  qint32 c = 0;
267 
268  //Select ROIs //change this use label info with a hash tabel
269  for(qint32 j = 0; j < vertno_labeled.rows(); ++j)
270  {
271  if(vertno_labeled[j] == label_ids[i])
272  {
273  idcs[c] = j;
274  ++c;
275  }
276  }
277  idcs.conservativeResize(c);
278 
279  //get selected G
280  MatrixXd t_G(this->sol->data.rows(), idcs.rows()*3);
281  MatrixXd t_G_Whitened_Roi(t_G_Whitened.rows(), idcs.rows()*3);
282 
283  for(qint32 j = 0; j < idcs.rows(); ++j)
284  {
285  t_G.block(0, j*3, t_G.rows(), 3) = this->sol->data.block(0, (idcs[j]+offset)*3, t_G.rows(), 3);
286  if(t_bUseWhitened)
287  t_G_Whitened_Roi.block(0, j*3, t_G_Whitened_Roi.rows(), 3) = t_G_Whitened.block(0, (idcs[j]+offset)*3, t_G_Whitened_Roi.rows(), 3);
288  }
289 
290  qint32 nSens = t_G.rows();
291  qint32 nSources = t_G.cols()/3;
292 
293  if (nSources > 0)
294  {
295  RegionData t_sensG;
296 
297  t_sensG.idcs = idcs;
298  t_sensG.iLabelIdxIn = i;
299  t_sensG.nClusters = ceil((double)nSources/(double)p_iClusterSize);
300 
301  t_sensG.matRoiGOrig = t_G;
302 // if(t_bUseWhitened)
303 // t_sensG.matRoiGOrigWhitened = t_G_Whitened_Roi;
304 
305  printf("%d Cluster(s)... ", t_sensG.nClusters);
306 
307  // Reshape Input data -> sources rows; sensors columns
308  t_sensG.matRoiG = MatrixXd(t_G.cols()/3, 3*nSens);
309  if(t_bUseWhitened)
310  t_sensG.matRoiGWhitened = MatrixXd(t_G_Whitened_Roi.cols()/3, 3*nSens);
311 
312  for(qint32 j = 0; j < nSens; ++j)
313  {
314  for(qint32 k = 0; k < t_sensG.matRoiG.rows(); ++k)
315  t_sensG.matRoiG.block(k,j*3,1,3) = t_G.block(j,k*3,1,3);
316  if(t_bUseWhitened)
317  for(qint32 k = 0; k < t_sensG.matRoiGWhitened.rows(); ++k)
318  t_sensG.matRoiGWhitened.block(k,j*3,1,3) = t_G_Whitened_Roi.block(j,k*3,1,3);
319  }
320 
321  t_sensG.bUseWhitened = t_bUseWhitened;
322 
323  t_sensG.sDistMeasure = p_sMethod;
324 
325  m_qListRegionDataIn.append(t_sensG);
326 
327  printf("[added]\n");
328  }
329  else
330  {
331  printf("failed! Label contains no sources.\n");
332  }
333  }
334  }
335 
336 
337  //
338  // Calculate clusters
339  //
340  printf("Clustering... ");
341  QFuture< RegionDataOut > res;
342  res = QtConcurrent::mapped(m_qListRegionDataIn, &RegionData::cluster);
343  res.waitForFinished();
344 
345  //
346  // Assign results
347  //
348  MatrixXd t_G_partial;
349 
350  qint32 nClusters;
351  qint32 nSens;
352  QList<RegionData>::const_iterator itIn;
353  itIn = m_qListRegionDataIn.begin();
354  QFuture<RegionDataOut>::const_iterator itOut;
355  for (itOut = res.constBegin(); itOut != res.constEnd(); ++itOut)
356  {
357  nClusters = itOut->ctrs.rows();
358  nSens = itOut->ctrs.cols()/3;
359  t_G_partial = MatrixXd::Zero(nSens, nClusters*3);
360 
361 // std::cout << "Number of Clusters: " << nClusters << " x " << nSens << std::endl;//itOut->iLabelIdcsOut << std::endl;
362 
363  //
364  // Assign the centroid for each cluster to the partial G
365  //
366  //ToDo change this use indeces found with whitened data
367  for(qint32 j = 0; j < nSens; ++j)
368  for(qint32 k = 0; k < nClusters; ++k)
369  t_G_partial.block(j, k*3, 1, 3) = itOut->ctrs.block(k,j*3,1,3);
370 
371  //
372  // Get cluster indizes and its distances to the centroid
373  //
374  for(qint32 j = 0; j < nClusters; ++j)
375  {
376  VectorXi clusterIdcs = VectorXi::Zero(itOut->roiIdx.rows());
377  VectorXd clusterDistance = VectorXd::Zero(itOut->roiIdx.rows());
378  MatrixX3f clusterSource_rr = MatrixX3f::Zero(itOut->roiIdx.rows(), 3);
379  qint32 nClusterIdcs = 0;
380  for(qint32 k = 0; k < itOut->roiIdx.rows(); ++k)
381  {
382  if(itOut->roiIdx[k] == j)
383  {
384  clusterIdcs[nClusterIdcs] = itIn->idcs[k];
385 
386  qint32 offset = h == 0 ? 0 : this->src[0].nuse;
387  clusterSource_rr.row(nClusterIdcs) = this->source_rr.row(offset + itIn->idcs[k]);
388  clusterDistance[nClusterIdcs] = itOut->D(k,j);
389  ++nClusterIdcs;
390  }
391  }
392  clusterIdcs.conservativeResize(nClusterIdcs);
393  clusterSource_rr.conservativeResize(nClusterIdcs,3);
394  clusterDistance.conservativeResize(nClusterIdcs);
395 
396  VectorXi clusterVertnos = VectorXi::Zero(clusterIdcs.size());
397  for(qint32 k = 0; k < clusterVertnos.size(); ++k)
398  clusterVertnos(k) = this->src[h].vertno[clusterIdcs(k)];
399 
400 
401  p_fwdOut.src[h].cluster_info.clusterVertnos.append(clusterVertnos);
402  p_fwdOut.src[h].cluster_info.clusterSource_rr.append(clusterSource_rr);
403  p_fwdOut.src[h].cluster_info.clusterDistances.append(clusterDistance);
404  p_fwdOut.src[h].cluster_info.clusterLabelIds.append(label_ids[itOut->iLabelIdxOut]);
405  p_fwdOut.src[h].cluster_info.clusterLabelNames.append(t_CurrentColorTable.getNames()[itOut->iLabelIdxOut]);
406  }
407 
408 
409  //
410  // Assign partial G to new LeadField
411  //
412  if(t_G_partial.rows() > 0 && t_G_partial.cols() > 0)
413  {
414  t_G_new.conservativeResize(t_G_partial.rows(), t_G_new.cols() + t_G_partial.cols());
415  t_G_new.block(0, t_G_new.cols() - t_G_partial.cols(), t_G_new.rows(), t_G_partial.cols()) = t_G_partial;
416 
417  // Map the centroids to the closest rr
418  for(qint32 k = 0; k < nClusters; ++k)
419  {
420  qint32 j = 0;
421 
422  double sqec = sqrt((itIn->matRoiGOrig.block(0, j*3, itIn->matRoiGOrig.rows(), 3) - t_G_partial.block(0, k*3, t_G_partial.rows(), 3)).array().pow(2).sum());
423  double sqec_min = sqec;
424  qint32 j_min = 0;
425 // MatrixXd matGainDiff;
426  for(qint32 j = 1; j < itIn->idcs.rows(); ++j)
427  {
428  sqec = sqrt((itIn->matRoiGOrig.block(0, j*3, itIn->matRoiGOrig.rows(), 3) - t_G_partial.block(0, k*3, t_G_partial.rows(), 3)).array().pow(2).sum());
429 
430  if(sqec < sqec_min)
431  {
432  sqec_min = sqec;
433  j_min = j;
434 // matGainDiff = itIn->matRoiGOrig.block(0, j*3, itIn->matRoiGOrig.rows(), 3) - t_G_partial.block(0, k*3, t_G_partial.rows(), 3);
435  }
436  }
437 
438 // qListGainDist.append(matGainDiff);
439 
440  // Take the closest coordinates
441  qint32 sel_idx = itIn->idcs[j_min];
442 
443  p_fwdOut.src[h].cluster_info.centroidVertno.append(this->src[h].vertno[sel_idx]);
444  p_fwdOut.src[h].cluster_info.centroidSource_rr.append(this->src[h].rr.row(sel_idx));
445 // p_fwdOut.src[h].nn.row(count) = MatrixXd::Zero(1,3);
446 
447 // // Option 1 closest vertno
448 // p_fwdOut.src[h].vertno[count] = this->src[h].vertno[sel_idx]; //ToDo resizing necessary?
449  // Option 2 label ID
450  p_fwdOut.src[h].vertno[count] = p_fwdOut.src[h].cluster_info.clusterLabelIds[count];
451 
452 
453 
454 // //vertices
455 // std::cout << this->src[h].vertno[sel_idx] << ", ";
456 
457  ++count;
458  }
459  }
460 
461  ++itIn;
462  }
463 
464  //
465  // Assemble new hemisphere information
466  //
467 // p_fwdOut.src[h].rr.conservativeResize(count, 3);
468 // p_fwdOut.src[h].nn.conservativeResize(count, 3);
469  p_fwdOut.src[h].vertno.conservativeResize(count);
470 
471  printf("[done]\n");
472  }
473 
474 
475  //
476  // Cluster operator D (sources x clusters)
477  //
478  qint32 totalNumOfClust = 0;
479  for (qint32 h = 0; h < 2; ++h)
480  totalNumOfClust += p_fwdOut.src[h].cluster_info.clusterVertnos.size();
481 
482  if(this->isFixedOrient())
483  p_D = MatrixXd::Zero(this->sol->data.cols(), totalNumOfClust);
484  else
485  p_D = MatrixXd::Zero(this->sol->data.cols(), totalNumOfClust*3);
486 
487  QList<VectorXi> t_vertnos = this->src.get_vertno();
488 
489 // qDebug() << "Size: " << t_vertnos[0].size() << t_vertnos[1].size();
490 // qDebug() << "this->sol->data.cols(): " << this->sol->data.cols();
491 
492  qint32 currentCluster = 0;
493  for (qint32 h = 0; h < 2; ++h)
494  {
495  int hemiOffset = h == 0 ? 0 : t_vertnos[0].size();
496  for(qint32 i = 0; i < p_fwdOut.src[h].cluster_info.clusterVertnos.size(); ++i)
497  {
498  VectorXi idx_sel;
499  MNEMath::intersect(t_vertnos[h], p_fwdOut.src[h].cluster_info.clusterVertnos[i], idx_sel);
500 
501 // std::cout << "\nVertnos:\n" << t_vertnos[h] << std::endl;
502 
503 // std::cout << "clusterVertnos[i]:\n" << p_fwdOut.src[h].cluster_info.clusterVertnos[i] << std::endl;
504 
505  idx_sel.array() += hemiOffset;
506 
507 // std::cout << "idx_sel]:\n" << idx_sel << std::endl;
508 
509 
510 
511  double selectWeight = 1.0/idx_sel.size();
512  if(this->isFixedOrient())
513  {
514  for(qint32 j = 0; j < idx_sel.size(); ++j)
515  p_D.col(currentCluster)[idx_sel(j)] = selectWeight;
516  }
517  else
518  {
519  qint32 clustOffset = currentCluster*3;
520  for(qint32 j = 0; j < idx_sel.size(); ++j)
521  {
522  qint32 idx_sel_Offset = idx_sel(j)*3;
523  //x
524  p_D(idx_sel_Offset,clustOffset) = selectWeight;
525  //y
526  p_D(idx_sel_Offset+1, clustOffset+1) = selectWeight;
527  //z
528  p_D(idx_sel_Offset+2, clustOffset+2) = selectWeight;
529  }
530  }
531  ++currentCluster;
532  }
533  }
534 
535 // std::cout << "D:\n" << D.row(0) << std::endl << D.row(1) << std::endl << D.row(2) << std::endl << D.row(3) << std::endl << D.row(4) << std::endl << D.row(5) << std::endl;
536 
537 
538 // //
539 // // get mean and std of original gain matrix
540 // //
541 // double mean = 0;
542 // qint32 c = 0;
543 // for(qint32 i = 0; i < this->sol->data.rows(); ++i)
544 // {
545 // if(i % 3 == 1 || i%3 == 2)
546 // {
547 // for(qint32 j = 0; j < this->sol->data.cols(); ++j)
548 // {
549 // mean += this->sol->data(i,j);
550 // ++c;
551 // }
552 // }
553 // }
554 
555 // mean /= c;
556 // double var = 0;
557 // c = 0;
558 // for(qint32 i = 0; i < this->sol->data.rows(); ++i)
559 // {
560 // if(i % 3 == 1 || i%3 == 2)
561 // {
562 // for(qint32 j = 0; j < this->sol->data.cols(); ++j)
563 // {
564 // var += pow(this->sol->data(i,j) - mean,2);
565 // ++c;
566 // }
567 // }
568 // }
569 // var /= c - 1;
570 // var = sqrt(var);
571 // std::cout << "Original gain matrix (gradiometer):\n mean: " << mean << "\n var: " << var << std::endl;
572 
573 // mean = 0;
574 // c = 0;
575 // for(qint32 i = 0; i < this->sol->data.rows(); ++i)
576 // {
577 // if(i % 3 == 0)
578 // {
579 // for(qint32 j = 0; j < this->sol->data.cols(); ++j)
580 // {
581 // mean += this->sol->data(i,j);
582 // ++c;
583 // }
584 // }
585 // }
586 // mean /= c;
587 
588 // var = 0;
589 // c = 0;
590 // for(qint32 i = 0; i < this->sol->data.rows(); ++i)
591 // {
592 // if(i % 3 == 0)
593 // {
594 // for(qint32 j = 0; j < this->sol->data.cols(); ++j)
595 // {
596 // var += pow(this->sol->data(i,j) - mean,2);
597 // ++c;
598 // }
599 // }
600 // }
601 // var /= c - 1;
602 // var = sqrt(var);
603 // std::cout << "Original gain matrix (magnetometer):\n mean: " << mean << "\n var: " << var << std::endl;
604 
605 // //
606 // // get mean and std of original gain matrix mapping
607 // //
608 // mean = 0;
609 // c = 0;
610 // for(qint32 h = 0; h < qListGainDist.size(); ++h)
611 // {
612 // for(qint32 i = 0; i < qListGainDist[h].rows(); ++i)
613 // {
614 // if(i % 3 == 1 || i%3 == 2)
615 // {
616 // for(qint32 j = 0; j < qListGainDist[h].cols(); ++j)
617 // {
618 // mean += qListGainDist[h](i,j);
619 // ++c;
620 // }
621 // }
622 // }
623 // }
624 // mean /= c;
625 
626 // var = 0;
627 // c = 0;
628 // for(qint32 h = 0; h < qListGainDist.size(); ++h)
629 // {
630 // for(qint32 i = 0; i < qListGainDist[h].rows(); ++i)
631 // {
632 // if(i % 3 == 1 || i%3 == 2)
633 // {
634 // for(qint32 j = 0; j < qListGainDist[h].cols(); ++j)
635 // {
636 // var += pow(qListGainDist[i](i,j) - mean,2);
637 // ++c;
638 // }
639 // }
640 // }
641 // }
642 
643 // var /= c - 1;
644 // var = sqrt(var);
645 
646 // std::cout << "Gain matrix offset mapping (gradiometer):\n mean: " << mean << "\n var: " << var << std::endl;
647 
648 // mean = 0;
649 // c = 0;
650 // for(qint32 h = 0; h < qListGainDist.size(); ++h)
651 // {
652 // for(qint32 i = 0; i < qListGainDist[h].rows(); ++i)
653 // {
654 // if(i % 3 == 0)
655 // {
656 // for(qint32 j = 0; j < qListGainDist[h].cols(); ++j)
657 // {
658 // mean += qListGainDist[h](i,j);
659 // ++c;
660 // }
661 // }
662 // }
663 // }
664 // mean /= c;
665 
666 // var = 0;
667 // c = 0;
668 // for(qint32 h = 0; h < qListGainDist.size(); ++h)
669 // {
670 // for(qint32 i = 0; i < qListGainDist[h].rows(); ++i)
671 // {
672 // if(i % 3 == 0)
673 // {
674 // for(qint32 j = 0; j < qListGainDist[h].cols(); ++j)
675 // {
676 // var += pow(qListGainDist[i](i,j) - mean,2);
677 // ++c;
678 // }
679 // }
680 // }
681 // }
682 
683 // var /= c - 1;
684 // var = sqrt(var);
685 
686 // std::cout << "Gain matrix offset mapping (magnetometer):\n mean: " << mean << "\n var: " << var << std::endl;
687 
688 
689  //
690  // Put it all together
691  //
692  p_fwdOut.sol->data = t_G_new;
693  p_fwdOut.sol->ncol = t_G_new.cols();
694 
695  p_fwdOut.nsource = p_fwdOut.sol->ncol/3;
696 
697  return p_fwdOut;
698 }
699 
700 
701 //*************************************************************************************************************
702 
703 MNEForwardSolution MNEForwardSolution::reduce_forward_solution(qint32 p_iNumDipoles, MatrixXd& p_D) const
704 {
705  MNEForwardSolution p_fwdOut = MNEForwardSolution(*this);
706 
707  bool isFixed = p_fwdOut.isFixedOrient();
708  qint32 np = isFixed ? p_fwdOut.sol->data.cols() : p_fwdOut.sol->data.cols()/3;
709 
710  if(p_iNumDipoles > np)
711  return p_fwdOut;
712 
713  VectorXi sel(p_iNumDipoles);
714 
715  float t_fStep = (float)np/(float)p_iNumDipoles;
716 
717  for(qint32 i = 0; i < p_iNumDipoles; ++i)
718  {
719  float t_fCurrent = ((float)i)*t_fStep;
720  sel[i] = (quint32)floor(t_fCurrent);
721  }
722 
723  if(isFixed)
724  {
725  p_D = MatrixXd::Zero(p_fwdOut.sol->data.cols(), p_iNumDipoles);
726  for(qint32 i = 0; i < p_iNumDipoles; ++i)
727  p_D(sel[i], i) = 1;
728  }
729  else
730  {
731  p_D = MatrixXd::Zero(p_fwdOut.sol->data.cols(), p_iNumDipoles*3);
732  for(qint32 i = 0; i < p_iNumDipoles; ++i)
733  for(qint32 j = 0; j < 3; ++j)
734  p_D((sel[i]*3)+j, (i*3)+j) = 1;
735  }
736 
737 
738 // //find idx of hemi switch
739 // qint32 vertno_size = this->src[0].nuse;
740 // qint32 hIdx = 0;
741 
742 // //LH
743 // VectorXi vertnosLH(this->src[0].nuse);
744 // for(qint32 i = 0; i < sel.size(); ++i)
745 // {
746 // if(sel[i] >= vertno_size)
747 // {
748 // hIdx = i;
749 // break;
750 // }
751 // vertnosLH[i] = this->src[0].vertno(sel[i]);
752 // }
753 // vertnosLH.conservativeResize(hIdx);
754 
755 // QFile file_centroids_LH("./centroids_LH_sel.txt");
756 // file_centroids_LH.open(QIODevice::WriteOnly | QIODevice::Text);
757 // QTextStream out_centroids_LH(&file_centroids_LH);
758 // for(qint32 i = 0; i < vertnosLH.size(); ++i)
759 // out_centroids_LH << vertnosLH[i] << ", ";
760 // file_centroids_LH.close();
761 
762 // QFile file_centroids_LH_full("./centroids_LH_full.txt");
763 // file_centroids_LH_full.open(QIODevice::WriteOnly | QIODevice::Text);
764 // QTextStream out_centroids_LH_full(&file_centroids_LH_full);
765 // for(qint32 i = 0; i < this->src[0].vertno.size(); ++i)
766 // out_centroids_LH_full << this->src[0].vertno[i] << ", ";
767 // file_centroids_LH_full.close();
768 
769 // //RH
770 // VectorXi vertnosRH(sel.size() - hIdx);
771 // for(qint32 i = hIdx; i < sel.size(); ++i)
772 // {
773 // vertnosRH[i - hIdx] = this->src[1].vertno(sel[i] - hIdx);
774 // }
775 
776 // QFile file_centroids_RH("./centroids_RH_sel.txt");
777 // file_centroids_RH.open(QIODevice::WriteOnly | QIODevice::Text);
778 // QTextStream out_centroids_RH(&file_centroids_RH);
779 // for(qint32 i = 0; i < vertnosRH.size(); ++i)
780 // out_centroids_RH << vertnosRH[i] << ", ";
781 // file_centroids_RH.close();
782 // //vertno end
783 
784  // New gain matrix
785  p_fwdOut.sol->data = this->sol->data * p_D;
786 
787  MatrixX3f rr(p_iNumDipoles,3);
788 
789  MatrixX3f nn(p_iNumDipoles,3);
790 
791 
792  for(qint32 i = 0; i < p_iNumDipoles; ++i)
793  {
794  rr.row(i) = this->source_rr.row(sel(i));
795  nn.row(i) = this->source_nn.row(sel(i));
796  }
797 
798  p_fwdOut.source_rr = rr;
799  p_fwdOut.source_nn = nn;
800 
801  p_fwdOut.sol->ncol = p_fwdOut.sol->data.cols();
802 
803  p_fwdOut.nsource = p_iNumDipoles;
804 
805  return p_fwdOut;
806 }
807 
808 
809 //*************************************************************************************************************
810 
811 FiffCov MNEForwardSolution::compute_depth_prior(const MatrixXd &Gain, const FiffInfo &gain_info, bool is_fixed_ori, double exp, double limit, const MatrixXd &patch_areas, bool limit_depth_chs)
812 {
813  printf("\tCreating the depth weighting matrix...\n");
814 
815  MatrixXd G(Gain);
816  // If possible, pick best depth-weighting channels
817  if(limit_depth_chs)
819 
820  VectorXd d;
821  // Compute the gain matrix
822  if(is_fixed_ori)
823  {
824  d = (G.array().square()).rowwise().sum(); //ToDo: is this correct - is G squared?
825 // d = np.sum(G ** 2, axis=0)
826  }
827  else
828  {
829  qint32 n_pos = G.cols() / 3;
830  d = VectorXd::Zero(n_pos);
831  MatrixXd Gk;
832  for (qint32 k = 0; k < n_pos; ++k)
833  {
834  Gk = G.block(0,3*k, G.rows(), 3);
835  JacobiSVD<MatrixXd> svd(Gk.transpose()*Gk);
836  d[k] = svd.singularValues().maxCoeff();
837  }
838  }
839 
840  // ToDo Currently the fwd solns never have "patch_areas" defined
841  if(patch_areas.cols() > 0)
842  {
843 // d /= patch_areas ** 2
844  printf("\tToDo!!!!! >>> Patch areas taken into account in the depth weighting\n");
845  }
846 
847  qint32 n_limit;
848  VectorXd w = d.cwiseInverse();
849  VectorXd ws = w;
850  VectorXd wpp;
851  MNEMath::sort<double>(ws, false);
852  double weight_limit = pow(limit, 2);
853  if (!limit_depth_chs)
854  {
855  // match old mne-python behavor
856  qint32 ind;
857  ws.minCoeff(&ind);
858  n_limit = ind;
859  limit = ws[ind] * weight_limit;
860  }
861  else
862  {
863  // match C code behavior
864  limit = ws[ws.size()-1];
865  qint32 ind;
866  n_limit = d.size();
867  if (ws[ws.size()-1] > weight_limit * ws[0])
868  {
869  double th = weight_limit * ws[0];
870  for(qint32 i = 0; i < ws.size(); ++i)
871  {
872  if(ws[i] > th)
873  {
874  ind = i;
875  break;
876  }
877  }
878  limit = ws[ind];
879  n_limit = ind;
880  }
881  }
882 
883  printf("\tlimit = %d/%li = %f", n_limit + 1, d.size(), sqrt(limit / ws[0]));
884  double scale = 1.0 / limit;
885  printf("\tscale = %g exp = %g", scale, exp);
886 
887  VectorXd t_w = w.array() / limit;
888  for(qint32 i = 0; i < t_w.size(); ++i)
889  t_w[i] = t_w[i] > 1 ? 1 : t_w[i];
890  wpp = t_w.array().pow(exp);
891 
892  FiffCov depth_prior;
893  if(is_fixed_ori)
894  depth_prior.data = wpp;
895  else
896  {
897  depth_prior.data.resize(wpp.rows()*3, 1);
898  qint32 idx = 0;
899  double v;
900  for(qint32 i = 0; i < wpp.rows(); ++i)
901  {
902  idx = i*3;
903  v = wpp[i];
904  depth_prior.data(idx, 0) = v;
905  depth_prior.data(idx+1, 0) = v;
906  depth_prior.data(idx+2, 0) = v;
907  }
908  }
909 
910  depth_prior.kind = FIFFV_MNE_DEPTH_PRIOR_COV;
911  depth_prior.diag = true;
912  depth_prior.dim = depth_prior.data.rows();
913  depth_prior.nfree = 1;
914 
915  return depth_prior;
916 }
917 
918 
919 //*************************************************************************************************************
920 
922 {
923  bool is_fixed_ori = this->isFixedOrient();
924  qint32 n_sources = this->sol->data.cols();
925 
926  if (0 <= loose && loose <= 1)
927  {
928  if(loose < 1 && !this->surf_ori)
929  {
930  printf("\tForward operator is not oriented in surface coordinates. loose parameter should be None not %f.", loose);//ToDo Throw here
931  loose = 1;
932  printf("\tSetting loose to %f.\n", loose);
933  }
934 
935  if(is_fixed_ori)
936  {
937  printf("\tIgnoring loose parameter with forward operator with fixed orientation.\n");
938  loose = 0.0;
939  }
940  }
941  else
942  {
943  if(loose < 0 || loose > 1)
944  {
945  qWarning("Warning: Loose value should be in interval [0,1] not %f.\n", loose);
946  loose = loose > 1 ? 1 : 0;
947  printf("Setting loose to %f.\n", loose);
948  }
949  }
950 
951  FiffCov orient_prior;
952  orient_prior.data = VectorXd::Ones(n_sources);
953  if(!is_fixed_ori && (0 <= loose && loose <= 1))
954  {
955  printf("\tApplying loose dipole orientations. Loose value of %f.\n", loose);
956  for(qint32 i = 0; i < n_sources; i+=3)
957  orient_prior.data.block(i,0,2,1).array() *= loose;
958 
959  orient_prior.kind = FIFFV_MNE_ORIENT_PRIOR_COV;
960  orient_prior.diag = true;
961  orient_prior.dim = orient_prior.data.size();
962  orient_prior.nfree = 1;
963  }
964  return orient_prior;
965 }
966 
967 
968 //*************************************************************************************************************
969 
970 MNEForwardSolution MNEForwardSolution::pick_channels(const QStringList& include, const QStringList& exclude) const
971 {
972  MNEForwardSolution fwd(*this);
973 
974  if(include.size() == 0 && exclude.size() == 0)
975  return fwd;
976 
977  RowVectorXi sel = FiffInfo::pick_channels(fwd.sol->row_names, include, exclude);
978 
979  // Do we have something?
980  quint32 nuse = sel.size();
981 
982  if (nuse == 0)
983  {
984  printf("Nothing remains after picking. Returning original forward solution.\n");
985  return fwd;
986  }
987  printf("\t%d out of %d channels remain after picking\n", nuse, fwd.nchan);
988 
989  // Pick the correct rows of the forward operator
990  MatrixXd newData(nuse, fwd.sol->data.cols());
991  for(quint32 i = 0; i < nuse; ++i)
992  newData.row(i) = fwd.sol->data.row(sel[i]);
993 
994  fwd.sol->data = newData;
995  fwd.sol->nrow = nuse;
996 
997  QStringList ch_names;
998  for(qint32 i = 0; i < sel.cols(); ++i)
999  ch_names << fwd.sol->row_names[sel(i)];
1000  fwd.nchan = nuse;
1001  fwd.sol->row_names = ch_names;
1002 
1003  QList<FiffChInfo> chs;
1004  for(qint32 i = 0; i < sel.cols(); ++i)
1005  chs.append(fwd.info.chs[sel(i)]);
1006  fwd.info.chs = chs;
1007  fwd.info.nchan = nuse;
1008 
1009  QStringList bads;
1010  for(qint32 i = 0; i < fwd.info.bads.size(); ++i)
1011  if(ch_names.contains(fwd.info.bads[i]))
1012  bads.append(fwd.info.bads[i]);
1013  fwd.info.bads = bads;
1014 
1015  if(!fwd.sol_grad->isEmpty())
1016  {
1017  newData.resize(nuse, fwd.sol_grad->data.cols());
1018  for(quint32 i = 0; i < nuse; ++i)
1019  newData.row(i) = fwd.sol_grad->data.row(sel[i]);
1020  fwd.sol_grad->data = newData;
1021  fwd.sol_grad->nrow = nuse;
1022  QStringList row_names;
1023  for(qint32 i = 0; i < sel.cols(); ++i)
1024  row_names << fwd.sol_grad->row_names[sel(i)];
1025  fwd.sol_grad->row_names = row_names;
1026  }
1027 
1028  return fwd;
1029 }
1030 
1031 
1032 //*************************************************************************************************************
1033 
1034 MNEForwardSolution MNEForwardSolution::pick_regions(const QList<Label> &p_qListLabels) const
1035 {
1036  VectorXi selVertices;
1037 
1038  qint32 iSize = 0;
1039  for(qint32 i = 0; i < p_qListLabels.size(); ++i)
1040  {
1041  VectorXi currentSelection;
1042  this->src.label_src_vertno_sel(p_qListLabels[i], currentSelection);
1043 
1044  selVertices.conservativeResize(iSize+currentSelection.size());
1045  selVertices.block(iSize,0,currentSelection.size(),1) = currentSelection;
1046  iSize = selVertices.size();
1047  }
1048 
1049  MNEMath::sort(selVertices, false);
1050 
1051  MNEForwardSolution selectedFwd(*this);
1052 
1053  MatrixX3f rr(selVertices.size(),3);
1054  MatrixX3f nn(selVertices.size(),3);
1055 
1056  for(qint32 i = 0; i < selVertices.size(); ++i)
1057  {
1058  rr.block(i, 0, 1, 3) = selectedFwd.source_rr.row(selVertices[i]);
1059  nn.block(i, 0, 1, 3) = selectedFwd.source_nn.row(selVertices[i]);
1060  }
1061 
1062  selectedFwd.source_rr = rr;
1063  selectedFwd.source_nn = nn;
1064 
1065  VectorXi selSolIdcs = tripletSelection(selVertices);
1066  MatrixXd G(selectedFwd.sol->data.rows(),selSolIdcs.size());
1067 // selectedFwd.sol_grad; //ToDo
1068  qint32 rows = G.rows();
1069 
1070  for(qint32 i = 0; i < selSolIdcs.size(); ++i)
1071  G.block(0, i, rows, 1) = selectedFwd.sol->data.col(selSolIdcs[i]);
1072 
1073  selectedFwd.sol->data = G;
1074  selectedFwd.sol->nrow = selectedFwd.sol->data.rows();
1075  selectedFwd.sol->ncol = selectedFwd.sol->data.cols();
1076  selectedFwd.nsource = selectedFwd.sol->ncol / 3;
1077 
1078  selectedFwd.src = selectedFwd.src.pick_regions(p_qListLabels);
1079 
1080  return selectedFwd;
1081 }
1082 
1083 
1084 //*************************************************************************************************************
1085 
1086 MNEForwardSolution MNEForwardSolution::pick_types(bool meg, bool eeg, const QStringList& include, const QStringList& exclude) const
1087 {
1088  RowVectorXi sel = info.pick_types(meg, eeg, false, include, exclude);
1089 
1090  QStringList include_ch_names;
1091  for(qint32 i = 0; i < sel.cols(); ++i)
1092  include_ch_names << info.ch_names[sel[i]];
1093 
1094  return this->pick_channels(include_ch_names);
1095 }
1096 
1097 
1098 //*************************************************************************************************************
1099 
1100 void MNEForwardSolution::prepare_forward(const FiffInfo &p_info, const FiffCov &p_noise_cov, bool p_pca, FiffInfo &p_outFwdInfo, MatrixXd &gain, FiffCov &p_outNoiseCov, MatrixXd &p_outWhitener, qint32 &p_outNumNonZero) const
1101 {
1102  QStringList fwd_ch_names, ch_names;
1103  for(qint32 i = 0; i < this->info.chs.size(); ++i)
1104  fwd_ch_names << this->info.chs[i].ch_name;
1105 
1106  ch_names.clear();
1107  for(qint32 i = 0; i < p_info.chs.size(); ++i)
1108  if( !p_info.bads.contains(p_info.chs[i].ch_name)
1109  && !p_noise_cov.bads.contains(p_info.chs[i].ch_name)
1110  && fwd_ch_names.contains(p_info.chs[i].ch_name))
1111  ch_names << p_info.chs[i].ch_name;
1112 
1113  qint32 n_chan = ch_names.size();
1114  printf("Computing inverse operator with %d channels.\n", n_chan);
1115 
1116  //
1117  // Handle noise cov
1118  //
1119  p_outNoiseCov = p_noise_cov.prepare_noise_cov(p_info, ch_names);
1120 
1121  // Omit the zeroes due to projection
1122  p_outNumNonZero = 0;
1123  VectorXi t_vecNonZero = VectorXi::Zero(n_chan);
1124  for(qint32 i = 0; i < p_outNoiseCov.eig.rows(); ++i)
1125  {
1126  if(p_outNoiseCov.eig[i] > 0)
1127  {
1128  t_vecNonZero[p_outNumNonZero] = i;
1129  ++p_outNumNonZero;
1130  }
1131  }
1132  if(p_outNumNonZero > 0)
1133  t_vecNonZero.conservativeResize(p_outNumNonZero);
1134 
1135  if(p_outNumNonZero > 0)
1136  {
1137  if (p_pca)
1138  {
1139  qWarning("Warning in MNEForwardSolution::prepare_forward: if (p_pca) havent been debugged.");
1140  p_outWhitener = MatrixXd::Zero(n_chan, p_outNumNonZero);
1141  // Rows of eigvec are the eigenvectors
1142  for(qint32 i = 0; i < p_outNumNonZero; ++i)
1143  p_outWhitener.col(t_vecNonZero[i]) = p_outNoiseCov.eigvec.col(t_vecNonZero[i]).array() / sqrt(p_outNoiseCov.eig(t_vecNonZero[i]));
1144  printf("\tReducing data rank to %d.\n", p_outNumNonZero);
1145  }
1146  else
1147  {
1148  printf("Creating non pca whitener.\n");
1149  p_outWhitener = MatrixXd::Zero(n_chan, n_chan);
1150  for(qint32 i = 0; i < p_outNumNonZero; ++i)
1151  p_outWhitener(t_vecNonZero[i],t_vecNonZero[i]) = 1.0 / sqrt(p_outNoiseCov.eig(t_vecNonZero[i]));
1152  // Cols of eigvec are the eigenvectors
1153  p_outWhitener *= p_outNoiseCov.eigvec;
1154  }
1155  }
1156 
1157  VectorXi fwd_idx = VectorXi::Zero(ch_names.size());
1158  VectorXi info_idx = VectorXi::Zero(ch_names.size());
1159  qint32 idx;
1160  qint32 count_fwd_idx = 0;
1161  qint32 count_info_idx = 0;
1162  for(qint32 i = 0; i < ch_names.size(); ++i)
1163  {
1164  idx = fwd_ch_names.indexOf(ch_names[i]);
1165  if(idx > -1)
1166  {
1167  fwd_idx[count_fwd_idx] = idx;
1168  ++count_fwd_idx;
1169  }
1170  idx = p_info.ch_names.indexOf(ch_names[i]);
1171  if(idx > -1)
1172  {
1173  info_idx[count_info_idx] = idx;
1174  ++count_info_idx;
1175  }
1176  }
1177  fwd_idx.conservativeResize(count_fwd_idx);
1178  info_idx.conservativeResize(count_info_idx);
1179 
1180  gain.resize(count_fwd_idx, this->sol->data.cols());
1181  for(qint32 i = 0; i < count_fwd_idx; ++i)
1182  gain.row(i) = this->sol->data.row(fwd_idx[i]);
1183 
1184  p_outFwdInfo = p_info.pick_info(info_idx);
1185 
1186  printf("\tTotal rank is %d\n", p_outNumNonZero);
1187 }
1188 
1189 
1190 //*************************************************************************************************************
1191 
1192 bool MNEForwardSolution::read(QIODevice& p_IODevice, MNEForwardSolution& fwd, bool force_fixed, bool surf_ori, const QStringList& include, const QStringList& exclude, bool bExcludeBads)
1193 {
1194  FiffStream::SPtr t_pStream(new FiffStream(&p_IODevice));
1195  FiffDirTree t_Tree;
1196  QList<FiffDirEntry> t_Dir;
1197 
1198  printf("Reading forward solution from %s...\n", t_pStream->streamName().toUtf8().constData());
1199  if(!t_pStream->open(t_Tree, t_Dir))
1200  return false;
1201  //
1202  // Find all forward solutions
1203  //
1204  QList<FiffDirTree> fwds = t_Tree.dir_tree_find(FIFFB_MNE_FORWARD_SOLUTION);
1205 
1206  if (fwds.size() == 0)
1207  {
1208  t_pStream->device()->close();
1209  std::cout << "No forward solutions in " << t_pStream->streamName().toUtf8().constData(); // ToDo throw error
1210  return false;
1211  }
1212  //
1213  // Parent MRI data
1214  //
1215  QList<FiffDirTree> parent_mri = t_Tree.dir_tree_find(FIFFB_MNE_PARENT_MRI_FILE);
1216  if (parent_mri.size() == 0)
1217  {
1218  t_pStream->device()->close();
1219  std::cout << "No parent MRI information in " << t_pStream->streamName().toUtf8().constData(); // ToDo throw error
1220  return false;
1221  }
1222 
1223  MNESourceSpace t_SourceSpace;// = NULL;
1224  if(!MNESourceSpace::readFromStream(t_pStream, true, t_Tree, t_SourceSpace))
1225  {
1226  t_pStream->device()->close();
1227  std::cout << "Could not read the source spaces\n"; // ToDo throw error
1228  //ToDo error(me,'Could not read the source spaces (%s)',mne_omit_first_line(lasterr));
1229  return false;
1230  }
1231 
1232  for(qint32 k = 0; k < t_SourceSpace.size(); ++k)
1233  t_SourceSpace[k].id = MNESourceSpace::find_source_space_hemi(t_SourceSpace[k]);
1234 
1235  //
1236  // Bad channel list
1237  //
1238  QStringList bads;
1239  if(bExcludeBads)
1240  {
1241  bads = t_pStream->read_bad_channels(t_Tree);
1242  if(bads.size() > 0)
1243  {
1244  printf("\t%d bad channels ( ",bads.size());
1245  for(qint32 i = 0; i < bads.size(); ++i)
1246  printf("\"%s\" ", bads[i].toLatin1().constData());
1247  printf(") read\n");
1248  }
1249  }
1250 
1251  //
1252  // Locate and read the forward solutions
1253  //
1254  FiffTag::SPtr t_pTag;
1255  FiffDirTree megnode;
1256  FiffDirTree eegnode;
1257  for(qint32 k = 0; k < fwds.size(); ++k)
1258  {
1259  if(!fwds[k].find_tag(t_pStream.data(), FIFF_MNE_INCLUDED_METHODS, t_pTag))
1260  {
1261  t_pStream->device()->close();
1262  std::cout << "Methods not listed for one of the forward solutions\n"; // ToDo throw error
1263  return false;
1264  }
1265  if (*t_pTag->toInt() == FIFFV_MNE_MEG)
1266  {
1267  printf("MEG solution found\n");
1268  megnode = fwds[k];
1269  }
1270  else if(*t_pTag->toInt() == FIFFV_MNE_EEG)
1271  {
1272  printf("EEG solution found\n");
1273  eegnode = fwds.at(k);
1274  }
1275  }
1276 
1277  MNEForwardSolution megfwd;
1278  QString ori;
1279  if (read_one(t_pStream.data(), megnode, megfwd))
1280  {
1281  if (megfwd.source_ori == FIFFV_MNE_FIXED_ORI)
1282  ori = QString("fixed");
1283  else
1284  ori = QString("free");
1285  printf("\tRead MEG forward solution (%d sources, %d channels, %s orientations)\n", megfwd.nsource,megfwd.nchan,ori.toUtf8().constData());
1286  }
1287  MNEForwardSolution eegfwd;
1288  if (read_one(t_pStream.data(), eegnode, eegfwd))
1289  {
1290  if (eegfwd.source_ori == FIFFV_MNE_FIXED_ORI)
1291  ori = QString("fixed");
1292  else
1293  ori = QString("free");
1294  printf("\tRead EEG forward solution (%d sources, %d channels, %s orientations)\n", eegfwd.nsource,eegfwd.nchan,ori.toUtf8().constData());
1295  }
1296 
1297  //
1298  // Merge the MEG and EEG solutions together
1299  //
1300  fwd.clear();
1301 
1302  if (!megfwd.isEmpty() && !eegfwd.isEmpty())
1303  {
1304  if (megfwd.sol->data.cols() != eegfwd.sol->data.cols() ||
1305  megfwd.source_ori != eegfwd.source_ori ||
1306  megfwd.nsource != eegfwd.nsource ||
1307  megfwd.coord_frame != eegfwd.coord_frame)
1308  {
1309  t_pStream->device()->close();
1310  std::cout << "The MEG and EEG forward solutions do not match\n"; // ToDo throw error
1311  return false;
1312  }
1313 
1314  fwd = MNEForwardSolution(megfwd);
1315  fwd.sol->data = MatrixXd(megfwd.sol->nrow + eegfwd.sol->nrow, megfwd.sol->ncol);
1316 
1317  fwd.sol->data.block(0,0,megfwd.sol->nrow,megfwd.sol->ncol) = megfwd.sol->data;
1318  fwd.sol->data.block(megfwd.sol->nrow,0,eegfwd.sol->nrow,eegfwd.sol->ncol) = eegfwd.sol->data;
1319  fwd.sol->nrow = megfwd.sol->nrow + eegfwd.sol->nrow;
1320  fwd.sol->row_names.append(eegfwd.sol->row_names);
1321 
1322  if (!fwd.sol_grad->isEmpty())
1323  {
1324  fwd.sol_grad->data.resize(megfwd.sol_grad->data.rows() + eegfwd.sol_grad->data.rows(), megfwd.sol_grad->data.cols());
1325 
1326  fwd.sol->data.block(0,0,megfwd.sol_grad->data.rows(),megfwd.sol_grad->data.cols()) = megfwd.sol_grad->data;
1327  fwd.sol->data.block(megfwd.sol_grad->data.rows(),0,eegfwd.sol_grad->data.rows(),eegfwd.sol_grad->data.cols()) = eegfwd.sol_grad->data;
1328 
1329  fwd.sol_grad->nrow = megfwd.sol_grad->nrow + eegfwd.sol_grad->nrow;
1330  fwd.sol_grad->row_names.append(eegfwd.sol_grad->row_names);
1331  }
1332  fwd.nchan = megfwd.nchan + eegfwd.nchan;
1333  printf("\tMEG and EEG forward solutions combined\n");
1334  }
1335  else if (!megfwd.isEmpty())
1336  fwd = megfwd; //new MNEForwardSolution(megfwd);//not copied for the sake of speed
1337  else
1338  fwd = eegfwd; //new MNEForwardSolution(eegfwd);//not copied for the sake of speed
1339 
1340  //
1341  // Get the MRI <-> head coordinate transformation
1342  //
1343  if(!parent_mri[0].find_tag(t_pStream.data(), FIFF_COORD_TRANS, t_pTag))
1344  {
1345  t_pStream->device()->close();
1346  std::cout << "MRI/head coordinate transformation not found\n"; // ToDo throw error
1347  return false;
1348  }
1349  else
1350  {
1351  fwd.mri_head_t = t_pTag->toCoordTrans();
1352 
1353  if (fwd.mri_head_t.from != FIFFV_COORD_MRI || fwd.mri_head_t.to != FIFFV_COORD_HEAD)
1354  {
1356  if (fwd.mri_head_t.from != FIFFV_COORD_MRI || fwd.mri_head_t.to != FIFFV_COORD_HEAD)
1357  {
1358  t_pStream->device()->close();
1359  std::cout << "MRI/head coordinate transformation not found\n"; // ToDo throw error
1360  return false;
1361  }
1362  }
1363  }
1364 
1365  //
1366  // get parent MEG info -> from python package
1367  //
1368  t_pStream->read_meas_info_base(t_Tree, fwd.info);
1369 
1370 
1371  t_pStream->device()->close();
1372 
1373  //
1374  // Transform the source spaces to the correct coordinate frame
1375  // if necessary
1376  //
1377  if (fwd.coord_frame != FIFFV_COORD_MRI && fwd.coord_frame != FIFFV_COORD_HEAD)
1378  {
1379  std::cout << "Only forward solutions computed in MRI or head coordinates are acceptable";
1380  return false;
1381  }
1382 
1383  //
1384  qint32 nuse = 0;
1385  t_SourceSpace.transform_source_space_to(fwd.coord_frame,fwd.mri_head_t);
1386  for(qint32 k = 0; k < t_SourceSpace.size(); ++k)
1387  nuse += t_SourceSpace[k].nuse;
1388 
1389  if (nuse != fwd.nsource)
1390  throw("Source spaces do not match the forward solution.\n");
1391 
1392  printf("\tSource spaces transformed to the forward solution coordinate frame\n");
1393  fwd.src = t_SourceSpace; //not new MNESourceSpace(t_SourceSpace); for sake of speed
1394  //
1395  // Handle the source locations and orientations
1396  //
1397  if (fwd.isFixedOrient() || force_fixed == true)
1398  {
1399  nuse = 0;
1400  fwd.source_rr = MatrixXf::Zero(fwd.nsource,3);
1401  fwd.source_nn = MatrixXf::Zero(fwd.nsource,3);
1402  for(qint32 k = 0; k < t_SourceSpace.size();++k)
1403  {
1404  for(qint32 q = 0; q < t_SourceSpace[k].nuse; ++q)
1405  {
1406  fwd.source_rr.block(q,0,1,3) = t_SourceSpace[k].rr.block(t_SourceSpace[k].vertno(q),0,1,3);
1407  fwd.source_nn.block(q,0,1,3) = t_SourceSpace[k].nn.block(t_SourceSpace[k].vertno(q),0,1,3);
1408  }
1409  nuse += t_SourceSpace[k].nuse;
1410  }
1411  //
1412  // Modify the forward solution for fixed source orientations
1413  //
1414  if (fwd.source_ori != FIFFV_MNE_FIXED_ORI)
1415  {
1416  printf("\tChanging to fixed-orientation forward solution...");
1417 
1418  MatrixXd tmp = fwd.source_nn.transpose().cast<double>();
1419  SparseMatrix<double>* fix_rot = MNEMath::make_block_diag(tmp,1);
1420  fwd.sol->data *= (*fix_rot);
1421  fwd.sol->ncol = fwd.nsource;
1422  fwd.source_ori = FIFFV_MNE_FIXED_ORI;
1423 
1424  if (!fwd.sol_grad->isEmpty())
1425  {
1426  SparseMatrix<double> t_matKron;
1427  SparseMatrix<double> t_eye(3,3);
1428  for (qint32 i = 0; i < 3; ++i)
1429  t_eye.insert(i,i) = 1.0f;
1430  t_matKron = kroneckerProduct(*fix_rot,t_eye);//kron(fix_rot,eye(3));
1431  fwd.sol_grad->data *= t_matKron;
1432  fwd.sol_grad->ncol = 3*fwd.nsource;
1433  }
1434  delete fix_rot;
1435  printf("[done]\n");
1436  }
1437  }
1438  else if (surf_ori)
1439  {
1440  //
1441  // Rotate the local source coordinate systems
1442  //
1443  printf("\tConverting to surface-based source orientations...");
1444 
1445  bool use_ave_nn = false;
1446  if(t_SourceSpace[0].patch_inds.size() > 0)
1447  {
1448  use_ave_nn = true;
1449  printf("\tAverage patch normals will be employed in the rotation to the local surface coordinates...\n");
1450  }
1451 
1452  nuse = 0;
1453  qint32 pp = 0;
1454  fwd.source_rr = MatrixXf::Zero(fwd.nsource,3);
1455  fwd.source_nn = MatrixXf::Zero(fwd.nsource*3,3);
1456 
1457  qWarning("Warning source_ori: Rotating the source coordinate system haven't been verified --> Singular Vectors U are different from MATLAB!");
1458 
1459  for(qint32 k = 0; k < t_SourceSpace.size();++k)
1460  {
1461 
1462  for (qint32 q = 0; q < t_SourceSpace[k].nuse; ++q)
1463  fwd.source_rr.block(q+nuse,0,1,3) = t_SourceSpace[k].rr.block(t_SourceSpace[k].vertno(q),0,1,3);
1464 
1465  for (qint32 p = 0; p < t_SourceSpace[k].nuse; ++p)
1466  {
1467  //
1468  // Project out the surface normal and compute SVD
1469  //
1470  Vector3f nn;
1471  if(use_ave_nn)
1472  {
1473  VectorXi t_vIdx = t_SourceSpace[k].pinfo[t_SourceSpace[k].patch_inds[p]];
1474  Matrix3Xf t_nn(3, t_vIdx.size());
1475  for(qint32 i = 0; i < t_vIdx.size(); ++i)
1476  t_nn.col(i) = t_SourceSpace[k].nn.block(t_vIdx[i],0,1,3).transpose();
1477  nn = t_nn.rowwise().sum();
1478  nn.array() /= nn.norm();
1479  }
1480  else
1481  nn = t_SourceSpace[k].nn.block(t_SourceSpace[k].vertno(p),0,1,3).transpose();
1482 
1483  Matrix3f tmp = Matrix3f::Identity(nn.rows(), nn.rows()) - nn*nn.transpose();
1484 
1485  JacobiSVD<MatrixXf> t_svd(tmp, Eigen::ComputeThinU);
1486  //Sort singular values and singular vectors
1487  VectorXf t_s = t_svd.singularValues();
1488  MatrixXf U = t_svd.matrixU();
1489  MNEMath::sort<float>(t_s, U);
1490 
1491  //
1492  // Make sure that ez is in the direction of nn
1493  //
1494  if ((nn.transpose() * U.block(0,2,3,1))(0,0) < 0)
1495  U *= -1;
1496  fwd.source_nn.block(pp, 0, 3, 3) = U.transpose();
1497  pp += 3;
1498  }
1499  nuse += t_SourceSpace[k].nuse;
1500  }
1501  MatrixXd tmp = fwd.source_nn.transpose().cast<double>();
1502  SparseMatrix<double>* surf_rot = MNEMath::make_block_diag(tmp,3);
1503 
1504  fwd.sol->data *= *surf_rot;
1505 
1506  if (!fwd.sol_grad->isEmpty())
1507  {
1508  SparseMatrix<double> t_matKron;
1509  SparseMatrix<double> t_eye(3,3);
1510  for (qint32 i = 0; i < 3; ++i)
1511  t_eye.insert(i,i) = 1.0f;
1512  t_matKron = kroneckerProduct(*surf_rot,t_eye);//kron(surf_rot,eye(3));
1513  fwd.sol_grad->data *= t_matKron;
1514  }
1515  delete surf_rot;
1516  printf("[done]\n");
1517  }
1518  else
1519  {
1520  printf("\tCartesian source orientations...");
1521  nuse = 0;
1522  fwd.source_rr = MatrixXf::Zero(fwd.nsource,3);
1523  for(qint32 k = 0; k < t_SourceSpace.size(); ++k)
1524  {
1525  for (qint32 q = 0; q < t_SourceSpace[k].nuse; ++q)
1526  fwd.source_rr.block(q+nuse,0,1,3) = t_SourceSpace[k].rr.block(t_SourceSpace[k].vertno(q),0,1,3);
1527 
1528  nuse += t_SourceSpace[k].nuse;
1529  }
1530 
1531  MatrixXf t_ones = MatrixXf::Ones(fwd.nsource,1);
1532  Matrix3f t_eye = Matrix3f::Identity();
1533  fwd.source_nn = kroneckerProduct(t_ones,t_eye);
1534 
1535  printf("[done]\n");
1536  }
1537 
1538  //
1539  // Do the channel selection
1540  //
1541  QStringList exclude_bads = exclude;
1542  if (bads.size() > 0)
1543  {
1544  for(qint32 k = 0; k < bads.size(); ++k)
1545  if(!exclude_bads.contains(bads[k],Qt::CaseInsensitive))
1546  exclude_bads << bads[k];
1547  }
1548 
1549  fwd.surf_ori = surf_ori;
1550  fwd = fwd.pick_channels(include, exclude_bads);
1551 
1552 // //
1553 // // Do the channel selection - OLD VERSION
1554 // //
1555 // if (include.size() > 0 || exclude.size() > 0 || bads.size() > 0)
1556 // {
1557 // //
1558 // // First do the channels to be included
1559 // //
1560 // RowVectorXi pick;
1561 // if (include.size() == 0)
1562 // pick = RowVectorXi::Ones(fwd.nchan);
1563 // else
1564 // {
1565 // pick = RowVectorXi::Zero(fwd.nchan);
1566 // for(qint32 k = 0; k < include.size(); ++k)
1567 // {
1568 // QList<int> c;
1569 // for(qint32 l = 0; l < fwd.sol->row_names.size(); ++l)
1570 // if(fwd.sol->row_names.at(l).contains(include.at(k),Qt::CaseInsensitive))
1571 // pick(l) = 1;
1572 
1577 // }
1578 // }
1579 // //
1580 // // Then exclude what needs to be excluded
1581 // //
1582 // if (exclude.size() > 0)
1583 // {
1584 // for(qint32 k = 0; k < exclude.size(); ++k)
1585 // {
1586 // QList<int> c;
1587 // for(qint32 l = 0; l < fwd.sol->row_names.size(); ++l)
1588 // if(fwd.sol->row_names.at(l).contains(exclude.at(k),Qt::CaseInsensitive))
1589 // pick(l) = 0;
1590 
1595 // }
1596 // }
1597 // if (bads.size() > 0)
1598 // {
1599 // for(qint32 k = 0; k < bads.size(); ++k)
1600 // {
1601 // QList<int> c;
1602 // for(qint32 l = 0; l < fwd.sol->row_names.size(); ++l)
1603 // if(fwd.sol->row_names.at(l).contains(bads.at(k),Qt::CaseInsensitive))
1604 // pick(l) = 0;
1605 
1610 // }
1611 // }
1612 // //
1613 // // Do we have something?
1614 // //
1615 // nuse = pick.sum();
1616 // if (nuse == 0)
1617 // throw("Nothing remains after picking");
1618 // //
1619 // // Create the selector
1620 // //
1621 // qint32 p = 0;
1622 // MatrixXd t_solData(nuse,fwd.sol->data.cols());
1623 // QStringList t_solRowNames;
1624 
1625 // MatrixXd t_solGradData;// = NULL;
1626 // QStringList t_solGradRowNames;
1627 
1628 // QList<FiffChInfo> chs;
1629 
1630 
1631 // if (!fwd.sol_grad->isEmpty())
1632 // t_solGradData.resize(nuse, fwd.sol_grad->data.cols());
1633 
1634 // for(qint32 k = 0; k < fwd.nchan; ++k)
1635 // {
1636 // if(pick(k) > 0)
1637 // {
1638 // t_solData.block(p, 0, 1, fwd.sol->data.cols()) = fwd.sol->data.block(k, 0, 1, fwd.sol->data.cols());
1639 // t_solRowNames.append(fwd.sol->row_names[k]);
1640 
1641 // if (!fwd.sol_grad->isEmpty())
1642 // {
1643 // t_solGradData.block(p, 0, 1, fwd.sol_grad->data.cols()) = fwd.sol_grad->data.block(k, 0, 1, fwd.sol_grad->data.cols());
1644 // t_solGradRowNames.append(fwd.sol_grad->row_names[k]);
1645 // }
1646 
1647 // chs.append(fwd.info.chs[k]);
1648 
1649 // ++p;
1650 // }
1651 // }
1652 // printf("\t%d out of %d channels remain after picking\n",nuse,fwd.nchan);
1653 // //
1654 // // Pick the correct rows of the forward operator
1655 // //
1656 // fwd.nchan = nuse;
1659 // fwd.sol->data = t_solData;
1660 // fwd.sol->nrow = nuse;
1661 // fwd.sol->row_names = t_solRowNames;
1662 
1663 // if (!fwd.sol_grad->isEmpty())
1664 // {
1665 // fwd.sol_grad->data = t_solGradData;
1666 // fwd.sol_grad->nrow = nuse;
1667 // fwd.sol_grad->row_names = t_solGradRowNames;
1668 // }
1669 
1670 // fwd.info.chs = chs;
1671 // }
1672 
1673  //garbage collecting
1674  t_pStream->device()->close();
1675 
1676  return true;
1677 }
1678 
1679 
1680 //*************************************************************************************************************
1681 
1682 bool MNEForwardSolution::read_one(FiffStream* p_pStream, const FiffDirTree& p_Node, MNEForwardSolution& one)
1683 {
1684  //
1685  // Read all interesting stuff for one forward solution
1686  //
1687  if(p_Node.isEmpty())
1688  return false;
1689 
1690  one.clear();
1691  FiffTag::SPtr t_pTag;
1692 
1693  if(!p_Node.find_tag(p_pStream, FIFF_MNE_SOURCE_ORIENTATION, t_pTag))
1694  {
1695  p_pStream->device()->close();
1696  std::cout << "Source orientation tag not found."; //ToDo: throw error.
1697  return false;
1698  }
1699 
1700  one.source_ori = *t_pTag->toInt();
1701 
1702  if(!p_Node.find_tag(p_pStream, FIFF_MNE_COORD_FRAME, t_pTag))
1703  {
1704  p_pStream->device()->close();
1705  std::cout << "Coordinate frame tag not found."; //ToDo: throw error.
1706  return false;
1707  }
1708 
1709  one.coord_frame = *t_pTag->toInt();
1710 
1711  if(!p_Node.find_tag(p_pStream, FIFF_MNE_SOURCE_SPACE_NPOINTS, t_pTag))
1712  {
1713  p_pStream->device()->close();
1714  std::cout << "Number of sources not found."; //ToDo: throw error.
1715  return false;
1716  }
1717 
1718  one.nsource = *t_pTag->toInt();
1719 
1720  if(!p_Node.find_tag(p_pStream, FIFF_NCHAN, t_pTag))
1721  {
1722  p_pStream->device()->close();
1723  printf("Number of channels not found."); //ToDo: throw error.
1724  return false;
1725  }
1726 
1727  one.nchan = *t_pTag->toInt();
1728 
1729  if(p_pStream->read_named_matrix(p_Node, FIFF_MNE_FORWARD_SOLUTION, *one.sol.data()))
1730  one.sol->transpose_named_matrix();
1731  else
1732  {
1733  p_pStream->device()->close();
1734  printf("Forward solution data not found ."); //ToDo: throw error.
1735  //error(me,'Forward solution data not found (%s)',mne_omit_first_line(lasterr));
1736  return false;
1737  }
1738 
1739  if(p_pStream->read_named_matrix(p_Node, FIFF_MNE_FORWARD_SOLUTION_GRAD, *one.sol_grad.data()))
1740  one.sol_grad->transpose_named_matrix();
1741  else
1742  one.sol_grad->clear();
1743 
1744 
1745  if (one.sol->data.rows() != one.nchan ||
1746  (one.sol->data.cols() != one.nsource && one.sol->data.cols() != 3*one.nsource))
1747  {
1748  p_pStream->device()->close();
1749  printf("Forward solution matrix has wrong dimensions.\n"); //ToDo: throw error.
1750  //error(me,'Forward solution matrix has wrong dimensions');
1751  return false;
1752  }
1753  if (!one.sol_grad->isEmpty())
1754  {
1755  if (one.sol_grad->data.rows() != one.nchan ||
1756  (one.sol_grad->data.cols() != 3*one.nsource && one.sol_grad->data.cols() != 3*3*one.nsource))
1757  {
1758  p_pStream->device()->close();
1759  printf("Forward solution gradient matrix has wrong dimensions.\n"); //ToDo: throw error.
1760  //error(me,'Forward solution gradient matrix has wrong dimensions');
1761  }
1762  }
1763  return true;
1764 }
1765 
1766 
1767 //*************************************************************************************************************
1768 
1770 {
1771  // Figure out which ones have been used
1772  if(info.chs.size() != G.rows())
1773  {
1774  printf("Error G.rows() and length of info.chs do not match: %li != %i", G.rows(), info.chs.size()); //ToDo throw
1775  return;
1776  }
1777 
1778  RowVectorXi sel = info.pick_types(QString("grad"));
1779  if(sel.size() > 0)
1780  {
1781  for(qint32 i = 0; i < sel.size(); ++i)
1782  G.row(i) = G.row(sel[i]);
1783  G.conservativeResize(sel.size(), G.cols());
1784  printf("\t%li planar channels", sel.size());
1785  }
1786  else
1787  {
1788  sel = info.pick_types(QString("mag"));
1789  if (sel.size() > 0)
1790  {
1791  for(qint32 i = 0; i < sel.size(); ++i)
1792  G.row(i) = G.row(sel[i]);
1793  G.conservativeResize(sel.size(), G.cols());
1794  printf("\t%li magnetometer or axial gradiometer channels", sel.size());
1795  }
1796  else
1797  {
1798  sel = info.pick_types(false, true);
1799  if(sel.size() > 0)
1800  {
1801  for(qint32 i = 0; i < sel.size(); ++i)
1802  G.row(i) = G.row(sel[i]);
1803  G.conservativeResize(sel.size(), G.cols());
1804  printf("\t%li EEG channels\n", sel.size());
1805  }
1806  else
1807  printf("Could not find MEG or EEG channels\n");
1808  }
1809  }
1810 }
1811 
1812 
1813 //*************************************************************************************************************
1814 
1816 {
1817  if(!this->surf_ori || this->isFixedOrient())
1818  {
1819  qWarning("Warning: Only surface-oriented, free-orientation forward solutions can be converted to fixed orientaton.\n");//ToDo: Throw here//qCritical//qFatal
1820  return;
1821  }
1822  qint32 count = 0;
1823  for(qint32 i = 2; i < this->sol->data.cols(); i += 3)
1824  this->sol->data.col(count) = this->sol->data.col(i);//ToDo: is this right? - just take z?
1825  this->sol->data.conservativeResize(this->sol->data.rows(), count);
1826  this->sol->ncol = this->sol->ncol / 3;
1827  this->source_ori = FIFFV_MNE_FIXED_ORI;
1828  printf("\tConverted the forward solution into the fixed-orientation mode.\n");
1829 }
static qint32 find_source_space_hemi(MNEHemisphere &p_Hemisphere)
#define FIFF_MNE_SOURCE_SPACE_NPOINTS
MatrixXd data
Definition: fiff_cov.h:211
MNEForwardSolution pick_types(bool meg, bool eeg, const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList) const
#define FIFF_MNE_SOURCE_ORIENTATION
FIFF measurement file information.
Definition: fiff_info.h:96
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)
fiff_int_t kind
Definition: fiff_cov.h:207
MNEForwardSolution cluster_forward_solution(const AnnotationSet &p_AnnotationSet, qint32 p_iClusterSize, MatrixXd &p_D=defaultD, const FiffCov &p_pNoise_cov=defaultCov, const FiffInfo &p_pInfo=defaultInfo, QString p_sMethod="cityblock") const
QList< FiffDirTree > dir_tree_find(fiff_int_t p_kind) const
QSharedDataPointer< FiffNamedMatrix > SDPtr
VectorXi getLabelIds() const
Definition: colortable.h:133
FiffInfo pick_info(const RowVectorXi &sel=defaultVectorXi) const
Definition: fiff_info.cpp:300
fiff_int_t nfree
Definition: fiff_cov.h:214
QSharedPointer< FiffTag > SPtr
Definition: fiff_tag.h:166
Colortable class declaration.
static void restrict_gain_matrix(MatrixXd &G, const FiffInfo &info)
QStringList bads
Definition: fiff_cov.h:213
VectorXi tripletSelection(const VectorXi &p_vecIdxSelection) const
static bool readFromStream(FiffStream::SPtr &p_pStream, bool add_geom, FiffDirTree &p_Tree, MNESourceSpace &p_SourceSpace)
QStringList struct_names
Definition: colortable.h:124
static SparseMatrix< double > * make_block_diag(const MatrixXd &A, qint32 n)
Definition: mnemath.cpp:295
QSharedPointer< FiffStream > SPtr
Definition: fiff_stream.h:132
bool read_named_matrix(const FiffDirTree &p_Node, fiff_int_t matkind, FiffNamedMatrix &mat)
static FiffCov compute_depth_prior(const MatrixXd &Gain, const FiffInfo &gain_info, bool is_fixed_ori, double exp=0.8, double limit=10.0, const MatrixXd &patch_areas=defaultConstMatrixXd, bool limit_depth_chs=false)
FiffNamedMatrix::SDPtr sol_grad
MNEForwardSolution class declaration, which provides the forward solution including the source space ...
static VectorXi intersect(const VectorXi &v1, const VectorXi &v2, VectorXi &idx_sel)
Definition: mnemath.cpp:158
Annotation set.
Definition: annotationset.h:96
bool isEmpty() const
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 FIFFV_MNE_ORIENT_PRIOR_COV
bool isEmpty() const
Definition: fiff_cov.h:225
QList< FiffChInfo > chs
fiff_int_t dim
Definition: fiff_cov.h:209
covariance data
Definition: fiff_cov.h:94
FiffNamedMatrix::SDPtr sol
MatrixXd eigvec
Definition: fiff_cov.h:216
#define FIFF_MNE_COORD_FRAME
MNEForwardSolution pick_regions(const QList< Label > &p_qListLabels) const
MNEForwardSolution pick_channels(const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList) const
RowVectorXi pick_types(const QString meg, bool eeg=false, bool stim=false, const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList) const
Vertices label based lookup table.
Definition: colortable.h:79
FiffCov prepare_noise_cov(const FiffInfo &p_info, const QStringList &p_chNames) const
Definition: fiff_cov.cpp:186
static bool read(QIODevice &p_IODevice, MNEForwardSolution &fwd, bool force_fixed=false, bool surf_ori=false, const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList, bool bExcludeBads=true)
KMeans class declaration.
QList< VectorXi > get_vertno() const
MNESourceSpace pick_regions(const QList< Label > &p_qListLabels) const
void prepare_forward(const FiffInfo &p_info, const FiffCov &p_noise_cov, bool p_pca, FiffInfo &p_outFwdInfo, MatrixXd &gain, FiffCov &p_outNoiseCov, MatrixXd &p_outWhitener, qint32 &p_outNumNonZero) const
FiffCov compute_orient_prior(float loose=0.2)
Label class declaration.
MNEMath class declaration.
MNEForwardSolution reduce_forward_solution(qint32 p_iNumDipoles, MatrixXd &p_D) const
FIFF File I/O routines.
Definition: fiff_stream.h:129
VectorXd eig
Definition: fiff_cov.h:215
static RowVectorXi pick_channels(const QStringList &ch_names, const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList)
#define FIFFV_MNE_DEPTH_PRIOR_COV
QStringList getNames() const
Definition: colortable.h:145
QList< VectorXi > label_src_vertno_sel(const Label &p_label, VectorXi &src_sel) const