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