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>
 
//=============================================================================================================
// 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   
 
//=============================================================================================================
typedef struct Pair
{
    int x1; 
    int x2; 
} Pair;
 
//=============================================================================================================
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,
                                Pair** p_ppPairIdxCombinations) 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;    
 
    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