inv_rap_music.h

Go to the documentation of this file.

//=============================================================================================================
 
#ifndef INV_RAP_MUSIC_H
#define INV_RAP_MUSIC_H
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "../inv_global.h"
 
#include "inv_dipole.h"
 
#include <mne/mne_forward_solution.h>
#include <inv/inv_source_estimate.h>
#include <time.h>
 
#include <QVector>
 
#include <vector>
 
//=============================================================================================================
// EIGEN INCLUDES
//=============================================================================================================
 
#include <Eigen/Core>
#include <Eigen/SVD>
#include <Eigen/LU>
 
//=============================================================================================================
// DEFINE NAMESPACE INVLIB
//=============================================================================================================
 
namespace INVLIB
{
 
//=============================================================================================================
// SOME DEFINES
//=============================================================================================================
 
#define NOT_TRANSPOSED 0 
#define IS_TRANSPOSED 1  
 
//=============================================================================================================
struct Pair
{
    int x1; 
    int x2; 
};
 
//=============================================================================================================
class INVSHARED_EXPORT InvRapMusic
{
public:
    typedef QSharedPointer<InvRapMusic> SPtr;            
    typedef QSharedPointer<const InvRapMusic> ConstSPtr; 
 
    //*********************************************************************************************************
    //=========================================================================================================
    // TYPEDEFS
    //=========================================================================================================
 
    typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> MatrixXT; 
    typedef Eigen::Matrix<double, Eigen::Dynamic, 6> MatrixX6T;             
    typedef Eigen::Matrix<double, 6, Eigen::Dynamic> Matrix6XT;             
    typedef Eigen::Matrix<double, 6, 6> Matrix6T;                           
    typedef Eigen::Matrix<double, Eigen::Dynamic, 1> VectorXT;              
    typedef Eigen::Matrix<double, 6, 1> Vector6T;                           
 
    //=========================================================================================================
    InvRapMusic();
 
    //=========================================================================================================
    InvRapMusic(MNELIB::MNEForwardSolution& p_pFwd, bool p_bSparsed, int p_iN = 2, double p_dThr = 0.5);
 
    virtual ~InvRapMusic();
 
    //=========================================================================================================
    bool init(MNELIB::MNEForwardSolution& p_pFwd, bool p_bSparsed = false, int p_iN = 2, double p_dThr = 0.5);
 
    virtual InvSourceEstimate calculateInverse(const FIFFLIB::FiffEvoked& p_fiffEvoked, bool pick_normal = false);
 
    virtual InvSourceEstimate calculateInverse(const Eigen::MatrixXd& data, float tmin, float tstep, bool pick_normal = false) const;
 
    virtual InvSourceEstimate calculateInverse(const Eigen::MatrixXd& p_matMeasurement, QList<InvDipolePair<double>>& p_RapDipoles) const;
 
    virtual const char* getName() const;
 
    virtual const MNELIB::MNESourceSpaces& getSourceSpace() const;
 
    //=========================================================================================================
    void setStcAttr(int p_iSampStcWin, float p_fStcOverlap);
 
protected:
    //=========================================================================================================
    int calcPhi_s(const MatrixXT& p_matMeasurement, MatrixXT*& p_pMatPhi_s) const;
 
    //=========================================================================================================
    static double subcorr(MatrixX6T& p_matProj_G, const MatrixXT& p_pMatU_B);
 
    //=========================================================================================================
    static double subcorr(MatrixX6T& p_matProj_G, const MatrixXT& p_matU_B, Vector6T& p_vec_phi_k_1);
 
    //=========================================================================================================
    static void calcA_k_1(const MatrixX6T& p_matG_k_1,
                          const Vector6T& p_matPhi_k_1,
                          const int p_iIdxk_1,
                          MatrixXT& p_matA_k_1);
 
    //=========================================================================================================
    void calcOrthProj(const MatrixXT& p_matA_k_1, MatrixXT& p_matOrthProj) const;
 
    //=========================================================================================================
    void calcPairCombinations(const int p_iNumPoints,
                              const int p_iNumCombinations,
                              std::vector<Pair>& p_pairIdxCombinations) const;
 
    //=========================================================================================================
    static void getPointPair(const int p_iPoints, const int p_iCurIdx, int& p_iIdx1, int& p_iIdx2);
 
    //=========================================================================================================
    static void getGainMatrixPair(const MatrixXT& p_matGainMarix,
                                  MatrixX6T& p_matGainMarix_Pair,
                                  int p_iIdx1, int p_iIdx2);
 
    //=========================================================================================================
    static void insertSource(int p_iDipoleIdx1, int p_iDipoleIdx2,
                             const Vector6T& p_vec_phi_k_1,
                             double p_valCor,
                             QList<InvDipolePair<double>>& p_RapDipoles);
 
    MNELIB::MNEForwardSolution m_ForwardSolution; 
 
    int m_iN;            
    double m_dThreshold; 
 
    int m_iNumGridPoints;            
    int m_iNumChannels;              
    int m_iNumLeadFieldCombinations; 
 
    std::vector<Pair> m_ppPairIdxCombinations; 
 
    int m_iMaxNumThreads; 
 
    bool m_bIsInit; 
 
    //Stc stuff
    int m_iSamplesStcWindow; 
    float m_fStcOverlap;     
 
    //=========================================================================================================
    static inline int getRank(const MatrixXT& p_matSigma);
 
    //=========================================================================================================
    static inline int useFullRank(const MatrixXT& p_Mat,
                                  const MatrixXT& p_matSigma_src,
                                  MatrixXT& p_matFull_Rank,
                                  int type = NOT_TRANSPOSED);
 
    //=========================================================================================================
    static inline MatrixXT makeSquareMat(const MatrixXT& p_matF);
};
 
//=============================================================================================================
// INLINE DEFINITIONS
//=============================================================================================================
 
inline int InvRapMusic::getRank(const MatrixXT& p_matSigma)
{
    int t_iRank;
    //if once a singularvalue is smaller than epsilon = 10^-5 the following values are also smaller
    // -> because Singular values are ordered
    for (t_iRank = p_matSigma.rows() - 1; t_iRank > 0; t_iRank--)
        if (p_matSigma(t_iRank, t_iRank) > 0.00001)
            break;
 
    t_iRank++; //rank corresponding to epsilon
 
    return t_iRank;
}
 
//=============================================================================================================
 
inline int InvRapMusic::useFullRank(const MatrixXT& p_Mat,
                                    const MatrixXT& p_matSigma_src,
                                    MatrixXT& p_matFull_Rank,
                                    int type)
{
    int rank = getRank(p_matSigma_src);
 
    if (type == NOT_TRANSPOSED)
        p_matFull_Rank = p_Mat.block(0, 0, p_Mat.rows(), rank);
    else
        p_matFull_Rank = p_Mat.block(0, 0, rank, p_Mat.cols());
 
    return rank;
}
 
//=============================================================================================================
 
inline InvRapMusic::MatrixXT InvRapMusic::makeSquareMat(const MatrixXT& p_matF)
{
    //Make rectangular - p_matF*p_matF^T
    //MatrixXT FFT = p_matF*p_matF.transpose();
 
    MatrixXT mat = p_matF.transpose();
 
    return p_matF * mat;
}
} //NAMESPACE
 
#endif // INV_RAP_MUSIC_H