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