MNE-CPP  0.1.9
A Framework for Electrophysiology
rapmusic.h
Go to the documentation of this file.
1 //=============================================================================================================
36 #ifndef RAPMUSIC_H
37 #define RAPMUSIC_H
38 
39 //=============================================================================================================
40 // INCLUDES
41 //=============================================================================================================
42 
43 #include "../inverse_global.h"
44 #include "../IInverseAlgorithm.h"
45 
46 #include "dipole.h"
47 
49 #include <mne/mne_sourceestimate.h>
50 #include <time.h>
51 
52 #include <QVector>
53 
54 //=============================================================================================================
55 // EIGEN INCLUDES
56 //=============================================================================================================
57 
58 #include <Eigen/Core>
59 #include <Eigen/SVD>
60 #include <Eigen/LU>
61 
62 //=============================================================================================================
63 // DEFINE NAMESPACE INVERSELIB
64 //=============================================================================================================
65 
66 namespace INVERSELIB
67 {
68 
69 //=============================================================================================================
70 // SOME DEFINES
71 //=============================================================================================================
72 
73 #define NOT_TRANSPOSED 0
74 #define IS_TRANSPOSED 1
76 //=============================================================================================================
77 
80 typedef struct Pair
81 {
82  int x1;
83  int x2;
84 } Pair;
85 
86 //=============================================================================================================
93 {
94 public:
95  typedef QSharedPointer<RapMusic> SPtr;
96  typedef QSharedPointer<const RapMusic> ConstSPtr;
98  //*********************************************************************************************************
99  //=========================================================================================================
100  // TYPEDEFS
101  //=========================================================================================================
102 
103  typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> MatrixXT;
105  typedef Eigen::Matrix<double, Eigen::Dynamic, 6> MatrixX6T;
107  typedef Eigen::Matrix<double, 6, Eigen::Dynamic> Matrix6XT;
109  typedef Eigen::Matrix<double, 6, 6> Matrix6T;
111  typedef Eigen::Matrix<double, Eigen::Dynamic, 1> VectorXT;
113  typedef Eigen::Matrix<double, 6, 1> Vector6T;
116  //=========================================================================================================
120  RapMusic();
121 
122  //=========================================================================================================
132  RapMusic(MNELIB::MNEForwardSolution& p_pFwd, bool p_bSparsed, int p_iN = 2, double p_dThr = 0.5);
133 
134  virtual ~RapMusic();
135 
136  //=========================================================================================================
147  bool init(MNELIB::MNEForwardSolution& p_pFwd, bool p_bSparsed = false, int p_iN = 2, double p_dThr = 0.5);
148 
149  virtual MNELIB::MNESourceEstimate calculateInverse(const FIFFLIB::FiffEvoked &p_fiffEvoked, bool pick_normal = false);
150 
151  virtual MNELIB::MNESourceEstimate calculateInverse(const Eigen::MatrixXd &data, float tmin, float tstep, bool pick_normal = false) const;
152 
153  virtual MNELIB::MNESourceEstimate calculateInverse(const Eigen::MatrixXd& p_matMeasurement, QList< DipolePair<double> > &p_RapDipoles) const;
154 
155  virtual const char* getName() const;
156 
157  virtual const MNELIB::MNESourceSpace& getSourceSpace() const;
158 
159  //=========================================================================================================
166  void setStcAttr(int p_iSampStcWin, float p_fStcOverlap);
167 
168 protected:
169  //=========================================================================================================
178  int calcPhi_s(const MatrixXT& p_matMeasurement, MatrixXT* &p_pMatPhi_s) const;
179 
180  //=========================================================================================================
193  static double subcorr(MatrixX6T& p_matProj_G, const MatrixXT& p_pMatU_B);
194 
195  //=========================================================================================================
211  static double subcorr(MatrixX6T& p_matProj_G, const MatrixXT& p_matU_B, Vector6T& p_vec_phi_k_1);
212 
213  //=========================================================================================================
223  static void calcA_k_1( const MatrixX6T& p_matG_k_1,
224  const Vector6T& p_matPhi_k_1,
225  const int p_iIdxk_1,
226  MatrixXT& p_matA_k_1);
227 
228  //=========================================================================================================
235  void calcOrthProj(const MatrixXT& p_matA_k_1, MatrixXT& p_matOrthProj) const;
236 
237  //=========================================================================================================
248  void calcPairCombinations( const int p_iNumPoints,
249  const int p_iNumCombinations,
250  Pair** p_ppPairIdxCombinations) const;
251 
252  //=========================================================================================================
268  static void getPointPair(const int p_iPoints, const int p_iCurIdx, int &p_iIdx1, int &p_iIdx2);
269 
270  //=========================================================================================================
279  static void getGainMatrixPair( const MatrixXT& p_matGainMarix,
280  MatrixX6T& p_matGainMarix_Pair,
281  int p_iIdx1, int p_iIdx2);
282 
283  //=========================================================================================================
293  static void insertSource( int p_iDipoleIdx1, int p_iDipoleIdx2,
294  const Vector6T &p_vec_phi_k_1,
295  double p_valCor,
296  QList< DipolePair<double> > &p_RapDipoles);
297 
300  int m_iN;
301  double m_dThreshold;
313  bool m_bIsInit;
315  //Stc stuff
319  //=========================================================================================================
327  static inline int getRank(const MatrixXT& p_matSigma);
328 
329  //=========================================================================================================
339  static inline int useFullRank( const MatrixXT& p_Mat,
340  const MatrixXT& p_matSigma_src,
341  MatrixXT& p_matFull_Rank,
342  int type = NOT_TRANSPOSED);
343 
344  //=========================================================================================================
351  static inline MatrixXT makeSquareMat(const MatrixXT& p_matF);
352 };
353 
354 //=============================================================================================================
355 // INLINE DEFINITIONS
356 //=============================================================================================================
357 
358 inline int RapMusic::getRank(const MatrixXT& p_matSigma)
359 {
360  int t_iRank;
361  //if once a singularvalue is smaller than epsilon = 10^-5 the following values are also smaller
362  // -> because Singular values are ordered
363  for(t_iRank = p_matSigma.rows()-1; t_iRank > 0; t_iRank--)
364  if (p_matSigma(t_iRank, t_iRank) > 0.00001)
365  break;
366 
367  t_iRank++;//rank corresponding to epsilon
368 
369  return t_iRank;
370 }
371 
372 //=============================================================================================================
373 
374 inline int RapMusic::useFullRank( const MatrixXT& p_Mat,
375  const MatrixXT& p_matSigma_src,
376  MatrixXT& p_matFull_Rank,
377  int type)
378 {
379  int rank = getRank(p_matSigma_src);
380 
381  if (type == NOT_TRANSPOSED)
382  p_matFull_Rank = p_Mat.block(0,0,p_Mat.rows(),rank);
383  else
384  p_matFull_Rank = p_Mat.block(0,0,rank,p_Mat.cols());
385 
386  return rank;
387 }
388 
389 //=============================================================================================================
390 
392 {
393  //Make rectangular - p_matF*p_matF^T
394  //MatrixXT FFT = p_matF*p_matF.transpose();
395 
396  MatrixXT mat = p_matF.transpose();
397 
398  return p_matF*mat;
399 }
400 } //NAMESPACE
401 
402 #endif // RAPMUSIC_H
Eigen::Matrix< double, 6, 6 > Matrix6T
Definition: rapmusic.h:109
static int useFullRank(const MatrixXT &p_Mat, const MatrixXT &p_matSigma_src, MatrixXT &p_matFull_Rank, int type=NOT_TRANSPOSED)
Definition: rapmusic.h:374
#define INVERSESHARED_EXPORT
Eigen::Matrix< double, 6, Eigen::Dynamic > Matrix6XT
Definition: rapmusic.h:107
MNESourceEstimate class declaration.
QSharedPointer< RapMusic > SPtr
Definition: rapmusic.h:95
int m_iNumLeadFieldCombinations
Definition: rapmusic.h:307
ToDo Documentation...
The RapMusic class provides the RAP MUSIC Algorithm CPU implementation. ToDo: Paper references...
Definition: rapmusic.h:92
static int getRank(const MatrixXT &p_matSigma)
Definition: rapmusic.h:358
Source Space descritpion.
evoked data
Definition: fiff_evoked.h:77
Inverse algorithm interface.
#define NOT_TRANSPOSED
Definition: rapmusic.h:73
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixXT
Definition: rapmusic.h:103
MNEForwardSolution class declaration, which provides the forward solution including the source space ...
static MatrixXT makeSquareMat(const MatrixXT &p_matF)
Definition: rapmusic.h:391
Eigen::Matrix< double, Eigen::Dynamic, 6 > MatrixX6T
Definition: rapmusic.h:105
Eigen::Matrix< double, 6, 1 > Vector6T
Definition: rapmusic.h:113
MNELIB::MNEForwardSolution m_ForwardSolution
Definition: rapmusic.h:298
QSharedPointer< const RapMusic > ConstSPtr
Definition: rapmusic.h:96
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorXT
Definition: rapmusic.h:111
Pair ** m_ppPairIdxCombinations
Definition: rapmusic.h:309