mne_forward_solution.h

Go to the documentation of this file.

//=============================================================================================================
 
#ifndef MNE_FORWARD_SOLUTION_H
#define MNE_FORWARD_SOLUTION_H
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "mne_global.h"
 
#include <mne/mne_source_spaces.h>
 
#include <math/kmeans.h>
 
#include <fs/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>
#include <memory>
 
//=============================================================================================================
// EIGEN INCLUDES
//=============================================================================================================
 
#include <Eigen/Core>
 
//=============================================================================================================
// QT INCLUDES
//=============================================================================================================
 
#include <QSharedPointer>
 
//=============================================================================================================
// 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; 
 
    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);
 
            // Centroids of the original gain: all channels, not only those the whitener kept
            Eigen::MatrixXd newCtrs = Eigen::MatrixXd::Zero(p_RegionDataOut.ctrs.rows(), this->matRoiG.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; 
    typedef std::unique_ptr<MNEForwardSolution> UPtr;           
 
    //=========================================================================================================
    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& operator=(const MNEForwardSolution& other);
 
    //=========================================================================================================
    MNEForwardSolution(MNEForwardSolution&& other) = default;
 
    //=========================================================================================================
    MNEForwardSolution& operator=(MNEForwardSolution&& other) = default;
 
    //=========================================================================================================
    ~MNEForwardSolution();
 
    //=========================================================================================================
    void clear();
 
    //=========================================================================================================
    MNEForwardSolution cluster_forward_solution(const FSLIB::FsAnnotationSet& 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);
 
    //=========================================================================================================
    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::FsLabel>& 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;
 
    //=========================================================================================================
    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 = false);
 
    //=========================================================================================================
    bool write(QIODevice& p_IODevice) const;
 
    //=========================================================================================================
    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();
 
    //=========================================================================================================
    void convert_to_surf_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::FsLabel>& lPickedLabels,
                                               const FSLIB::FsSurfaceSet& 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;       
    QString mri_filename;                     
    FIFFLIB::FiffId mri_id;                   
    MNELIB::MNESourceSpaces src;              
    Eigen::MatrixX3f source_rr;               
    Eigen::MatrixX3f source_nn;               
};
 
//=============================================================================================================
// INLINE DEFINITIONS
//=============================================================================================================
 
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_FORWARD_SOLUTION_H