rapmusic.h

Go to the documentation of this file.

//=============================================================================================================
#ifndef RAPMUSIC_H
#define RAPMUSIC_H

//=============================================================================================================
// INCLUDES
//=============================================================================================================

#include "../inverse_global.h"
#include "../IInverseAlgorithm.h"

#include "dipole.h"

#include <mne/mne_forwardsolution.h>
#include <mne/mne_sourceestimate.h>
#include <time.h>

#include <QVector>

//=============================================================================================================
// EIGEN INCLUDES
//=============================================================================================================

#include <Eigen/Core>
#include <Eigen/SVD>
#include <Eigen/LU>

//=============================================================================================================
// DEFINE NAMESPACE INVERSELIB
//=============================================================================================================

namespace INVERSELIB
{

//=============================================================================================================
// SOME DEFINES
//=============================================================================================================

#define NOT_TRANSPOSED   0  
#define IS_TRANSPOSED   1   
//=============================================================================================================

typedef struct Pair
{
    int x1; 
    int x2; 
} Pair;

//=============================================================================================================
class INVERSESHARED_EXPORT RapMusic : public IInverseAlgorithm
{
public:
    typedef QSharedPointer<RapMusic> SPtr;             
    typedef QSharedPointer<const RapMusic> 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;                            
    //=========================================================================================================
    RapMusic();

    //=========================================================================================================
    RapMusic(MNELIB::MNEForwardSolution& p_pFwd, bool p_bSparsed, int p_iN = 2, double p_dThr = 0.5);

    virtual ~RapMusic();

    //=========================================================================================================
    bool init(MNELIB::MNEForwardSolution& p_pFwd, bool p_bSparsed = false, int p_iN = 2, double p_dThr = 0.5);

    virtual MNELIB::MNESourceEstimate calculateInverse(const FIFFLIB::FiffEvoked &p_fiffEvoked, bool pick_normal = false);

    virtual MNELIB::MNESourceEstimate calculateInverse(const Eigen::MatrixXd &data, float tmin, float tstep, bool pick_normal = false) const;

    virtual MNELIB::MNESourceEstimate calculateInverse(const Eigen::MatrixXd& p_matMeasurement, QList< DipolePair<double> > &p_RapDipoles) const;

    virtual const char* getName() const;

    virtual const MNELIB::MNESourceSpace& 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< DipolePair<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 RapMusic::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 RapMusic::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 RapMusic::MatrixXT RapMusic::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 // RAPMUSIC_H