MNE-CPP  beta 1.0
rapmusic.h
Go to the documentation of this file.
1 //=============================================================================================================
35 #ifndef RAPMUSIC_H
36 #define RAPMUSIC_H
37 
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 
56 //*************************************************************************************************************
57 //=============================================================================================================
58 // EIGEN INCLUDES
59 //=============================================================================================================
60 
61 #include <Eigen/Core>
62 #include <Eigen/SVD>
63 #include <Eigen/LU>
64 
65 
66 //*************************************************************************************************************
67 //=============================================================================================================
68 // DEFINE NAMESPACE INVERSELIB
69 //=============================================================================================================
70 
71 namespace INVERSELIB
72 {
73 
74 //*************************************************************************************************************
75 //=============================================================================================================
76 // USED NAMESPACES
77 //=============================================================================================================
78 
79 using namespace MNELIB;
80 
81 
82 //*************************************************************************************************************
83 //=============================================================================================================
84 // SOME DEFINES
85 //=============================================================================================================
86 
87 #define NOT_TRANSPOSED 0
88 #define IS_TRANSPOSED 1
91 //=============================================================================================================
92 
95 typedef struct Pair
96 {
97  int x1;
98  int x2;
99 } Pair;
100 
101 
102 
103 //=============================================================================================================
110 {
111 public:
112  typedef QSharedPointer<RapMusic> SPtr;
113  typedef QSharedPointer<const RapMusic> ConstSPtr;
115  //*********************************************************************************************************
116  //=========================================================================================================
117  // TYPEDEFS
118  //=========================================================================================================
119 
120  typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> MatrixXT;
122  typedef Eigen::Matrix<double, Eigen::Dynamic, 6> MatrixX6T;
124  typedef Eigen::Matrix<double, 6, Eigen::Dynamic> Matrix6XT;
126  typedef Eigen::Matrix<double, 6, 6> Matrix6T;
128  typedef Eigen::Matrix<double, Eigen::Dynamic, 1> VectorXT;
130  typedef Eigen::Matrix<double, 6, 1> Vector6T;
134  //=========================================================================================================
138  RapMusic();
139 
140  //=========================================================================================================
150  RapMusic(MNEForwardSolution& p_pFwd, bool p_bSparsed, int p_iN = 2, double p_dThr = 0.5);
151 
152  virtual ~RapMusic();
153 
154  //=========================================================================================================
165  bool init(MNEForwardSolution& p_pFwd, bool p_bSparsed = false, int p_iN = 2, double p_dThr = 0.5);
166 
167  virtual MNESourceEstimate calculateInverse(const FiffEvoked &p_fiffEvoked, bool pick_normal = false);
168 
169  virtual MNESourceEstimate calculateInverse(const MatrixXd &data, float tmin, float tstep) const;
170 
171  virtual MNESourceEstimate calculateInverse(const MatrixXd& p_matMeasurement, QList< DipolePair<double> > &p_RapDipoles) const;
172 
173  virtual const char* getName() const;
174 
175  virtual const MNESourceSpace& getSourceSpace() const;
176 
177  //=========================================================================================================
184  void setStcAttr(int p_iSampStcWin, float p_fStcOverlap);
185 
186 protected:
187  //=========================================================================================================
196  int calcPhi_s(const MatrixXT& p_matMeasurement, MatrixXT* &p_pMatPhi_s) const;
197 
198  //=========================================================================================================
211  static double subcorr(MatrixX6T& p_matProj_G, const MatrixXT& p_pMatU_B);
212 
213  //=========================================================================================================
229  static double subcorr(MatrixX6T& p_matProj_G, const MatrixXT& p_matU_B, Vector6T& p_vec_phi_k_1);
230 
231  //=========================================================================================================
241  static void calcA_k_1( const MatrixX6T& p_matG_k_1,
242  const Vector6T& p_matPhi_k_1,
243  const int p_iIdxk_1,
244  MatrixXT& p_matA_k_1);
245 
246  //=========================================================================================================
253  void calcOrthProj(const MatrixXT& p_matA_k_1, MatrixXT& p_matOrthProj) const;
254 
255  //=========================================================================================================
266  void calcPairCombinations( const int p_iNumPoints,
267  const int p_iNumCombinations,
268  Pair** p_ppPairIdxCombinations) const;
269 
270  //=========================================================================================================
286  static void getPointPair(const int p_iPoints, const int p_iCurIdx, int &p_iIdx1, int &p_iIdx2);
287 
288  //=========================================================================================================
297  static void getGainMatrixPair( const MatrixXT& p_matGainMarix,
298  MatrixX6T& p_matGainMarix_Pair,
299  int p_iIdx1, int p_iIdx2);
300 
301  //=========================================================================================================
311  static void insertSource( int p_iDipoleIdx1, int p_iDipoleIdx2,
312  const Vector6T &p_vec_phi_k_1,
313  double p_valCor,
314  QList< DipolePair<double> > &p_RapDipoles);
315 
318  int m_iN;
319  double m_dThreshold;
331  bool m_bIsInit;
333  //Stc stuff
337  //=========================================================================================================
345  static inline int getRank(const MatrixXT& p_matSigma);
346 
347  //=========================================================================================================
357  static inline int useFullRank( const MatrixXT& p_Mat,
358  const MatrixXT& p_matSigma_src,
359  MatrixXT& p_matFull_Rank,
360  int type = NOT_TRANSPOSED);
361 
362  //=========================================================================================================
369  static inline MatrixXT makeSquareMat(const MatrixXT& p_matF);
370 };
371 
372 //*************************************************************************************************************
373 //=============================================================================================================
374 // INLINE DEFINITIONS
375 //=============================================================================================================
376 
377 inline int RapMusic::getRank(const MatrixXT& p_matSigma)
378 {
379  int t_iRank;
380  //if once a singularvalue is smaller than epsilon = 10^-5 the following values are also smaller
381  // -> because Singular values are ordered
382  for(t_iRank = p_matSigma.rows()-1; t_iRank > 0; t_iRank--)
383  if (p_matSigma(t_iRank, t_iRank) > 0.00001)
384  break;
385 
386  t_iRank++;//rank corresponding to epsilon
387 
388  return t_iRank;
389 }
390 
391 
392 //*************************************************************************************************************
393 
394 inline int RapMusic::useFullRank( const MatrixXT& p_Mat,
395  const MatrixXT& p_matSigma_src,
396  MatrixXT& p_matFull_Rank,
397  int type)
398 {
399  int rank = getRank(p_matSigma_src);
400 
401  if (type == NOT_TRANSPOSED)
402  p_matFull_Rank = p_Mat.block(0,0,p_Mat.rows(),rank);
403  else
404  p_matFull_Rank = p_Mat.block(0,0,rank,p_Mat.cols());
405 
406  return rank;
407 }
408 
409 
410 //*************************************************************************************************************
411 
413 {
414  //Make rectangular - p_matF*p_matF^T
415  //MatrixXT FFT = p_matF*p_matF.transpose();
416 
417  return p_matF*p_matF.transpose();
418 }
419 
420 } //NAMESPACE
421 
422 #endif // RAPMUSIC_H
Eigen::Matrix< double, 6, 1 > Vector6T
Definition: rapmusic.h:130
static int useFullRank(const MatrixXT &p_Mat, const MatrixXT &p_matSigma_src, MatrixXT &p_matFull_Rank, int type=NOT_TRANSPOSED)
Definition: rapmusic.h:394
Source Space descritpion.
#define NOT_TRANSPOSED
Definition: rapmusic.h:87
Eigen::Matrix< double, Eigen::Dynamic, 6 > MatrixX6T
Definition: rapmusic.h:122
QSharedPointer< RapMusic > SPtr
Definition: rapmusic.h:112
QSharedPointer< const RapMusic > ConstSPtr
Definition: rapmusic.h:113
Eigen::Matrix< double, 6, Eigen::Dynamic > Matrix6XT
Definition: rapmusic.h:124
int m_iNumLeadFieldCombinations
Definition: rapmusic.h:325
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixXT
Definition: rapmusic.h:120
evoked data
Definition: fiff_evoked.h:91
Eigen::Matrix< double, 6, 6 > Matrix6T
Definition: rapmusic.h:126
ToDo Documentation...
MNEForwardSolution class declaration, which provides the forward solution including the source space ...
#define INVERSESHARED_EXPORT
MNEForwardSolution m_ForwardSolution
Definition: rapmusic.h:316
static MatrixXT makeSquareMat(const MatrixXT &p_matF)
Definition: rapmusic.h:412
static int getRank(const MatrixXT &p_matSigma)
Definition: rapmusic.h:377
MNESourceEstimate class declaration.
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorXT
Definition: rapmusic.h:128
The RapMusic class provides the RAP MUSIC Algorithm CPU implementation. ToDo: Paper references...
Definition: rapmusic.h:109
Inverse algorithm interface.
Pair ** m_ppPairIdxCombinations
Definition: rapmusic.h:327