mne_forwardsolution.h

Go to the documentation of this file.

//=============================================================================================================
 
#ifndef MNE_FORWARDSOLUTION_H
#define MNE_FORWARDSOLUTION_H
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "mne_global.h"
#include "mne_source_spaces.h"
 
#include <utils/mnemath.h>
#include <utils/kmeans.h>
 
#include <fs/annotationset.h>
 
#include <fiff/fiff_constants.h>
#include <fiff/fiff_coord_trans.h>
#include <fiff/fiff_types.h>
#include <fiff/fiff_info_base.h>
#include <fiff/fiff_cov.h>
 
#include <math.h>
 
//=============================================================================================================
// EIGEN INCLUDES
//=============================================================================================================
 
#include <Eigen/Core>
 
//=============================================================================================================
// QT INCLUDES
//=============================================================================================================
 
#include <QFile>
#include <QSharedPointer>
#include <QDataStream>
 
//=============================================================================================================
// DEFINE NAMESPACE MNELIB
//=============================================================================================================
 
namespace MNELIB
{
 
//=========================================================================================================
struct RegionDataOut
{
    Eigen::VectorXi roiIdx;         
    Eigen::MatrixXd ctrs;           
    Eigen::VectorXd sumd;           
    Eigen::MatrixXd D;              
 
    qint32 iLabelIdxOut;            
};
 
//=========================================================================================================
struct RegionData
{
    Eigen::MatrixXd matRoiG;            
    Eigen::MatrixXd matRoiGWhitened;    
    bool bUseWhitened;                  
 
    Eigen::MatrixXd matRoiGOrig;        
//    Eigen::MatrixXd    matRoiGOrigWhitened;    /**< Whitened region gain matrix sensors x sources(x,y,z)*/
 
    qint32 nClusters;      
 
    Eigen::VectorXi idcs;           
    qint32 iLabelIdxIn;    
    QString sDistMeasure;   
 
    RegionDataOut cluster() const
    {
        QString t_sDistMeasure;
        if(sDistMeasure.isEmpty())
            t_sDistMeasure = QString("cityblock");
        else
            t_sDistMeasure = sDistMeasure;
 
        // Kmeans Reduction
        RegionDataOut p_RegionDataOut;
 
        UTILSLIB::KMeans t_kMeans(t_sDistMeasure, QString("sample"), 5);
 
        if(bUseWhitened)
        {
            t_kMeans.calculate(this->matRoiGWhitened, this->nClusters, p_RegionDataOut.roiIdx, p_RegionDataOut.ctrs, p_RegionDataOut.sumd, p_RegionDataOut.D);
 
            Eigen::MatrixXd newCtrs = Eigen::MatrixXd::Zero(p_RegionDataOut.ctrs.rows(), p_RegionDataOut.ctrs.cols());
            for(qint32 c = 0; c < p_RegionDataOut.ctrs.rows(); ++c)
            {
                qint32 num = 0;
 
                for(qint32 idx = 0; idx < p_RegionDataOut.roiIdx.size(); ++idx)
                {
                    if(c == p_RegionDataOut.roiIdx[idx])
                    {
                        newCtrs.row(c) += this->matRoiG.row(idx); //just take whitened to get indeces calculate centroids using the original matrix
                        ++num;
                    }
                }
 
                if(num > 0)
                    newCtrs.row(c) /= num;
            }
            p_RegionDataOut.ctrs = newCtrs; //Replace whitened with original
        }
        else
            t_kMeans.calculate(this->matRoiG, this->nClusters, p_RegionDataOut.roiIdx, p_RegionDataOut.ctrs, p_RegionDataOut.sumd, p_RegionDataOut.D);
 
        p_RegionDataOut.iLabelIdxOut = this->iLabelIdxIn;
 
        return p_RegionDataOut;
    }
};
 
const static FIFFLIB::FiffCov defaultCov;
const static FIFFLIB::FiffInfo defaultInfo;
static Eigen::MatrixXd defaultD;
 
//=============================================================================================================
class MNESHARED_EXPORT MNEForwardSolution
{
public:
    typedef QSharedPointer<MNEForwardSolution> SPtr;            
    typedef QSharedPointer<const MNEForwardSolution> ConstSPtr; 
 
    //=========================================================================================================
    MNEForwardSolution();
 
    //=========================================================================================================
    MNEForwardSolution(QIODevice &p_IODevice,
                       bool force_fixed = false,
                       bool surf_ori = false,
                       const QStringList& include = FIFFLIB::defaultQStringList,
                       const QStringList& exclude = FIFFLIB::defaultQStringList,
                       bool bExcludeBads = false);
 
    //=========================================================================================================
    MNEForwardSolution(const MNEForwardSolution &p_MNEForwardSolution);
 
    //=========================================================================================================
    ~MNEForwardSolution();
 
    //=========================================================================================================
    void clear();
 
    //=========================================================================================================
    MNEForwardSolution cluster_forward_solution(const FSLIB::AnnotationSet &p_AnnotationSet,
                                                qint32 p_iClusterSize,
                                                Eigen::MatrixXd& p_D = defaultD,
                                                const FIFFLIB::FiffCov &p_pNoise_cov = defaultCov,
                                                const FIFFLIB::FiffInfo &p_pInfo = defaultInfo,
                                                QString p_sMethod = "cityblock") const;
 
    //=========================================================================================================
    FIFFLIB::FiffCov compute_orient_prior(float loose = 0.2);
 
    //=========================================================================================================
    static FIFFLIB::FiffCov compute_depth_prior(const Eigen::MatrixXd &Gain,
                                                const FIFFLIB::FiffInfo &gain_info,
                                                bool is_fixed_ori,
                                                double exp = 0.8,
                                                double limit = 10.0,
                                                const Eigen::MatrixXd &patch_areas = FIFFLIB::defaultConstMatrixXd,
                                                bool limit_depth_chs = false);
 
    //=========================================================================================================
    inline bool isClustered() const;
 
    //=========================================================================================================
    inline bool isEmpty() const;
 
    //=========================================================================================================
    inline bool isFixedOrient() const;
 
    //=========================================================================================================
    MNEForwardSolution pick_channels(const QStringList& include = FIFFLIB::defaultQStringList,
                                     const QStringList& exclude = FIFFLIB::defaultQStringList) const;
 
    //=========================================================================================================
    MNEForwardSolution pick_regions(const QList<FSLIB::Label> &p_qListLabels) const;
 
    //=========================================================================================================
    MNEForwardSolution pick_types(bool meg,
                                  bool eeg,
                                  const QStringList& include = FIFFLIB::defaultQStringList,
                                  const QStringList& exclude = FIFFLIB::defaultQStringList) const;
 
    //=========================================================================================================
    void prepare_forward(const FIFFLIB::FiffInfo &p_info,
                         const FIFFLIB::FiffCov &p_noise_cov,
                         bool p_pca,
                         FIFFLIB::FiffInfo &p_outFwdInfo,
                         Eigen::MatrixXd &gain,
                         FIFFLIB::FiffCov &p_outNoiseCov,
                         Eigen::MatrixXd &p_outWhitener,
                         qint32 &p_outNumNonZero) const;
 
//    //=========================================================================================================
//    /**
//    * Prepares a forward solution, Bad channels, after clustering etc ToDo...
//    *
//    * @param[in] p_FiffInfo   Fif measurement info
//    *
//    */
//    void prepare_forward(const FiffInfo& p_FiffInfo)
//    {
//        QStringList fwdChNames = this->sol->row_names;
//        QStringList chNames;
//        for(qint32 i = 0; i < p_FiffInfo.ch_names.size(); ++i)
//        {
//            bool inBads = false;
//            bool inFwd = false;
 
//            for(qint32 j = 0; j < p_FiffInfo.bads.size(); ++j)
//            {
//                if(QString::compare(p_FiffInfo.bads[j], p_FiffInfo.ch_names[i]) == 0)
//                {
//                    inBads = true;
//                    break;
//                }
//            }
 
//            for(qint32 j = 0; j < fwdChNames.size(); ++j)
//            {
//                if(QString::compare(fwdChNames[j], p_FiffInfo.ch_names[i]) == 0)
//                {
//                    inFwd = true;
//                    break;
//                }
//            }
 
//            if(!inBads && inFwd)
//                chNames.append(p_FiffInfo.ch_names[i]);
//        }
 
//        qint32 nchan = chNames.size();
//    }
 
    //=========================================================================================================
    Eigen::VectorXi tripletSelection(const Eigen::VectorXi& p_vecIdxSelection) const
    {
        Eigen::MatrixXi triSelect = p_vecIdxSelection.transpose().replicate(3,1).array() * 3;//repmat((p_vecIdxSelection - 1) * 3 + 1, 3, 1);
        triSelect.row(1).array() += 1;
        triSelect.row(2).array() += 2;
        Eigen::VectorXi retTriSelect(triSelect.cols()*3);
        for(int i = 0; i < triSelect.cols(); ++i)
            retTriSelect.block(i*3,0,3,1) = triSelect.col(i);
        return retTriSelect;
    } // tripletSelection
 
    //=========================================================================================================
    static bool read(QIODevice& p_IODevice,
                     MNEForwardSolution& fwd,
                     bool force_fixed = false,
                     bool surf_ori = false,
                     const QStringList& include = FIFFLIB::defaultQStringList,
                     const QStringList& exclude = FIFFLIB::defaultQStringList,
                     bool bExcludeBads = true);
 
    //ToDo readFromStream
 
    //=========================================================================================================
    MNEForwardSolution reduce_forward_solution(qint32 p_iNumDipoles, Eigen::MatrixXd& p_D) const;
 
    //=========================================================================================================
    static void restrict_gain_matrix(Eigen::MatrixXd &G, const FIFFLIB::FiffInfo &info);
 
    //=========================================================================================================
    void to_fixed_ori();
 
    //=========================================================================================================
    friend std::ostream& operator<<(std::ostream& out, const MNELIB::MNEForwardSolution &p_MNEForwardSolution);

    friend bool operator== (const MNEForwardSolution &a, const MNEForwardSolution &b);
 
    //=========================================================================================================
    Eigen::MatrixX3f getSourcePositionsByLabel(const QList<FSLIB::Label> &lPickedLabels,
                                               const FSLIB::SurfaceSet& tSurfSetInflated);
 
private:
    //=========================================================================================================
    static bool read_one(FIFFLIB::FiffStream::SPtr& p_pStream,
                         const FIFFLIB::FiffDirNode::SPtr& p_Node,
                         MNEForwardSolution& one);
 
public:
    FIFFLIB::FiffInfoBase info;                 
    FIFFLIB::fiff_int_t source_ori;             
    bool surf_ori;                              
    FIFFLIB::fiff_int_t coord_frame;            
    FIFFLIB::fiff_int_t nsource;                
    FIFFLIB::fiff_int_t nchan;                  
    FIFFLIB::FiffNamedMatrix::SDPtr sol;        
    FIFFLIB::FiffNamedMatrix::SDPtr sol_grad;   
    FIFFLIB::FiffCoordTrans mri_head_t;         
    MNESourceSpaces src;                         
    Eigen::MatrixX3f source_rr;                 
    Eigen::MatrixX3f source_nn;                 
};
 
//=============================================================================================================
// INLINE DEFINITIONS
//=============================================================================================================
 
inline bool MNEForwardSolution::isClustered() const
{
    auto* hemi = src.hemisphereAt(0);
    return hemi && hemi->isClustered();
}
 
//=============================================================================================================
 
inline bool MNEForwardSolution::isEmpty() const
{
    return this->nchan <= 0;
}
 
//=============================================================================================================
 
inline bool MNEForwardSolution::isFixedOrient() const
{
    return this->source_ori == FIFFV_MNE_FIXED_ORI;
}
 
//=============================================================================================================
 
inline std::ostream& operator<<(std::ostream& out, const MNELIB::MNEForwardSolution &p_MNEForwardSolution)
{
    out << "#### MNE Forward Solution ####\n";
 
    out << "\n source_ori: " << p_MNEForwardSolution.source_ori << std::endl;
    out << "\n coord_frame: " << p_MNEForwardSolution.coord_frame << std::endl;
    out << "\n nsource: " << p_MNEForwardSolution.nsource << std::endl;
    out << "\n nchan: " << p_MNEForwardSolution.nchan << std::endl;
    out << "\n sol:\n\t" << *p_MNEForwardSolution.sol.data() << std::endl;
    out << "\n sol_grad:\n\t" << *p_MNEForwardSolution.sol_grad.data() << std::endl;
 
    return out;
}
 
//=============================================================================================================
 
inline bool operator== (const MNEForwardSolution &a, const MNEForwardSolution &b)
{
    return (a.info == b.info &&
            a.source_ori == b.source_ori &&
            a.surf_ori == b.surf_ori &&
            a.coord_frame == b.coord_frame &&
            a.nsource == b.nsource &&
            a.nchan == b.nchan &&
            *a.sol == *b.sol &&
            *a.sol_grad == *b.sol_grad &&
            a.mri_head_t == b.mri_head_t &&
            a.src == b.src &&
            a.source_rr.isApprox(b.source_rr, 0.0001f) &&
            a.source_nn.isApprox(b.source_nn, 0.0001f));
}
} // NAMESPACE
 
#endif // MNE_FORWARDSOLUTION_H