MNE-CPP  0.1.9
A Framework for Electrophysiology
mne_inverse_operator.cpp
1 //=============================================================================================================
38 //=============================================================================================================
39 // INCLUDES
40 //=============================================================================================================
41 
42 #include "mne_inverse_operator.h"
43 #include <fs/label.h>
44 
45 #include <iostream>
46 
47 //=============================================================================================================
48 // QT INCLUDES
49 //=============================================================================================================
50 
51 #include <QFuture>
52 #include <QtConcurrent>
53 
54 //=============================================================================================================
55 // EIGEN INCLUDES
56 //=============================================================================================================
57 
58 #include <Eigen/SVD>
59 
60 //=============================================================================================================
61 // USED NAMESPACES
62 //=============================================================================================================
63 
64 using namespace UTILSLIB;
65 using namespace MNELIB;
66 using namespace FIFFLIB;
67 using namespace FSLIB;
68 using namespace Eigen;
69 
70 //=============================================================================================================
71 // DEFINE MEMBER METHODS
72 //=============================================================================================================
73 
74 MNEInverseOperator::MNEInverseOperator()
75 : methods(-1)
76 , source_ori(-1)
77 , nsource(-1)
78 , nchan(-1)
79 , coord_frame(-1)
80 , eigen_leads_weighted(false)
81 , eigen_leads(new FiffNamedMatrix)
82 , eigen_fields(new FiffNamedMatrix)
83 , noise_cov(new FiffCov)
84 , source_cov(new FiffCov)
85 , orient_prior(new FiffCov)
86 , depth_prior(new FiffCov)
87 , fmri_prior(new FiffCov)
88 , nave(-1)
89 {
90  qRegisterMetaType<QSharedPointer<MNELIB::MNEInverseOperator> >("QSharedPointer<MNELIB::MNEInverseOperator>");
91  qRegisterMetaType<MNELIB::MNEInverseOperator>("MNELIB::MNEInverseOperator");
92 }
93 
94 //=============================================================================================================
95 
97 {
99  qRegisterMetaType<QSharedPointer<MNELIB::MNEInverseOperator> >("QSharedPointer<MNELIB::MNEInverseOperator>");
100  qRegisterMetaType<MNELIB::MNEInverseOperator>("MNELIB::MNEInverseOperator");
101 }
102 
103 //=============================================================================================================
104 
106  const MNEForwardSolution& forward,
107  const FiffCov& p_noise_cov,
108  float loose,
109  float depth,
110  bool fixed,
111  bool limit_depth_chs)
112 {
113  *this = MNEInverseOperator::make_inverse_operator(info, forward, p_noise_cov, loose, depth, fixed, limit_depth_chs);
114  qRegisterMetaType<QSharedPointer<MNELIB::MNEInverseOperator> >("QSharedPointer<MNELIB::MNEInverseOperator>");
115  qRegisterMetaType<MNELIB::MNEInverseOperator>("MNELIB::MNEInverseOperator");
116 }
117 
118 //=============================================================================================================
119 
121 : info(p_MNEInverseOperator.info)
122 , methods(p_MNEInverseOperator.methods)
123 , source_ori(p_MNEInverseOperator.source_ori)
124 , nsource(p_MNEInverseOperator.nsource)
125 , nchan(p_MNEInverseOperator.nchan)
126 , coord_frame(p_MNEInverseOperator.coord_frame)
127 , source_nn(p_MNEInverseOperator.source_nn)
128 , sing(p_MNEInverseOperator.sing)
129 , eigen_leads_weighted(p_MNEInverseOperator.eigen_leads_weighted)
130 , eigen_leads(p_MNEInverseOperator.eigen_leads)
131 , eigen_fields(p_MNEInverseOperator.eigen_fields)
132 , noise_cov(p_MNEInverseOperator.noise_cov)
133 , source_cov(p_MNEInverseOperator.source_cov)
134 , orient_prior(p_MNEInverseOperator.orient_prior)
135 , depth_prior(p_MNEInverseOperator.depth_prior)
136 , fmri_prior(p_MNEInverseOperator.fmri_prior)
137 , src(p_MNEInverseOperator.src)
138 , mri_head_t(p_MNEInverseOperator.mri_head_t)
139 , nave(p_MNEInverseOperator.nave)
140 , projs(p_MNEInverseOperator.projs)
141 , proj(p_MNEInverseOperator.proj)
142 , whitener(p_MNEInverseOperator.whitener)
143 , reginv(p_MNEInverseOperator.reginv)
144 , noisenorm(p_MNEInverseOperator.noisenorm)
145 {
146  qRegisterMetaType<QSharedPointer<MNELIB::MNEInverseOperator> >("QSharedPointer<MNELIB::MNEInverseOperator>");
147  qRegisterMetaType<MNELIB::MNEInverseOperator>("MNELIB::MNEInverseOperator");
148 }
149 
150 //=============================================================================================================
151 
153 {
154 }
155 
156 //=============================================================================================================
157 
159  QString method,
160  bool pick_normal,
161  MatrixXd &K,
162  SparseMatrix<double> &noise_norm,
163  QList<VectorXi> &vertno)
164 {
165  MatrixXd t_eigen_leads = this->eigen_leads->data;
166  MatrixXd t_source_cov = this->source_cov->data;
167  if(method.compare("MNE") != 0)
168  noise_norm = this->noisenorm;
169 
170  vertno = this->src.get_vertno();
171 
172  typedef Eigen::Triplet<double> T;
173  std::vector<T> tripletList;
174 
175  if(!label.isEmpty())
176  {
177  qDebug() << "ToDo: Label selection needs to be debugged - not done jet!";
178  VectorXi src_sel;
179  vertno = this->src.label_src_vertno_sel(label, src_sel);
180 
181  if(method.compare("MNE") != 0)
182  {
183  tripletList.clear();
184  tripletList.reserve(noise_norm.nonZeros());
185 
186  for (qint32 k = 0; k < noise_norm.outerSize(); ++k)
187  {
188  for (SparseMatrix<double>::InnerIterator it(noise_norm,k); it; ++it)
189  {
190  //ToDo bad search - increase speed by using stl search
191  qint32 row = -1;
192  for(qint32 i = 0; i < src_sel.size(); ++i)
193  if(src_sel[i] == it.row())
194  row = i;
195  if(row != -1)
196  tripletList.push_back(T(it.row(), it.col(), it.value()));
197  }
198  }
199 
200  noise_norm = SparseMatrix<double>(src_sel.size(),noise_norm.cols());
201  noise_norm.setFromTriplets(tripletList.begin(), tripletList.end());
202  }
203 
204  if(this->source_ori == FIFFV_MNE_FREE_ORI)
205  {
206  VectorXi src_sel_new(src_sel.size()*3);
207 
208  for(qint32 i = 0; i < src_sel.size(); ++i)
209  {
210  src_sel_new[i*3] = src_sel[i]*3;
211  src_sel_new[i*3+1] = src_sel[i]*3+1;
212  src_sel_new[i*3+2] = src_sel[i]*3+2;
213  }
214  src_sel = src_sel_new;
215  }
216 
217  for(qint32 i = 0; i < src_sel.size(); ++i)
218  {
219  t_eigen_leads.row(i) = t_eigen_leads.row(src_sel[i]);
220  t_source_cov = t_source_cov.row(src_sel[i]);
221  }
222  t_eigen_leads.conservativeResize(src_sel.size(), t_eigen_leads.cols());
223  t_source_cov.conservativeResize(src_sel.size(), t_source_cov.cols());
224  }
225 
226  if(pick_normal)
227  {
228  if(this->source_ori != FIFFV_MNE_FREE_ORI)
229  {
230  qWarning("Warning: Pick normal can only be used with a free orientation inverse operator.\n");
231  return false;
232  }
233 
234  bool is_loose = ((0 < this->orient_prior->data(0,0)) && (this->orient_prior->data(0,0) < 1)) ? true : false;
235  if(!is_loose)
236  {
237  qWarning("The pick_normal parameter is only valid when working with loose orientations.\n");
238  return false;
239  }
240 
241  // keep only the normal components
242  qint32 count = 0;
243  for(qint32 i = 2; i < t_eigen_leads.rows(); i+=3)
244  {
245  t_eigen_leads.row(count) = t_eigen_leads.row(i);
246  ++count;
247  }
248  t_eigen_leads.conservativeResize(count, t_eigen_leads.cols());
249 
250  count = 0;
251  for(qint32 i = 2; i < t_source_cov.rows(); i+=3)
252  {
253  t_source_cov.row(count) = t_source_cov.row(i);
254  ++count;
255  }
256  t_source_cov.conservativeResize(count, t_source_cov.cols());
257  }
258 
259  tripletList.clear();
260  tripletList.reserve(reginv.rows());
261  for(qint32 i = 0; i < reginv.rows(); ++i)
262  tripletList.push_back(T(i, i, reginv(i,0)));
263  SparseMatrix<double> t_reginv(reginv.rows(),reginv.rows());
264  t_reginv.setFromTriplets(tripletList.begin(), tripletList.end());
265 
266  MatrixXd trans = t_reginv*eigen_fields->data*whitener*proj;
267  //
268  // Transformation into current distributions by weighting the eigenleads
269  // with the weights computed above
270  //
272  {
273  //
274  // R^0.5 has been already factored in
275  //
276  printf("(eigenleads already weighted)...\n");
277  K = t_eigen_leads*trans;
278  }
279  else
280  {
281  //
282  // R^0.5 has to factored in
283  //
284  printf("(eigenleads need to be weighted)...\n");
285 
286  std::vector<T> tripletList2;
287  tripletList2.reserve(t_source_cov.rows());
288  for(qint32 i = 0; i < t_source_cov.rows(); ++i)
289  tripletList2.push_back(T(i, i, sqrt(t_source_cov(i,0))));
290  SparseMatrix<double> t_sourceCov(t_source_cov.rows(),t_source_cov.rows());
291  t_sourceCov.setFromTriplets(tripletList2.begin(), tripletList2.end());
292 
293  K = t_sourceCov*t_eigen_leads*trans;
294  }
295 
296  if(method.compare("MNE") == 0)
297  noise_norm = SparseMatrix<double>();
298 
299  //store assembled kernel
300  m_K = K;
301 
302  return true;
303 }
304 
305 //=============================================================================================================
306 
308 {
309  QStringList inv_ch_names = this->eigen_fields->col_names;
310 
311  bool t_bContains = true;
312  if(this->eigen_fields->col_names.size() != this->noise_cov->names.size())
313  t_bContains = false;
314  else
315  {
316  for(qint32 i = 0; i < this->noise_cov->names.size(); ++i)
317  {
318  if(inv_ch_names[i] != this->noise_cov->names[i])
319  {
320  t_bContains = false;
321  break;
322  }
323  }
324  }
325 
326  if(!t_bContains)
327  {
328  qCritical("Channels in inverse operator eigen fields do not match noise covariance channels.");
329  return false;
330  }
331 
332  QStringList data_ch_names = info.ch_names;
333 
334  QStringList missing_ch_names;
335  for(qint32 i = 0; i < inv_ch_names.size(); ++i)
336  if(!data_ch_names.contains(inv_ch_names[i]))
337  missing_ch_names.append(inv_ch_names[i]);
338 
339  qint32 n_missing = missing_ch_names.size();
340 
341  if(n_missing > 0)
342  {
343  qCritical() << n_missing << "channels in inverse operator are not present in the data (" << missing_ch_names << ")";
344  return false;
345  }
346 
347  return true;
348 }
349 
350 //=============================================================================================================
351 
352 MatrixXd MNEInverseOperator::cluster_kernel(const AnnotationSet &p_AnnotationSet, qint32 p_iClusterSize, MatrixXd& p_D, QString p_sMethod) const
353 {
354  printf("Cluster kernel using %s.\n", p_sMethod.toUtf8().constData());
355 
356  MatrixXd p_outMT = this->m_K.transpose();
357 
358  QList<MNEClusterInfo> t_qListMNEClusterInfo;
359  MNEClusterInfo t_MNEClusterInfo;
360  t_qListMNEClusterInfo.append(t_MNEClusterInfo);
361  t_qListMNEClusterInfo.append(t_MNEClusterInfo);
362 
363  //
364  // Check consisty
365  //
366  if(this->isFixedOrient())
367  {
368  printf("Error: Fixed orientation not implemented jet!\n");
369  return p_outMT;
370  }
371 
372 // qDebug() << "p_outMT" << p_outMT.rows() << "x" << p_outMT.cols();
373 
374 // MatrixXd t_G_Whitened(0,0);
375 // bool t_bUseWhitened = false;
376 // //
377 // //Whiten gain matrix before clustering -> cause diffenerent units Magnetometer, Gradiometer and EEG
378 // //
379 // if(!p_pNoise_cov.isEmpty() && !p_pInfo.isEmpty())
380 // {
381 // FiffInfo p_outFwdInfo;
382 // FiffCov p_outNoiseCov;
383 // MatrixXd p_outWhitener;
384 // qint32 p_outNumNonZero;
385 // //do whitening with noise cov
386 // this->prepare_forward(p_pInfo, p_pNoise_cov, false, p_outFwdInfo, t_G_Whitened, p_outNoiseCov, p_outWhitener, p_outNumNonZero);
387 // printf("\tWhitening the forward solution.\n");
388 
389 // t_G_Whitened = p_outWhitener*t_G_Whitened;
390 // t_bUseWhitened = true;
391 // }
392 
393  //
394  // Assemble input data
395  //
396  qint32 count;
397  qint32 offset;
398 
399  MatrixXd t_MT_new;
400 
401  for(qint32 h = 0; h < this->src.size(); ++h )
402  {
403 
404  count = 0;
405  offset = 0;
406 
407  // Offset for continuous indexing;
408  if(h > 0)
409  for(qint32 j = 0; j < h; ++j)
410  offset += this->src[j].nuse;
411 
412  if(h == 0)
413  printf("Cluster Left Hemisphere\n");
414  else
415  printf("Cluster Right Hemisphere\n");
416 
417  Colortable t_CurrentColorTable = p_AnnotationSet[h].getColortable();
418  VectorXi label_ids = t_CurrentColorTable.getLabelIds();
419 
420  // Get label ids for every vertex
421  VectorXi vertno_labeled = VectorXi::Zero(this->src[h].vertno.rows());
422 
423  //ToDo make this more universal -> using Label instead of annotations - obsolete when using Labels
424  for(qint32 i = 0; i < vertno_labeled.rows(); ++i)
425  vertno_labeled[i] = p_AnnotationSet[h].getLabelIds()[this->src[h].vertno[i]];
426 
427  //Qt Concurrent List
428  QList<RegionMT> m_qListRegionMTIn;
429 
430  //
431  // Generate cluster input data
432  //
433  for (qint32 i = 0; i < label_ids.rows(); ++i)
434  {
435  if (label_ids[i] != 0)
436  {
437  QString curr_name = t_CurrentColorTable.struct_names[i];//obj.label2AtlasName(label(i));
438  printf("\tCluster %d / %ld %s...", i+1, label_ids.rows(), curr_name.toUtf8().constData());
439 
440  //
441  // Get source space indeces
442  //
443  VectorXi idcs = VectorXi::Zero(vertno_labeled.rows());
444  qint32 c = 0;
445 
446  //Select ROIs //change this use label info with a hash tabel
447  for(qint32 j = 0; j < vertno_labeled.rows(); ++j)
448  {
449  if(vertno_labeled[j] == label_ids[i])
450  {
451  idcs[c] = j;
452  ++c;
453  }
454  }
455  idcs.conservativeResize(c);
456 
457  //get selected MT
458  MatrixXd t_MT(p_outMT.rows(), idcs.rows()*3);
459 
460  for(qint32 j = 0; j < idcs.rows(); ++j)
461  t_MT.block(0, j*3, t_MT.rows(), 3) = p_outMT.block(0, (idcs[j]+offset)*3, t_MT.rows(), 3);
462 
463  qint32 nSens = t_MT.rows();
464  qint32 nSources = t_MT.cols()/3;
465 
466  if (nSources > 0)
467  {
468  RegionMT t_sensMT;
469 
470  t_sensMT.idcs = idcs;
471  t_sensMT.iLabelIdxIn = i;
472  t_sensMT.nClusters = ceil((double)nSources/(double)p_iClusterSize);
473 
474  t_sensMT.matRoiMTOrig = t_MT;
475 
476  printf("%d Cluster(s)... ", t_sensMT.nClusters);
477 
478  // Reshape Input data -> sources rows; sensors columns
479  t_sensMT.matRoiMT = MatrixXd(t_MT.cols()/3, 3*nSens);
480 
481  for(qint32 j = 0; j < nSens; ++j)
482  for(qint32 k = 0; k < t_sensMT.matRoiMT.rows(); ++k)
483  t_sensMT.matRoiMT.block(k,j*3,1,3) = t_MT.block(j,k*3,1,3);
484 
485  t_sensMT.sDistMeasure = p_sMethod;
486 
487  m_qListRegionMTIn.append(t_sensMT);
488 
489  printf("[added]\n");
490  }
491  else
492  {
493  printf("failed! Label contains no sources.\n");
494  }
495  }
496  }
497 
498  //
499  // Calculate clusters
500  //
501  printf("Clustering... ");
502  QFuture< RegionMTOut > res;
503  res = QtConcurrent::mapped(m_qListRegionMTIn, &RegionMT::cluster);
504  res.waitForFinished();
505 
506  //
507  // Assign results
508  //
509  MatrixXd t_MT_partial;
510 
511  qint32 nClusters;
512  qint32 nSens;
513  QList<RegionMT>::const_iterator itIn;
514  itIn = m_qListRegionMTIn.begin();
515  QFuture<RegionMTOut>::const_iterator itOut;
516  for (itOut = res.constBegin(); itOut != res.constEnd(); ++itOut)
517  {
518  nClusters = itOut->ctrs.rows();
519  nSens = itOut->ctrs.cols()/3;
520  t_MT_partial = MatrixXd::Zero(nSens, nClusters*3);
521 
522 // std::cout << "Number of Clusters: " << nClusters << " x " << nSens << std::endl;//itOut->iLabelIdcsOut << std::endl;
523 
524  //
525  // Assign the centroid for each cluster to the partial G
526  //
527  //ToDo change this use indeces found with whitened data
528  for(qint32 j = 0; j < nSens; ++j)
529  for(qint32 k = 0; k < nClusters; ++k)
530  t_MT_partial.block(j, k*3, 1, 3) = itOut->ctrs.block(k,j*3,1,3);
531 
532  //
533  // Get cluster indizes and its distances to the centroid
534  //
535  for(qint32 j = 0; j < nClusters; ++j)
536  {
537  VectorXi clusterIdcs = VectorXi::Zero(itOut->roiIdx.rows());
538  VectorXd clusterDistance = VectorXd::Zero(itOut->roiIdx.rows());
539  qint32 nClusterIdcs = 0;
540  for(qint32 k = 0; k < itOut->roiIdx.rows(); ++k)
541  {
542  if(itOut->roiIdx[k] == j)
543  {
544  clusterIdcs[nClusterIdcs] = itIn->idcs[k];
545 // qint32 offset = h == 0 ? 0 : this->src[0].nuse;
546 // Q_UNUSED(offset)
547  clusterDistance[nClusterIdcs] = itOut->D(k,j);
548  ++nClusterIdcs;
549  }
550  }
551  clusterIdcs.conservativeResize(nClusterIdcs);
552  clusterDistance.conservativeResize(nClusterIdcs);
553 
554  VectorXi clusterVertnos = VectorXi::Zero(clusterIdcs.size());
555  for(qint32 k = 0; k < clusterVertnos.size(); ++k)
556  clusterVertnos(k) = this->src[h].vertno[clusterIdcs(k)];
557 
558  t_qListMNEClusterInfo[h].clusterVertnos.append(clusterVertnos);
559 
560  }
561 
562  //
563  // Assign partial G to new LeadField
564  //
565  if(t_MT_partial.rows() > 0 && t_MT_partial.cols() > 0)
566  {
567  t_MT_new.conservativeResize(t_MT_partial.rows(), t_MT_new.cols() + t_MT_partial.cols());
568  t_MT_new.block(0, t_MT_new.cols() - t_MT_partial.cols(), t_MT_new.rows(), t_MT_partial.cols()) = t_MT_partial;
569 
570  // Map the centroids to the closest rr
571  for(qint32 k = 0; k < nClusters; ++k)
572  {
573  qint32 j = 0;
574 
575  double sqec = sqrt((itIn->matRoiMTOrig.block(0, j*3, itIn->matRoiMTOrig.rows(), 3) - t_MT_partial.block(0, k*3, t_MT_partial.rows(), 3)).array().pow(2).sum());
576  double sqec_min = sqec;
577  qint32 j_min = 0;
578 // MatrixXd matGainDiff;
579  for(qint32 j = 1; j < itIn->idcs.rows(); ++j)
580  {
581  sqec = sqrt((itIn->matRoiMTOrig.block(0, j*3, itIn->matRoiMTOrig.rows(), 3) - t_MT_partial.block(0, k*3, t_MT_partial.rows(), 3)).array().pow(2).sum());
582 
583  if(sqec < sqec_min)
584  {
585  sqec_min = sqec;
586  j_min = j;
587 // matGainDiff = itIn->matRoiGOrig.block(0, j*3, itIn->matRoiGOrig.rows(), 3) - t_G_partial.block(0, k*3, t_G_partial.rows(), 3);
588  }
589  }
590 
591 // qListGainDist.append(matGainDiff);
592 
593  // Take the closest coordinates
594 // qint32 sel_idx = itIn->idcs[j_min];
595 // Q_UNUSED(sel_idx)
596 
597 // //vertices
598 // std::cout << this->src[h].vertno[sel_idx] << ", ";
599 
600  ++count;
601  }
602  }
603 
604  ++itIn;
605  }
606 
607  printf("[done]\n");
608  }
609 
610  //
611  // Cluster operator D (sources x clusters)
612  //
613  qint32 totalNumOfClust = 0;
614  for (qint32 h = 0; h < 2; ++h)
615  totalNumOfClust += t_qListMNEClusterInfo[h].clusterVertnos.size();
616 
617  if(this->isFixedOrient())
618  p_D = MatrixXd::Zero(p_outMT.cols(), totalNumOfClust);
619  else
620  p_D = MatrixXd::Zero(p_outMT.cols(), totalNumOfClust*3);
621 
622  QList<VectorXi> t_vertnos = this->src.get_vertno();
623 
624 // qDebug() << "Size: " << t_vertnos[0].size() << t_vertnos[1].size();
625 // qDebug() << "this->sol->data.cols(): " << this->sol->data.cols();
626 
627  qint32 currentCluster = 0;
628  for (qint32 h = 0; h < 2; ++h)
629  {
630  int hemiOffset = h == 0 ? 0 : t_vertnos[0].size();
631  for(qint32 i = 0; i < t_qListMNEClusterInfo[h].clusterVertnos.size(); ++i)
632  {
633  VectorXi idx_sel;
634  MNEMath::intersect(t_vertnos[h], t_qListMNEClusterInfo[h].clusterVertnos[i], idx_sel);
635 
636 // std::cout << "\nVertnos:\n" << t_vertnos[h] << std::endl;
637 
638 // std::cout << "clusterVertnos[i]:\n" << t_qListMNEClusterInfo[h].clusterVertnos[i] << std::endl;
639 
640  idx_sel.array() += hemiOffset;
641 
642 // std::cout << "idx_sel]:\n" << idx_sel << std::endl;
643 
644  double selectWeight = 1.0/idx_sel.size();
645  if(this->isFixedOrient())
646  {
647  for(qint32 j = 0; j < idx_sel.size(); ++j)
648  p_D.col(currentCluster)[idx_sel(j)] = selectWeight;
649  }
650  else
651  {
652  qint32 clustOffset = currentCluster*3;
653  for(qint32 j = 0; j < idx_sel.size(); ++j)
654  {
655  qint32 idx_sel_Offset = idx_sel(j)*3;
656  //x
657  p_D(idx_sel_Offset,clustOffset) = selectWeight;
658  //y
659  p_D(idx_sel_Offset+1, clustOffset+1) = selectWeight;
660  //z
661  p_D(idx_sel_Offset+2, clustOffset+2) = selectWeight;
662  }
663  }
664  ++currentCluster;
665  }
666  }
667 
668 // 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;
669 
670  //
671  // Put it all together
672  //
673  p_outMT = t_MT_new;
674 
675  return p_outMT;
676 }
677 
678 //=============================================================================================================
679 
681  MNEForwardSolution forward,
682  const FiffCov &p_noise_cov,
683  float loose,
684  float depth,
685  bool fixed,
686  bool limit_depth_chs)
687 {
688  bool is_fixed_ori = forward.isFixedOrient();
689  MNEInverseOperator p_MNEInverseOperator;
690 
691  std::cout << "ToDo MNEInverseOperator::make_inverse_operator: do surf_ori check" << std::endl;
692 
693  //Check parameters
694  if(fixed && loose > 0)
695  {
696  qWarning("Warning: When invoking make_inverse_operator with fixed = true, the loose parameter is ignored.\n");
697  loose = 0.0f;
698  }
699 
700  if(is_fixed_ori && !fixed)
701  {
702  qWarning("Warning: Setting fixed parameter = true. Because the given forward operator has fixed orientation and can only be used to make a fixed-orientation inverse operator.\n");
703  fixed = true;
704  }
705 
706  if(forward.source_ori == -1 && loose > 0)
707  {
708  qCritical("Error: Forward solution is not oriented in surface coordinates. loose parameter should be 0 not %f.\n", loose);
709  return p_MNEInverseOperator;
710  }
711 
712  if(loose < 0 || loose > 1)
713  {
714  qWarning("Warning: Loose value should be in interval [0,1] not %f.\n", loose);
715  loose = loose > 1 ? 1 : 0;
716  printf("Setting loose to %f.\n", loose);
717  }
718 
719  if(depth < 0 || depth > 1)
720  {
721  qWarning("Warning: Depth value should be in interval [0,1] not %f.\n", depth);
722  depth = depth > 1 ? 1 : 0;
723  printf("Setting depth to %f.\n", depth);
724  }
725 
726  //
727  // 1. Read the bad channels
728  // 2. Read the necessary data from the forward solution matrix file
729  // 3. Load the projection data
730  // 4. Load the sensor noise covariance matrix and attach it to the forward
731  //
732  FiffInfo gain_info;
733  MatrixXd gain;
734  MatrixXd whitener;
735  qint32 n_nzero;
736  FiffCov p_outNoiseCov;
737  forward.prepare_forward(info, p_noise_cov, false, gain_info, gain, p_outNoiseCov, whitener, n_nzero);
738 
739  //
740  // 5. Compose the depth weight matrix
741  //
742  FiffCov::SDPtr p_depth_prior;
743  MatrixXd patch_areas;
744  if(depth > 0)
745  {
746  std::cout << "ToDo: patch_areas" << std::endl;
747 // patch_areas = forward.get('patch_areas', None)
748  p_depth_prior = FiffCov::SDPtr(new FiffCov(MNEForwardSolution::compute_depth_prior(gain, gain_info, is_fixed_ori, depth, 10.0, patch_areas, limit_depth_chs)));
749  }
750  else
751  {
752  p_depth_prior->data = MatrixXd::Ones(gain.cols(), gain.cols());
753  p_depth_prior->kind = FIFFV_MNE_DEPTH_PRIOR_COV;
754  p_depth_prior->diag = true;
755  p_depth_prior->dim = gain.cols();
756  p_depth_prior->nfree = 1;
757  }
758 
759  // Deal with fixed orientation forward / inverse
760  if(fixed)
761  {
762  if(depth < 0 || depth > 1)
763  {
764  // Convert the depth prior into a fixed-orientation one
765  printf("\tToDo: Picked elements from a free-orientation depth-weighting prior into the fixed-orientation one.\n");
766  }
767  if(!is_fixed_ori)
768  {
769  // Convert to the fixed orientation forward solution now
770  qint32 count = 0;
771  for(qint32 i = 2; i < p_depth_prior->data.rows(); i+=3)
772  {
773  p_depth_prior->data.row(count) = p_depth_prior->data.row(i);
774  ++count;
775  }
776  p_depth_prior->data.conservativeResize(count, 1);
777 
778 // forward = deepcopy(forward)
779  forward.to_fixed_ori();
780  is_fixed_ori = forward.isFixedOrient();
781  forward.prepare_forward(info, p_outNoiseCov, false, gain_info, gain, p_outNoiseCov, whitener, n_nzero);
782  }
783  }
784  printf("\tComputing inverse operator with %d channels.\n", gain_info.ch_names.size());
785 
786  //
787  // 6. Compose the source covariance matrix
788  //
789  printf("\tCreating the source covariance matrix\n");
790  FiffCov::SDPtr p_source_cov = p_depth_prior;
791 
792  // apply loose orientations
793  FiffCov::SDPtr p_orient_prior;
794  if(!is_fixed_ori)
795  {
796  p_orient_prior = FiffCov::SDPtr(new FiffCov(forward.compute_orient_prior(loose)));
797  p_source_cov->data.array() *= p_orient_prior->data.array();
798  }
799 
800  // 7. Apply fMRI weighting (not done)
801 
802  //
803  // 8. Apply the linear projection to the forward solution
804  // 9. Apply whitening to the forward computation matrix
805  //
806  printf("\tWhitening the forward solution.\n");
807  gain = whitener*gain;
808 
809  // 10. Exclude the source space points within the labels (not done)
810 
811  //
812  // 11. Do appropriate source weighting to the forward computation matrix
813  //
814 
815  // Adjusting Source Covariance matrix to make trace of G*R*G' equal
816  // to number of sensors.
817  printf("\tAdjusting source covariance matrix.\n");
818  RowVectorXd source_std = p_source_cov->data.array().sqrt().transpose();
819 
820  for(qint32 i = 0; i < gain.rows(); ++i)
821  gain.row(i) = gain.row(i).array() * source_std.array();
822 
823  double trace_GRGT = (gain * gain.transpose()).trace();//pow(gain.norm(), 2);
824  double scaling_source_cov = (double)n_nzero / trace_GRGT;
825 
826  p_source_cov->data.array() *= scaling_source_cov;
827 
828  gain.array() *= sqrt(scaling_source_cov);
829 
830  // now np.trace(np.dot(gain, gain.T)) == n_nzero
831  // logger.info(np.trace(np.dot(gain, gain.T)), n_nzero)
832 
833  //
834  // 12. Decompose the combined matrix
835  //
836  printf("Computing SVD of whitened and weighted lead field matrix.\n");
837  JacobiSVD<MatrixXd> svd(gain, ComputeThinU | ComputeThinV);
838  std::cout << "ToDo Sorting Necessary?" << std::endl;
839  VectorXd p_sing = svd.singularValues();
840  MatrixXd t_U = svd.matrixU();
841  MNEMath::sort<double>(p_sing, t_U);
842  FiffNamedMatrix::SDPtr p_eigen_fields = FiffNamedMatrix::SDPtr(new FiffNamedMatrix( svd.matrixU().cols(),
843  svd.matrixU().rows(),
844  defaultQStringList,
845  gain_info.ch_names,
846  t_U.transpose() ));
847 
848  p_sing = svd.singularValues();
849  MatrixXd t_V = svd.matrixV();
850  MNEMath::sort<double>(p_sing, t_V);
851  FiffNamedMatrix::SDPtr p_eigen_leads = FiffNamedMatrix::SDPtr(new FiffNamedMatrix( svd.matrixV().rows(),
852  svd.matrixV().cols(),
853  defaultQStringList,
854  defaultQStringList,
855  t_V ));
856  printf("\tlargest singular value = %f\n", p_sing.maxCoeff());
857  printf("\tscaling factor to adjust the trace = %f\n", trace_GRGT);
858 
859  qint32 p_nave = 1.0;
860 
861  // Handle methods
862  bool has_meg = false;
863  bool has_eeg = false;
864 
865  RowVectorXd ch_idx(info.chs.size());
866  qint32 count = 0;
867  for(qint32 i = 0; i < info.chs.size(); ++i)
868  {
869  if(gain_info.ch_names.contains(info.chs[i].ch_name))
870  {
871  ch_idx[count] = i;
872  ++count;
873  }
874  }
875  ch_idx.conservativeResize(count);
876 
877  for(qint32 i = 0; i < ch_idx.size(); ++i)
878  {
879  QString ch_type = info.channel_type(ch_idx[i]);
880  if (ch_type == "eeg")
881  has_eeg = true;
882  if ((ch_type == "mag") || (ch_type == "grad"))
883  has_meg = true;
884  }
885 
886  qint32 p_iMethods;
887 
888  if(has_eeg && has_meg)
889  p_iMethods = FIFFV_MNE_MEG_EEG;
890  else if(has_meg)
891  p_iMethods = FIFFV_MNE_MEG;
892  else
893  p_iMethods = FIFFV_MNE_EEG;
894 
895  // We set this for consistency with mne C code written inverses
896  if(depth == 0)
897  p_depth_prior = FiffCov::SDPtr();
898 
899  p_MNEInverseOperator.eigen_fields = p_eigen_fields;
900  p_MNEInverseOperator.eigen_leads = p_eigen_leads;
901  p_MNEInverseOperator.sing = p_sing;
902  p_MNEInverseOperator.nave = p_nave;
903  p_MNEInverseOperator.depth_prior = p_depth_prior;
904  p_MNEInverseOperator.source_cov = p_source_cov;
905  p_MNEInverseOperator.noise_cov = FiffCov::SDPtr(new FiffCov(p_outNoiseCov));
906  p_MNEInverseOperator.orient_prior = p_orient_prior;
907  p_MNEInverseOperator.projs = info.projs;
908  p_MNEInverseOperator.eigen_leads_weighted = false;
909  p_MNEInverseOperator.source_ori = forward.source_ori;
910  p_MNEInverseOperator.mri_head_t = forward.mri_head_t;
911  p_MNEInverseOperator.methods = p_iMethods;
912  p_MNEInverseOperator.nsource = forward.nsource;
913  p_MNEInverseOperator.coord_frame = forward.coord_frame;
914  p_MNEInverseOperator.source_nn = forward.source_nn;
915  p_MNEInverseOperator.src = forward.src;
916  p_MNEInverseOperator.info = forward.info;
917  p_MNEInverseOperator.info.bads = info.bads;
918 
919  return p_MNEInverseOperator;
920 }
921 
922 //=============================================================================================================
923 
924 MNEInverseOperator MNEInverseOperator::prepare_inverse_operator(qint32 nave ,float lambda2, bool dSPM, bool sLORETA) const
925 {
926  if(nave <= 0)
927  {
928  printf("The number of averages should be positive\n");
929  return MNEInverseOperator();
930  }
931  printf("Preparing the inverse operator for use...\n");
932  MNEInverseOperator inv(*this);
933  //
934  // Scale some of the stuff
935  //
936  float scale = ((float)inv.nave)/((float)nave);
937  inv.noise_cov->data *= scale;
938  inv.noise_cov->eig *= scale;
939  inv.source_cov->data *= scale;
940  //
941  if (inv.eigen_leads_weighted)
942  inv.eigen_leads->data *= sqrt(scale);
943  //
944  printf("\tScaled noise and source covariance from nave = %d to nave = %d\n",inv.nave,nave);
945  inv.nave = nave;
946  //
947  // Create the diagonal matrix for computing the regularized inverse
948  //
949  VectorXd tmp = inv.sing.cwiseProduct(inv.sing) + VectorXd::Constant(inv.sing.size(), lambda2);
950 // if(inv.reginv)
951 // delete inv.reginv;
952  inv.reginv = VectorXd(inv.sing.cwiseQuotient(tmp));
953  printf("\tCreated the regularized inverter\n");
954  //
955  // Create the projection operator
956  //
957 
958  qint32 ncomp = FiffProj::make_projector(inv.projs, inv.noise_cov->names, inv.proj);
959  if (ncomp > 0)
960  printf("\tCreated an SSP operator (subspace dimension = %d)\n",ncomp);
961 
962  //
963  // Create the whitener
964  //
965 // if(inv.whitener)
966 // delete inv.whitener;
967  inv.whitener = MatrixXd::Zero(inv.noise_cov->dim, inv.noise_cov->dim);
968 
969  qint32 nnzero, k;
970  if (inv.noise_cov->diag == 0)
971  {
972  //
973  // Omit the zeroes due to projection
974  //
975  nnzero = 0;
976 
977  for (k = ncomp; k < inv.noise_cov->dim; ++k)
978  {
979  if (inv.noise_cov->eig[k] > 0)
980  {
981  inv.whitener(k,k) = 1.0/sqrt(inv.noise_cov->eig[k]);
982  ++nnzero;
983  }
984  }
985  //
986  // Rows of eigvec are the eigenvectors
987  //
988  inv.whitener *= inv.noise_cov->eigvec;
989  printf("\tCreated the whitener using a full noise covariance matrix (%d small eigenvalues omitted)\n", inv.noise_cov->dim - nnzero);
990  }
991  else
992  {
993  //
994  // No need to omit the zeroes due to projection
995  //
996  for (k = 0; k < inv.noise_cov->dim; ++k)
997  inv.whitener(k,k) = 1.0/sqrt(inv.noise_cov->data(k,0));
998 
999  printf("\tCreated the whitener using a diagonal noise covariance matrix (%d small eigenvalues discarded)\n",ncomp);
1000  }
1001  //
1002  // Finally, compute the noise-normalization factors
1003  //
1004  if (dSPM || sLORETA)
1005  {
1006  VectorXd noise_norm = VectorXd::Zero(inv.eigen_leads->nrow);
1007  VectorXd noise_weight;
1008  if (dSPM)
1009  {
1010  printf("\tComputing noise-normalization factors (dSPM)...");
1011  noise_weight = VectorXd(inv.reginv);
1012  }
1013  else
1014  {
1015  printf("\tComputing noise-normalization factors (sLORETA)...");
1016  VectorXd tmp = (VectorXd::Constant(inv.sing.size(), 1) + inv.sing.cwiseProduct(inv.sing)/lambda2);
1017  noise_weight = inv.reginv.cwiseProduct(tmp.cwiseSqrt());
1018  }
1019  VectorXd one;
1020  if (inv.eigen_leads_weighted)
1021  {
1022  for (k = 0; k < inv.eigen_leads->nrow; ++k)
1023  {
1024  one = inv.eigen_leads->data.block(k,0,1,inv.eigen_leads->data.cols()).cwiseProduct(noise_weight);
1025  noise_norm[k] = sqrt(one.dot(one));
1026  }
1027  }
1028  else
1029  {
1030 // qDebug() << 32;
1031  double c;
1032  for (k = 0; k < inv.eigen_leads->nrow; ++k)
1033  {
1034 // qDebug() << 321;
1035  c = sqrt(inv.source_cov->data(k,0));
1036 // qDebug() << 322;
1037 // qDebug() << "inv.eigen_leads->data" << inv.eigen_leads->data.rows() << "x" << inv.eigen_leads->data.cols();
1038 // qDebug() << "noise_weight" << noise_weight.rows() << "x" << noise_weight.cols();
1039  one = c*(inv.eigen_leads->data.row(k).transpose()).cwiseProduct(noise_weight);//ToDo eigenleads data -> pointer
1040  noise_norm[k] = sqrt(one.dot(one));
1041 // qDebug() << 324;
1042  }
1043  }
1044 
1045 // qDebug() << 4;
1046 
1047  //
1048  // Compute the final result
1049  //
1050  VectorXd noise_norm_new;
1051  if (inv.source_ori == FIFFV_MNE_FREE_ORI)
1052  {
1053  //
1054  // The three-component case is a little bit more involved
1055  // The variances at three consequtive entries must be squeared and
1056  // added together
1057  //
1058  // Even in this case return only one noise-normalization factor
1059  // per source location
1060  //
1061  VectorXd* t = MNEMath::combine_xyz(noise_norm.transpose());
1062  noise_norm_new = t->cwiseSqrt();//double otherwise values are getting too small
1063  delete t;
1064  //
1065  // This would replicate the same value on three consequtive
1066  // entries
1067  //
1068  // noise_norm = kron(sqrt(mne_combine_xyz(noise_norm)),ones(3,1));
1069  }
1070  VectorXd vOnes = VectorXd::Ones(noise_norm_new.size());
1071  VectorXd tmp = vOnes.cwiseQuotient(noise_norm_new.cwiseAbs());
1072 // if(inv.noisenorm)
1073 // delete inv.noisenorm;
1074 
1075  typedef Eigen::Triplet<double> T;
1076  std::vector<T> tripletList;
1077  tripletList.reserve(noise_norm_new.size());
1078  for(qint32 i = 0; i < noise_norm_new.size(); ++i)
1079  tripletList.push_back(T(i, i, tmp[i]));
1080 
1081  inv.noisenorm = SparseMatrix<double>(noise_norm_new.size(),noise_norm_new.size());
1082  inv.noisenorm.setFromTriplets(tripletList.begin(), tripletList.end());
1083 
1084  printf("[done]\n");
1085  }
1086  else
1087  {
1088 // if(inv.noisenorm)
1089 // delete inv.noisenorm;
1090  inv.noisenorm = SparseMatrix<double>();
1091  }
1092 
1093  return inv;
1094 }
1095 
1096 //=============================================================================================================
1097 
1099 {
1100  //
1101  // Open the file, create directory
1102  //
1103  FiffStream::SPtr t_pStream(new FiffStream(&p_IODevice));
1104  printf("Reading inverse operator decomposition from %s...\n",t_pStream->streamName().toUtf8().constData());
1105 
1106  if(!t_pStream->open())
1107  return false;
1108  //
1109  // Find all inverse operators
1110  //
1111  QList <FiffDirNode::SPtr> invs_list = t_pStream->dirtree()->dir_tree_find(FIFFB_MNE_INVERSE_SOLUTION);
1112  if ( invs_list.size()== 0)
1113  {
1114  printf("No inverse solutions in %s\n", t_pStream->streamName().toUtf8().constData());
1115  return false;
1116  }
1117  FiffDirNode::SPtr invs = invs_list[0];
1118  //
1119  // Parent MRI data
1120  //
1121  QList <FiffDirNode::SPtr> parent_mri = t_pStream->dirtree()->dir_tree_find(FIFFB_MNE_PARENT_MRI_FILE);
1122  if (parent_mri.size() == 0)
1123  {
1124  printf("No parent MRI information in %s", t_pStream->streamName().toUtf8().constData());
1125  return false;
1126  }
1127  printf("\tReading inverse operator info...");
1128  //
1129  // Methods and source orientations
1130  //
1131  FiffTag::SPtr t_pTag;
1132  if (!invs->find_tag(t_pStream, FIFF_MNE_INCLUDED_METHODS, t_pTag))
1133  {
1134  printf("Modalities not found\n");
1135  return false;
1136  }
1137 
1138  inv = MNEInverseOperator();
1139  inv.methods = *t_pTag->toInt();
1140  //
1141  if (!invs->find_tag(t_pStream, FIFF_MNE_SOURCE_ORIENTATION, t_pTag))
1142  {
1143  printf("Source orientation constraints not found\n");
1144  return false;
1145  }
1146  inv.source_ori = *t_pTag->toInt();
1147  //
1148  if (!invs->find_tag(t_pStream, FIFF_MNE_SOURCE_SPACE_NPOINTS, t_pTag))
1149  {
1150  printf("Number of sources not found\n");
1151  return false;
1152  }
1153  inv.nsource = *t_pTag->toInt();
1154  inv.nchan = 0;
1155  //
1156  // Coordinate frame
1157  //
1158  if (!invs->find_tag(t_pStream, FIFF_MNE_COORD_FRAME, t_pTag))
1159  {
1160  printf("Coordinate frame tag not found\n");
1161  return false;
1162  }
1163  inv.coord_frame = *t_pTag->toInt();
1164  //
1165  // The actual source orientation vectors
1166  //
1167  if (!invs->find_tag(t_pStream, FIFF_MNE_INVERSE_SOURCE_ORIENTATIONS, t_pTag))
1168  {
1169  printf("Source orientation information not found\n");
1170  return false;
1171  }
1172 
1173 // if(inv.source_nn)
1174 // delete inv.source_nn;
1175  inv.source_nn = t_pTag->toFloatMatrix();
1176  inv.source_nn.transposeInPlace();
1177 
1178  printf("[done]\n");
1179  //
1180  // The SVD decomposition...
1181  //
1182  printf("\tReading inverse operator decomposition...");
1183  if (!invs->find_tag(t_pStream, FIFF_MNE_INVERSE_SING, t_pTag))
1184  {
1185  printf("Singular values not found\n");
1186  return false;
1187  }
1188 
1189 // if(inv.sing)
1190 // delete inv.sing;
1191  inv.sing = Map<VectorXf>(t_pTag->toFloat(), t_pTag->size()/4).cast<double>();
1192  inv.nchan = inv.sing.rows();
1193  //
1194  // The eigenleads and eigenfields
1195  //
1196  inv.eigen_leads_weighted = false;
1197  if(!t_pStream->read_named_matrix(invs, FIFF_MNE_INVERSE_LEADS, *inv.eigen_leads.data()))
1198  {
1199  inv.eigen_leads_weighted = true;
1200  if(!t_pStream->read_named_matrix(invs, FIFF_MNE_INVERSE_LEADS_WEIGHTED, *inv.eigen_leads.data()))
1201  {
1202  printf("Error reading eigenleads named matrix.\n");
1203  return false;
1204  }
1205  }
1206  //
1207  // Having the eigenleads as columns is better for the inverse calculations
1208  //
1209  inv.eigen_leads->transpose_named_matrix();
1210 
1211  if(!t_pStream->read_named_matrix(invs, FIFF_MNE_INVERSE_FIELDS, *inv.eigen_fields.data()))
1212  {
1213  printf("Error reading eigenfields named matrix.\n");
1214  return false;
1215  }
1216  printf("[done]\n");
1217  //
1218  // Read the covariance matrices
1219  //
1220  if(t_pStream->read_cov(invs, FIFFV_MNE_NOISE_COV, *inv.noise_cov.data()))
1221  {
1222  printf("\tNoise covariance matrix read.\n");
1223  }
1224  else
1225  {
1226  printf("\tError: Not able to read noise covariance matrix.\n");
1227  return false;
1228  }
1229 
1230  if(t_pStream->read_cov(invs, FIFFV_MNE_SOURCE_COV, *inv.source_cov.data()))
1231  {
1232  printf("\tSource covariance matrix read.\n");
1233  }
1234  else
1235  {
1236  printf("\tError: Not able to read source covariance matrix.\n");
1237  return false;
1238  }
1239  //
1240  // Read the various priors
1241  //
1242  if(t_pStream->read_cov(invs, FIFFV_MNE_ORIENT_PRIOR_COV, *inv.orient_prior.data()))
1243  {
1244  printf("\tOrientation priors read.\n");
1245  }
1246  else
1247  inv.orient_prior->clear();
1248 
1249  if(t_pStream->read_cov(invs, FIFFV_MNE_DEPTH_PRIOR_COV, *inv.depth_prior.data()))
1250  {
1251  printf("\tDepth priors read.\n");
1252  }
1253  else
1254  {
1255  inv.depth_prior->clear();
1256  }
1257  if(t_pStream->read_cov(invs, FIFFV_MNE_FMRI_PRIOR_COV, *inv.fmri_prior.data()))
1258  {
1259  printf("\tfMRI priors read.\n");
1260  }
1261  else
1262  {
1263  inv.fmri_prior->clear();
1264  }
1265  //
1266  // Read the source spaces
1267  //
1268  if(!MNESourceSpace::readFromStream(t_pStream, false, inv.src))
1269  {
1270  printf("\tError: Could not read the source spaces.\n");
1271  return false;
1272  }
1273  for (qint32 k = 0; k < inv.src.size(); ++k)
1274  inv.src[k].id = MNESourceSpace::find_source_space_hemi(inv.src[k]);
1275  //
1276  // Get the MRI <-> head coordinate transformation
1277  //
1278  FiffCoordTrans mri_head_t;// = NULL;
1279  if (!parent_mri[0]->find_tag(t_pStream, FIFF_COORD_TRANS, t_pTag))
1280  {
1281  printf("MRI/head coordinate transformation not found\n");
1282  return false;
1283  }
1284  else
1285  {
1286  mri_head_t = t_pTag->toCoordTrans();
1287  if (mri_head_t.from != FIFFV_COORD_MRI || mri_head_t.to != FIFFV_COORD_HEAD)
1288  {
1289  mri_head_t.invert_transform();
1290  if (mri_head_t.from != FIFFV_COORD_MRI || mri_head_t.to != FIFFV_COORD_HEAD)
1291  {
1292  printf("MRI/head coordinate transformation not found");
1293 // if(mri_head_t)
1294 // delete mri_head_t;
1295  return false;
1296  }
1297  }
1298  }
1299  inv.mri_head_t = mri_head_t;
1300 
1301  //
1302  // get parent MEG info
1303  //
1304  t_pStream->read_meas_info_base(t_pStream->dirtree(), inv.info);
1305 
1306  //
1307  // Transform the source spaces to the correct coordinate frame
1308  // if necessary
1309  //
1310  if (inv.coord_frame != FIFFV_COORD_MRI && inv.coord_frame != FIFFV_COORD_HEAD)
1311  printf("Only inverse solutions computed in MRI or head coordinates are acceptable");
1312  //
1313  // Number of averages is initially one
1314  //
1315  inv.nave = 1;
1316  //
1317  // We also need the SSP operator
1318  //
1319  inv.projs = t_pStream->read_proj(t_pStream->dirtree());
1320  //
1321  // Some empty fields to be filled in later
1322  //
1323 // inv.proj = []; % This is the projector to apply to the data
1324 // inv.whitener = []; % This whitens the data
1325 // inv.reginv = []; % This the diagonal matrix implementing
1326 // % regularization and the inverse
1327 // inv.noisenorm = []; % These are the noise-normalization factors
1328  //
1329  if(!inv.src.transform_source_space_to(inv.coord_frame, mri_head_t))
1330  {
1331  printf("Could not transform source space.\n");
1332  }
1333  printf("\tSource spaces transformed to the inverse solution coordinate frame\n");
1334  //
1335  // Done!
1336  //
1337 
1338  return true;
1339 }
1340 
1341 //=============================================================================================================
1342 
1343 void MNEInverseOperator::write(QIODevice &p_IODevice)
1344 {
1345  //
1346  // Open the file, create directory
1347  //
1348 
1349  // Create the file and save the essentials
1350  FiffStream::SPtr t_pStream = FiffStream::start_file(p_IODevice);
1351  printf("Write inverse operator decomposition in %s...", t_pStream->streamName().toUtf8().constData());
1352  this->writeToStream(t_pStream.data());
1353  t_pStream->end_file();
1354 }
1355 
1356 //=============================================================================================================
1357 
1359 {
1360  p_pStream->start_block(FIFFB_MNE_INVERSE_SOLUTION);
1361 
1362  printf("\tWriting inverse operator info...\n");
1363 
1364  p_pStream->write_int(FIFF_MNE_INCLUDED_METHODS, &this->methods);
1365  p_pStream->write_int(FIFF_MNE_SOURCE_ORIENTATION, &this->source_ori);
1366  p_pStream->write_int(FIFF_MNE_SOURCE_SPACE_NPOINTS, &this->nsource);
1367  p_pStream->write_int(FIFF_MNE_COORD_FRAME, &this->coord_frame);
1369  VectorXf tmp_sing = this->sing.cast<float>();
1370  p_pStream->write_float(FIFF_MNE_INVERSE_SING, tmp_sing.data(), tmp_sing.size());
1371 
1372  //
1373  // The eigenleads and eigenfields
1374  //
1375  if(this->eigen_leads_weighted)
1376  {
1377  FiffNamedMatrix tmpMatrix(*this->eigen_leads.data());
1378  tmpMatrix.transpose_named_matrix();
1379  p_pStream->write_named_matrix(FIFF_MNE_INVERSE_LEADS_WEIGHTED, tmpMatrix);
1380  }
1381  else
1382  {
1383  FiffNamedMatrix tmpMatrix(*this->eigen_leads.data());
1384  tmpMatrix.transpose_named_matrix();
1385  p_pStream->write_named_matrix(FIFF_MNE_INVERSE_LEADS, tmpMatrix);
1386  }
1387 
1388  p_pStream->write_named_matrix(FIFF_MNE_INVERSE_FIELDS, *this->eigen_fields.data());
1389  printf("\t[done]\n");
1390  //
1391  // write the covariance matrices
1392  //
1393  printf("\tWriting noise covariance matrix.");
1394  p_pStream->write_cov(*this->noise_cov.data());
1395 
1396  printf("\tWriting source covariance matrix.\n");
1397  p_pStream->write_cov(*this->source_cov.data());
1398  //
1399  // write the various priors
1400  //
1401  printf("\tWriting orientation priors.\n");
1402  if(!this->orient_prior->isEmpty())
1403  p_pStream->write_cov(*this->orient_prior.data());
1404  if(!this->depth_prior->isEmpty())
1405  p_pStream->write_cov(*this->depth_prior.data());
1406  if(!this->fmri_prior->isEmpty())
1407  p_pStream->write_cov(*this->fmri_prior.data());
1408 
1409  //
1410  // Parent MRI data
1411  //
1412  p_pStream->start_block(FIFFB_MNE_PARENT_MRI_FILE);
1413  // write the MRI <-> head coordinate transformation
1414  p_pStream->write_coord_trans(this->mri_head_t);
1415  p_pStream->end_block(FIFFB_MNE_PARENT_MRI_FILE);
1416 
1417  //
1418  // Parent MEG measurement info
1419  //
1420  p_pStream->write_info_base(this->info);
1421 
1422  //
1423  // Write the source spaces
1424  //
1425  if(!src.isEmpty())
1426  this->src.writeToStream(p_pStream);
1427 
1428  //
1429  // We also need the SSP operator
1430  //
1431  p_pStream->write_proj(this->projs);
1432  //
1433  // Done!
1434  //
1435  p_pStream->end_block(FIFFB_MNE_INVERSE_SOLUTION);
1436 }
QSharedPointer< FiffStream > SPtr
Definition: fiff_stream.h:107
#define FIFFV_MNE_NOISE_COV
FIFFLIB::fiff_int_t coord_frame
Eigen::MatrixXd matRoiMT
FIFFLIB::FiffCov compute_orient_prior(float loose=0.2)
FIFFLIB::FiffCoordTrans mri_head_t
QStringList struct_names
Definition: colortable.h:112
#define FIFF_MNE_INVERSE_FIELDS
bool assemble_kernel(const FSLIB::Label &label, QString method, bool pick_normal, Eigen::MatrixXd &K, Eigen::SparseMatrix< double > &noise_norm, QList< Eigen::VectorXi > &vertno)
void write(QIODevice &p_IODevice)
fiff_long_t write_proj(const QList< FiffProj > &projs)
FIFFLIB::FiffInfoBase info
Vertices label based lookup table.
Definition: colortable.h:67
FIFFLIB::fiff_int_t source_ori
covariance data
Definition: fiff_cov.h:77
QSharedPointer< FiffDirNode > SPtr
Definition: fiff_dir_node.h:77
Eigen::MatrixXd matRoiMTOrig
fiff_long_t end_block(fiff_int_t kind, fiff_int_t next=FIFFV_NEXT_SEQ)
QList< FiffChInfo > chs
QString channel_type(qint32 idx) const
cluster information
QList< FIFFLIB::FiffProj > projs
#define FIFF_MNE_SOURCE_ORIENTATION
Annotation set.
Definition: annotationset.h:80
Coordinate transformation description.
QSharedDataPointer< FiffNamedMatrix > SDPtr
#define FIFFV_MNE_ORIENT_PRIOR_COV
FIFF measurement file information.
Definition: fiff_info.h:84
Freesurfer/MNE label.
Definition: label.h:80
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
#define FIFFV_MNE_DEPTH_PRIOR_COV
FIFFLIB::FiffCov::SDPtr noise_cov
void writeToStream(FIFFLIB::FiffStream *p_pStream)
FIFFLIB::FiffCov::SDPtr fmri_prior
Eigen::SparseMatrix< double > noisenorm
void writeToStream(FIFFLIB::FiffStream *p_pStream)
MNEInverseOperator prepare_inverse_operator(qint32 nave, float lambda2, bool dSPM, bool sLORETA=false) const
QList< FiffProj > projs
Definition: fiff_info.h:268
FIFF File I/O routines.
Definition: fiff_stream.h:104
QList< Eigen::VectorXi > label_src_vertno_sel(const FSLIB::Label &p_label, Eigen::VectorXi &src_sel) const
static bool read_inverse_operator(QIODevice &p_IODevice, MNEInverseOperator &inv)
QList< Eigen::VectorXi > get_vertno() const
#define FIFF_MNE_INVERSE_LEADS
Eigen::VectorXi idcs
#define FIFF_COORD_TRANS
Definition: fiff_file.h:475
FIFFLIB::FiffCov::SDPtr depth_prior
FIFFLIB::FiffCov::SDPtr source_cov
Eigen::VectorXi getLabelIds() const
Definition: colortable.h:120
bool transform_source_space_to(FIFFLIB::fiff_int_t dest, FIFFLIB::FiffCoordTrans &trans)
fiff_long_t write_int(fiff_int_t kind, const fiff_int_t *data, fiff_int_t nel=1, fiff_int_t next=FIFFV_NEXT_SEQ)
QSharedPointer< FiffTag > SPtr
Definition: fiff_tag.h:152
fiff_long_t write_float_matrix(fiff_int_t kind, const Eigen::MatrixXf &mat)
#define FIFF_MNE_INVERSE_SING
fiff_long_t write_cov(const FiffCov &p_FiffCov)
Label class declaration.
fiff_long_t write_info_base(const FiffInfoBase &p_FiffInfoBase)
QSharedDataPointer< FiffCov > SDPtr
Definition: fiff_cov.h:82
Eigen::MatrixXd cluster_kernel(const FSLIB::AnnotationSet &p_AnnotationSet, qint32 p_iClusterSize, Eigen::MatrixXd &p_D, QString p_sMethod="cityblock") const
FIFFLIB::FiffCoordTrans mri_head_t
fiff_long_t write_float(fiff_int_t kind, const float *data, fiff_int_t nel=1)
bool isEmpty() const
Definition: label.h:184
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)
bool check_ch_names(const FIFFLIB::FiffInfo &info) const
#define FIFF_MNE_SOURCE_SPACE_NPOINTS
fiff_long_t start_block(fiff_int_t kind)
FIFFLIB::FiffCov::SDPtr orient_prior
FIFFLIB::FiffNamedMatrix::SDPtr eigen_fields
#define FIFF_MNE_INVERSE_SOURCE_ORIENTATIONS
#define FIFF_MNE_COORD_FRAME
fiff_long_t write_coord_trans(const FiffCoordTrans &trans)
FIFFLIB::FiffNamedMatrix::SDPtr eigen_leads
static bool readFromStream(FIFFLIB::FiffStream::SPtr &p_pStream, bool add_geom, MNESourceSpace &p_SourceSpace)
fiff_long_t write_named_matrix(fiff_int_t kind, const FiffNamedMatrix &mat)
static qint32 find_source_space_hemi(MNEHemisphere &p_Hemisphere)
static MNEInverseOperator make_inverse_operator(const FIFFLIB::FiffInfo &info, MNEForwardSolution forward, const FIFFLIB::FiffCov &p_noise_cov, float loose=0.2f, float depth=0.8f, bool fixed=false, bool limit_depth_chs=true)
#define FIFF_MNE_INVERSE_LEADS_WEIGHTED