ex_connectivity

Example of using the MNE-CPP Connectivity library.

Libraries: Connectivity, Inverse, 3D Visualization, Visualization, Signal Processing, FreeSurfer, MNE Core, FIFF I/O, Utilities

Options

Option	Description
--fileIn	The input file <in>.
--eve	Path to the event <file>.
--fwd	Path to forwad solution <file>.
--subject	Selected subject <subject>.
--subjectPath	Selected subject path <subjectPath>.
--cov	Path to the covariance <file>.
--annotType	FsAnnotation <type> (for source level usage only).
--doSourceLoc	Do source localization (for source level usage only).
--doClust	Do clustering of source space (for source level usage only).
--sourceLocMethod	Inverse estimation <method> (for source level usage only), i.e., 'MNE', 'dSPM' or 'sLORETA'.
--connectMethod	Connectivity <method>, i.e., 'COR', 'XCOR.
--snr	The SNR <value> used for computation (for source level usage only).
--aveIdx	The average <index> to choose from the average file.
--chType	The channel <type> (for sensor level usage only), i.e. 'eeg' or 'meg'.
--coilType	The coil <type> (for sensor level usage only), i.e. 'grad' or 'mag'.
--tmin	The time minimum value for averaging in seconds relativ to the trigger onset.
--tmax	The time maximum value for averaging in seconds relativ to the trigger onset.

Build and run

cmake --build build/developer-dynamic --target ex_connectivity
./out/Release/examples/ex_connectivity --help

The default inputs come from the MNE sample dataset.

Used in the API reference

• AbstractMetric
• Connectivity
• ConnectivitySettings
• Network
• VisualizationInfo
• ConnectivitySettingsView
• ImageSc
• Plot
• TFplot
• RtConnectivity
• InvMinimumNorm
• MNE
• MNEEpochDataList
• MNEForwardSolution

Source

src/examples/ex_connectivity/connectivitysettingsmanager.h

//=============================================================================================================
/**
 * SPDX-License-Identifier: BSD-3-Clause
 * Copyright (c) 2018-2026 MNE-CPP Authors
 *
 * @file     connectivitysettingsmanager.h
 * @author   Gabriel Motta <gabrielbenmotta@gmail.com>;
 *           johaenns <j.vorw01@gmail.com>;
 *           Christoph Dinh <christoph.dinh@mne-cpp.org>;
 *           Lorenz Esch <lorenz.esch@tu-ilmenau.de>
 * @since    0.1.0
 * @date     September, 2018
 * @brief    Definition of the DraggableFramelessWidget Class.
 */

#ifndef CONNECTIVITYSETTINGSMANAGER_H
#define CONNECTIVITYSETTINGSMANAGER_H

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

#include <dsp/rt/rt_connectivity.h>

#include <connectivity/metrics/abstractmetric.h>
#include <connectivity/connectivitysettings.h>
#include <connectivity/network/network.h>
#include <connectivity/network/networknode.h>
#include <connectivity/network/networkedge.h>

#include <disp/plots/plot.h>
#include <disp/plots/tfplot.h>
#include <disp/plots/imagesc.h>
#include <dsp/spectrogram.h>
#include <mne/mne_epoch_data_list.h>

//=============================================================================================================
// QT INCLUDES
//=============================================================================================================

#include <QWidget>
#include <QDebug>
#include <QElapsedTimer>

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

//=============================================================================================================
// USED NAMESPACES
//=============================================================================================================

using namespace CONNECTIVITYLIB;
using namespace RTPROCESSINGLIB;
using namespace Eigen;
using namespace DISPLIB;
using namespace UTILSLIB;

//=============================================================================================================
// FORWARD DECLARATIONS
//=============================================================================================================

//=============================================================================================================
/**
 * DECLARE CLASS ConnectivitySettingsManager
 *
 * @brief The ConnectivitySettingsManager class provides a manager to handle connectivty computation for the ex_connectivity example.
 */
class ConnectivitySettingsManager : public QObject
{
    Q_OBJECT

public:
    ConnectivitySettingsManager(QObject* parent = 0)
    : QObject(parent)
    , m_pRtConnectivity(RtConnectivity::SPtr::create())
    {
        QObject::connect(m_pRtConnectivity.data(), &RtConnectivity::newConnectivityResultAvailable,
                         this, &ConnectivitySettingsManager::onNewConnectivityResultAvailable);

        // Default frequency range
        m_fFreqBandLow = 7.0f;
        m_fFreqBandHigh = 13.0f;

        //set coloring
        mColor.insert("COR", Vector4i(90, 26, 100, 1));
        mColor.insert("XCOR", Vector4i(255, 20, 80, 1));
        mColor.insert("PLI", Vector4i(255, 255, 0, 1));
        mColor.insert("COH", Vector4i(2, 89, 100, 1));
        mColor.insert("IMAGCOH", Vector4i(50, 255, 48, 1));
        mColor.insert("PLV", Vector4i(0, 255, 255, 1));
        mColor.insert("WPLI", Vector4i(255, 0, 100, 1));
        mColor.insert("USPLI", Vector4i(255, 89, 200, 1));
        mColor.insert("DSWPLI", Vector4i(25, 10, 255, 1));
    }

    ConnectivitySettings m_settings;
    RtConnectivity::SPtr m_pRtConnectivity;
    QList<Network> m_networkData;
    QMap<QString, Vector4i> mColor;

    float m_fFreqBandLow;
    float m_fFreqBandHigh;

    QVector<int> m_indexList;

    QList<ConnectivitySettings::IntermediateTrialData> m_dataListOriginal;

    DISPLIB::Plot* m_pSignalCoursePlot = nullptr;
    DISPLIB::Plot* m_pSpectrumPlot = nullptr;
    DISPLIB::Plot* m_pEvokedSignalCoursePlot = nullptr;
    DISPLIB::Plot* m_pEpochSignalCoursePlot = nullptr;
    DISPLIB::Plot* m_pEvokedSourceSignalCoursePlot = nullptr;
    DISPLIB::ImageSc* m_pImageConnWeights = nullptr;

    TFplot::SPtr m_pTfPlot;
    Eigen::MatrixXd m_matEvoked;
    Eigen::MatrixXd m_matEvokedSource;
    MNELIB::MNEEpochDataList epochs;

    void onConnectivityMetricChanged(const QString& sMetric)
    {
        if (m_settings.getConnectivityMethods().contains(sMetric)) {
            return;
        }

        m_pRtConnectivity->restart();

        m_settings.setConnectivityMethods(QStringList() << sMetric);

        m_pRtConnectivity->append(m_settings);
    }

    void onNumberTrialsChanged(int iNumberTrials)
    {
        //        QElapsedTimer timer;
        //        qint64 iTime = 0;
        //        timer.start();

        //The maximum number of trials will always be the number of orginal trials stored
        if (iNumberTrials > m_dataListOriginal.size()) {
            iNumberTrials = m_dataListOriginal.size();
        }

        //Pop data from connectivity settings
        int size = m_settings.size();

        if (size > iNumberTrials) {
            m_settings.removeLast(size - iNumberTrials);
        } else {
            while (m_settings.size() < iNumberTrials) {
                //            bool finish = false;
                //            int index = 0;

                //            while(!finish) {
                //                index = rand() % iNumberTrials;

                //                if(!m_indexList.contains(index)) {
                //                    m_indexList.append(index);
                //                    finish = true;
                //                }
                //            }

                m_settings.append(m_dataListOriginal.at(m_settings.size()));
            }
        }

        //qDebug() << "ConnectivitySettingsManager::onNumberTrialsChanged - m_indexList" << m_indexList;

        m_pRtConnectivity->append(m_settings);

        //        iTime = timer.elapsed();
        //        qDebug() << "Coherency::computeCoherencyImag timer - Preparation:" << iTime;
        //        timer.restart();
    }

    void onFreqBandChanged(float fFreqLow, float fFreqHigh)
    {
        if (m_settings.isEmpty()) {
            return;
        }

        // Convert to frequency bins
        m_fFreqBandLow = fFreqLow;
        m_fFreqBandHigh = fFreqHigh;

        onNewConnectivityResultAvailable(m_networkData, m_settings);
    }

    void onNewConnectivityResultAvailable(const QList<Network>& connectivityResults,
                                          const ConnectivitySettings& connectivitySettings)
    {
        m_settings = connectivitySettings;
        m_networkData = connectivityResults;

        for (int i = 0; i < connectivityResults.size(); ++i) {
            m_networkData[i].setFrequencyRange(m_fFreqBandLow, m_fFreqBandHigh);
            m_networkData[i].normalize();

            VisualizationInfo visInfo = m_networkData[i].getVisualizationInfo();
            visInfo.sMethod = "Color";
            visInfo.colNodes = mColor[m_networkData[i].getConnectivityMethod()];
            visInfo.colEdges = mColor[m_networkData[i].getConnectivityMethod()];
            m_networkData[i].setVisualizationInfo(visInfo);

            //            if(!m_networkData.isEmpty()) {
            //                Network network = m_networkData.first();
            //                Eigen::MatrixXd image;

            //                for(int i = 0; i < network.getNodes().size(); i++) {
            //                    for(int j = 0; j < network.getNodes().at(i)->getFullEdges().size(); j++) {
            //                        NetworkEdge::SPtr edge = network.getNodes().at(i)->getFullEdges().at(j);

            //                        if(edge->isActive()) {
            //                            if(image.cols() == 0) {
            //                                image = edge->getMatrixWeight();
            //                            } else {
            //                                image.conservativeResize(image.rows(),image.cols()+1);
            //                                image.col(image.cols()-1) = edge->getMatrixWeight();
            //                            }
            //                        }
            //                    }
            //                }

            //                image.conservativeResize(image.rows()-1,image.cols());

            //                if(!m_pImageConnWeights) {
            //                    m_pImageConnWeights = new DISPLIB::ImageSc(image);
            //                } else {
            //                    m_pImageConnWeights->updateData(image);
            //                }
            //                m_pImageConnWeights->show();
            //            }
        }

        if (!m_networkData.isEmpty()) {
            for (int i = 0; i < m_networkData.size(); ++i) {
                emit newConnectivityResultAvailable("sample",
                                                    "1",
                                                    m_networkData.at(i));
            }
        }
    }

    void plotTimeCourses(int iTrialNumber, int iRowNumber)
    {
        if (iTrialNumber < m_settings.size()) {
            if (iRowNumber < m_settings.at(iTrialNumber).matData.rows()) {
                Eigen::RowVectorXd plotVeca = m_settings.at(iTrialNumber).matData.row(iRowNumber).array() - m_settings.at(iTrialNumber).matData.row(iRowNumber).mean();
                Eigen::Map<Eigen::VectorXd> v1a(plotVeca.data(), plotVeca.size());
                Eigen::VectorXd tempa = v1a.array() - v1a.mean();
                if (!m_pSignalCoursePlot) {
                    m_pSignalCoursePlot = new DISPLIB::Plot(tempa);
                } else {
                    m_pSignalCoursePlot->updateData(tempa);
                }

                m_pSignalCoursePlot->setTitle(QString("Conn used signal for trial %1 and source %2").arg(QString::number(iTrialNumber)).arg(QString::number(iRowNumber)));
                m_pSignalCoursePlot->show();

                //            Eigen::MatrixXd dataSpectrum = Spectrogram::makeSpectrogram(plotVeca, m_settings.getSamplingFrequency()*0.05);

                //            m_pTfPlot = TFplot::SPtr::create(dataSpectrum, m_settings.getSamplingFrequency(), 2,50, ColorMaps::Jet);
                //            m_pTfPlot->show();
            }

            if (iRowNumber < epochs.at(iTrialNumber)->epoch.rows()) {
                Eigen::RowVectorXd plotVecc = epochs.at(iTrialNumber)->epoch.row(iRowNumber);
                Eigen::Map<Eigen::VectorXd> v1c(plotVecc.data(), plotVecc.size());
                Eigen::VectorXd tempc = v1c;
                if (!m_pEpochSignalCoursePlot) {
                    m_pEpochSignalCoursePlot = new DISPLIB::Plot(tempc);
                } else {
                    m_pEpochSignalCoursePlot->updateData(tempc);
                }

                m_pEpochSignalCoursePlot->setTitle(QString("Sensor signal for trial %1 and row %2").arg(QString::number(iTrialNumber)).arg(QString::number(iRowNumber)));
                m_pEpochSignalCoursePlot->show();
            }

            if (iRowNumber < m_settings.at(iTrialNumber).vecTapSpectra.size()) {
                Eigen::RowVectorXd plotVec = m_settings.at(iTrialNumber).vecTapSpectra.at(iRowNumber).cwiseAbs().row(0);
                Eigen::Map<Eigen::VectorXd> v1(plotVec.data(), plotVec.size());
                Eigen::VectorXd temp = v1;
                if (!m_pSpectrumPlot) {
                    m_pSpectrumPlot = new DISPLIB::Plot(temp);
                } else {
                    m_pSpectrumPlot->updateData(temp);
                }

                m_pSpectrumPlot->setTitle(QString("Spectrum of conn used signal for trial %1 and row %2").arg(QString::number(iTrialNumber)).arg(QString::number(iRowNumber)));
                m_pSpectrumPlot->show();
            }
        }

        if (iRowNumber < m_matEvoked.rows()) {
            Eigen::RowVectorXd plotVeca = m_matEvoked.row(iRowNumber);
            Eigen::Map<Eigen::VectorXd> v1a(plotVeca.data(), plotVeca.size());
            Eigen::VectorXd tempa = v1a;
            if (!m_pEvokedSignalCoursePlot) {
                m_pEvokedSignalCoursePlot = new DISPLIB::Plot(tempa);
            } else {
                m_pEvokedSignalCoursePlot->updateData(tempa);
            }

            m_pEvokedSignalCoursePlot->setTitle(QString("Evoked signal for row %1").arg(QString::number(iRowNumber)));
            m_pEvokedSignalCoursePlot->show();
        }

        if (iRowNumber < m_matEvokedSource.rows()) {
            Eigen::RowVectorXd plotVeca = m_matEvokedSource.row(iRowNumber);
            Eigen::Map<Eigen::VectorXd> v1a(plotVeca.data(), plotVeca.size());
            Eigen::VectorXd tempa = v1a;
            if (!m_pEvokedSourceSignalCoursePlot) {
                m_pEvokedSourceSignalCoursePlot = new DISPLIB::Plot(tempa);
            } else {
                m_pEvokedSourceSignalCoursePlot->updateData(tempa);
            }

            m_pEvokedSourceSignalCoursePlot->setTitle(QString("Evoked source signal for row %1").arg(QString::number(iRowNumber)));
            m_pEvokedSourceSignalCoursePlot->show();
        }
    }

signals:
    void newConnectivityResultAvailable(const QString& sSubject,
                                        const QString& sMeasurement,
                                        const Network& tNetworkData);
};

#endif // CONNECTIVITYSETTINGSMANAGER_H

src/examples/ex_connectivity/main.cpp

//=============================================================================================================
/**
 * SPDX-License-Identifier: BSD-3-Clause
 * Copyright (c) 2016-2026 MNE-CPP Authors
 *
 * @file     main.cpp
 * @author   Gabriel Motta <gabrielbenmotta@gmail.com>;
 *           johaenns <j.vorw01@gmail.com>;
 *           Christoph Dinh <christoph.dinh@mne-cpp.org>;
 *           Daniel Strohmeier <daniel.strohmeier@gmail.com>;
 *           Lorenz Esch <lorenz.esch@tu-ilmenau.de>;
 *           Faris Yahya <mfarisyahya@gmail.com>
 * @since    0.1.0
 * @date     July, 2016
 * @brief    Example of using the MNE-CPP Connectivity library
 */

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

#include "connectivitysettingsmanager.h"

#include <disp3D/view/brainview.h>
#include <disp3D/model/braintreemodel.h>
#include <disp3D/model/items/networktreeitem.h>

#include <connectivity/connectivity.h>
#include <connectivity/connectivitysettings.h>
#include <connectivity/network/network.h>

#include <fiff/fiff_raw_data.h>

#include <fs/fs_label.h>
#include <fs/fs_annotationset.h>
#include <fs/fs_surfaceset.h>

#include <inv/inv_source_estimate.h>
#include <mne/mne_epoch_data_list.h>
#include <mne/mne.h>

#include <inv/minimum_norm/inv_minimum_norm.h>

#include <utils/ioutils.h>
#include <utils/generics/mne_logger.h>

#include <disp/viewers/connectivitysettingsview.h>
#include <disp/viewers/minimumnormsettingsview.h>
#include <disp/viewers/tfsettingsview.h>


//=============================================================================================================
// QT INCLUDES
//=============================================================================================================

#include <QApplication>
#include <QMainWindow>
#include <QCommandLineParser>
#include <QDebug>
#include <QFile>
#include <QObject>
#include <QVariant>
#include <QDir>
#include <QMatrix4x4>

//=============================================================================================================
// USED NAMESPACES
//=============================================================================================================

using namespace DISPLIB;
using namespace INVLIB;
using namespace Eigen;
using namespace FIFFLIB;
using namespace CONNECTIVITYLIB;
using namespace Eigen;
using namespace UTILSLIB;
using namespace MNELIB;
using namespace FSLIB;
using namespace DISP3DLIB;

//=============================================================================================================
// MAIN
//=============================================================================================================

//=============================================================================================================
/**
 * The function main marks the entry point of the program.
 * By default, main has the storage class extern.
 *
 * @param[in] argc (argument count) is an integer that indicates how many arguments were entered on the command line when the program was started.
 * @param[in] argv (argument vector) is an array of pointers to arrays of character objects. The array objects are null-terminated strings, representing the arguments that were entered on the command line when the program was started.
 * @return the value that was set to exit() (which is 0 if exit() is called via quit()).
 */
int main(int argc, char* argv[])
{
#ifdef STATICBUILD
// Q_INIT_RESOURCE(mne_disp3d);
#endif
    qInstallMessageHandler(MNELogger::customLogWriter);
    QApplication a(argc, argv);
    QApplication::addLibraryPath(QApplication::applicationDirPath() + "/../lib");

    AbstractMetric::m_bStorageModeIsActive = false;
    AbstractMetric::m_iNumberBinStart = 0;
    AbstractMetric::m_iNumberBinAmount = 50;

    QCommandLineParser parser;
    parser.setApplicationDescription("Connectivity Example");
    parser.addHelpOption();

    QCommandLineOption rawFileOption("fileIn", "The input file <in>.", "in", QCoreApplication::applicationDirPath() + "/../resources/data/MNE-sample-data/MEG/sample/sample_audvis_raw.fif");
    QCommandLineOption eventsFileOption("eve", "Path to the event <file>.", "file", QCoreApplication::applicationDirPath() + "/../resources/data/MNE-sample-data/MEG/sample/sample_audvis_raw-eve.fif");
    QCommandLineOption fwdOption("fwd", "Path to forwad solution <file>.", "file", QCoreApplication::applicationDirPath() + "/../resources/data/MNE-sample-data/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif");
    QCommandLineOption subjectOption("subject", "Selected subject <subject>.", "subject", "sample");
    QCommandLineOption subjectPathOption("subjectPath", "Selected subject path <subjectPath>.", "subjectPath", QCoreApplication::applicationDirPath() + "/../resources/data/MNE-sample-data/subjects");
    QCommandLineOption covFileOption("cov", "Path to the covariance <file>.", "file", QCoreApplication::applicationDirPath() + "/../resources/data/MNE-sample-data/MEG/sample/sample_audvis-cov.fif");
    QCommandLineOption annotOption("annotType", "FsAnnotation <type> (for source level usage only).", "type", "aparc.a2009s");

    QCommandLineOption sourceLocOption("doSourceLoc", "Do source localization (for source level usage only).", "doSourceLoc", "true");
    QCommandLineOption clustOption("doClust", "Do clustering of source space (for source level usage only).", "doClust", "true");
    QCommandLineOption sourceLocMethodOption("sourceLocMethod", "Inverse estimation <method> (for source level usage only), i.e., 'MNE', 'dSPM' or 'sLORETA'.", "method", "dSPM");
    QCommandLineOption connectMethodOption("connectMethod", "Connectivity <method>, i.e., 'COR', 'XCOR.", "method", "IMAGCOH");
    QCommandLineOption snrOption("snr", "The SNR <value> used for computation (for source level usage only).", "value", "3.0");
    QCommandLineOption evokedIndexOption("aveIdx", "The average <index> to choose from the average file.", "index", "1");
    QCommandLineOption chTypeOption("chType", "The channel <type> (for sensor level usage only), i.e. 'eeg' or 'meg'.", "type", "meg");
    QCommandLineOption coilTypeOption("coilType", "The coil <type> (for sensor level usage only), i.e. 'grad' or 'mag'.", "type", "grad");
    QCommandLineOption tMinOption("tmin", "The time minimum value for averaging in seconds relativ to the trigger onset.", "value", "-0.1");
    QCommandLineOption tMaxOption("tmax", "The time maximum value for averaging in seconds relativ to the trigger onset.", "value", "0.4");

    parser.addOption(annotOption);
    parser.addOption(subjectOption);
    parser.addOption(subjectPathOption);
    parser.addOption(fwdOption);
    parser.addOption(sourceLocOption);
    parser.addOption(clustOption);
    parser.addOption(covFileOption);
    parser.addOption(connectMethodOption);
    parser.addOption(sourceLocMethodOption);
    parser.addOption(snrOption);
    parser.addOption(evokedIndexOption);
    parser.addOption(coilTypeOption);
    parser.addOption(chTypeOption);
    parser.addOption(eventsFileOption);
    parser.addOption(rawFileOption);
    parser.addOption(tMinOption);
    parser.addOption(tMaxOption);

    parser.process(a);

    // Init from arguments
    QString sConnectivityMethod = parser.value(connectMethodOption);
    QString sAnnotType = parser.value(annotOption);
    QString sSubj = parser.value(subjectOption);
    QString sSubjDir = parser.value(subjectPathOption);
    QString sFwd = parser.value(fwdOption);
    QString sCov = parser.value(covFileOption);
    QString sSourceLocMethod = parser.value(sourceLocMethodOption);
    QString sCoilType = parser.value(coilTypeOption);
    QString sChType = parser.value(chTypeOption);
    QString sEve = parser.value(eventsFileOption);
    QString sRaw = parser.value(rawFileOption);
    float fTMin = parser.value(tMinOption).toFloat();
    float fTMax = parser.value(tMaxOption).toFloat();
    double dSnr = parser.value(snrOption).toDouble();
    int iEvent = parser.value(evokedIndexOption).toInt();

    bool bDoSourceLoc = false;
    bool bDoClust = false;
    if (parser.value(sourceLocOption) == "false" || parser.value(sourceLocOption) == "0") {
        bDoSourceLoc = false;
    } else if (parser.value(sourceLocOption) == "true" || parser.value(sourceLocOption) == "1") {
        bDoSourceLoc = true;
    }

    if (parser.value(clustOption) == "false" || parser.value(clustOption) == "0") {
        bDoClust = false;
    } else if (parser.value(clustOption) == "true" || parser.value(clustOption) == "1") {
        bDoClust = true;
    }

    //Set parameters
    QList<MatrixXd> matDataList;
    MatrixXi events;
    RowVectorXi picks;

    MNEForwardSolution t_clusteredFwd;
    MNEForwardSolution t_Fwd;

    FsSurfaceSet tSurfSetInflated(sSubj, 2, "inflated", sSubjDir);
    FsAnnotationSet tAnnotSet(sSubj, 2, sAnnotType, sSubjDir);

    QFile coordTransfile(QCoreApplication::applicationDirPath() + "/../resources/data/MNE-sample-data/MEG/sample/all-trans.fif");

    bool keep_comp = false;
    fiff_int_t dest_comp = 0;

    // Create sensor level data
    QFile t_fileRaw(sRaw);
    FiffRawData raw(t_fileRaw);

    int samplesToCutOut = (abs(fTMin) + 0.0) * raw.info.sfreq;
    QSharedPointer<ConnectivitySettingsManager> pConnectivitySettingsManager = QSharedPointer<ConnectivitySettingsManager>::create();

    // Select bad channels
    //raw.info.bads << "MEG2412" << "MEG2413";

    // Setup compensators and projectors so they get applied while reading
    MNE::setup_compensators(raw,
                            dest_comp,
                            keep_comp);

    // Read the events
    MNE::read_events(sEve,
                     sRaw,
                     events);

    // Read the epochs and reject bad epochs. Note, that SSPs are automatically applied to the data if MNE::setup_compensators was called beforehand.
    QStringList exludeChs;
    //exludeChs << "EOG061" << "EOG063";

    QMap<QString, double> mapReject;
    mapReject.insert("eog", 100e-06);
    mapReject.insert("grad", 3000e-13);
    mapReject.insert("mag", 3.5e-12);

    MNEEpochDataList data = MNEEpochDataList::readEpochs(raw,
                                                         events, //events.block(0,0,400,events.cols()),
                                                         fTMin,
                                                         fTMax,
                                                         iEvent,
                                                         mapReject,
                                                         exludeChs);
    data.dropRejected();

    QPair<float, float> pair(fTMin, 0.0f);
    data.applyBaselineCorrection(pair);

    // Average epochs. Do not use SSPs.
    FiffEvoked evoked = data.average(raw.info,
                                     fTMin,
                                     fTMax);

    InvSourceEstimate sourceEstimateEvoked;
    VectorXi vDataIndices;

    QStringList exclude;
    exclude << raw.info.bads << raw.info.ch_names.filter("EOG");

    if (!bDoSourceLoc) {
        // Pick relevant channels
        if (sChType.contains("EEG", Qt::CaseInsensitive)) {
            picks = raw.info.pick_types(false, true, false, QStringList(), exclude);
        } else if (sCoilType.contains("grad", Qt::CaseInsensitive)) {
            // Only pick every second gradiometer which are not marked as bad.
            RowVectorXi picksTmp = raw.info.pick_types(QString("grad"), false, false, QStringList(), exclude);
            picks.resize(0);

            for (int i = 0; i < picksTmp.cols() - 1; i += 2) {
                picks.conservativeResize(picks.cols() + 1);
                picks(picks.cols() - 1) = picksTmp(i);
            }
        } else if (sCoilType.contains("mag", Qt::CaseInsensitive)) {
            picks = raw.info.pick_types(QString("mag"), false, false, QStringList(), exclude);
        }

        // Transform to a more generic data matrix list, pick only channels of interest and remove EOG channel
        MatrixXd matData;
        int iNumberRows = picks.cols(); //picks.cols() 32

        vDataIndices = VectorXi::LinSpaced(iNumberRows, 0, iNumberRows);

        for (int i = 0; i < data.size(); ++i) {
            matData.resize(iNumberRows, data.at(i)->epoch.cols());

            for (qint32 j = 0; j < iNumberRows; ++j) {
                matData.row(j) = data.at(i)->epoch.row(picks(j));
            }
            matDataList << matData;
        }

        // Generate network nodes
        pConnectivitySettingsManager->m_settings.setNodePositions(raw.info, picks);
    } else {
        //Create source level data
        QFile t_fileFwd(sFwd);
        t_Fwd = MNEForwardSolution(t_fileFwd, false, true);

        // Load data
        InvSourceEstimate sourceEstimate;

        double lambda2 = 1.0 / pow(dSnr, 2);
        QString method(sSourceLocMethod);

        // regularize noise covariance
        QFile t_fileCov(sCov);
        FiffCov noise_cov(t_fileCov);
        noise_cov = noise_cov.regularize(raw.info,
                                         0.05,
                                         0.05,
                                         0.1,
                                         true);

        // Cluster forward solution;
        if (bDoClust) {
            t_clusteredFwd = t_Fwd.cluster_forward_solution(tAnnotSet, 40);
        } else {
            t_clusteredFwd = t_Fwd;
        }

        MNEInverseOperator inverse_operator(raw.info,
                                            t_clusteredFwd,
                                            noise_cov,
                                            0.2f,
                                            0.8f);

        // Compute inverse solution
        InvMinimumNorm minimumNorm(inverse_operator, lambda2, method);
        minimumNorm.doInverseSetup(1, true);

        picks = raw.info.pick_types(QString("all"), true, false, QStringList(), exclude);
        data.pick_channels(picks);
        for (int i = 0; i < data.size(); i++) {
            sourceEstimate = minimumNorm.calculateInverse(data.at(i)->epoch,
                                                          evoked.times[0],
                                                          //0.0f,
                                                          1.0 / raw.info.sfreq,
                                                          true);

            if (sourceEstimate.isEmpty()) {
                printf("Source estimate is empty");
            } else {
                matDataList << sourceEstimate.data;
            }
        }

        InvMinimumNorm minimumNormEvoked(inverse_operator, lambda2, method);
        sourceEstimateEvoked = minimumNormEvoked.calculateInverse(evoked, false);
        pConnectivitySettingsManager->m_matEvoked = evoked.data;
        pConnectivitySettingsManager->m_matEvokedSource = sourceEstimateEvoked.data;

        //Generate network nodes and define ROIs
        /*QList<FsLabel> lLabels;
        QList<RowVector4i> qListLabelRGBAs;
        QStringList lWantedLabels;

        lWantedLabels << "G_frontal_inf-Opercular_part-lh"
                        << "G_postcentral-lh"
                        << "G_precentral-lh"
                        << "G_subcentral-lh"
                        << "S_central-lh"
                        << "G_frontal_inf-Opercular_part-rh"
                        << "G_postcentral-rh"
                        << "G_precentral-rh"
                        << "G_subcentral-rh"
                        << "S_central-rh"
                        << "G_occipital_middle-lh"
                        << "G_occipital_middle-rh";

        tAnnotSet.toLabels(tSurfSetInflated, lLabels, qListLabelRGBAs, lWantedLabels);

        //Get active source indices based on picked labels
        vDataIndices = sourceEstimateEvoked.getIndicesByLabel(lLabels, bDoClust);

        //Generate node vertices based on picked labels
        MatrixX3f matNodePositions = t_clusteredFwd.getSourcePositionsByLabel(lLabels, tSurfSetInflated);
        pConnectivitySettingsManager->m_settings.setNodePositions(matNodePositions);*/
        pConnectivitySettingsManager->m_settings.setNodePositions(t_clusteredFwd, tSurfSetInflated);
    }

    pConnectivitySettingsManager->epochs = data;

    //Do connectivity estimation and visualize results
    pConnectivitySettingsManager->m_settings.setConnectivityMethods(QStringList() << sConnectivityMethod);
    pConnectivitySettingsManager->m_settings.setSamplingFrequency(raw.info.sfreq);
    pConnectivitySettingsManager->m_settings.setWindowType("hanning");

    ConnectivitySettings::IntermediateTrialData connectivityData;
    int iNumTrials = matDataList.size();
    for (int i = 0; i < iNumTrials; i++) {
        /*connectivityData.matData.resize(vDataIndices.rows(),matDataList.at(i).cols()-samplesToCutOut);

        for(int j = 0; j < vDataIndices.rows(); j++) {
            // Only calculate connectivity for post stim
            connectivityData.matData.row(j) = matDataList.at(i).row(vDataIndices(j)).segment(samplesToCutOut, matDataList.at(i).cols()-samplesToCutOut);
        }*/
        connectivityData.matData = matDataList.at(i).block(0,
                                                           samplesToCutOut,
                                                           matDataList.at(i).rows(),
                                                           matDataList.at(i).cols() - samplesToCutOut);

        pConnectivitySettingsManager->m_settings.append(connectivityData);
        pConnectivitySettingsManager->m_dataListOriginal.append(connectivityData);
    }

    //Create NetworkView
    // Create BrainView for visualization
    BrainView* pBrainView = new BrainView();
    BrainTreeModel* pModel = new BrainTreeModel();
    pBrainView->setModel(pModel);

    // Create ConnectivitySettingsView as standalone widget
    QSharedPointer<ConnectivitySettingsView> pConnSettingsView = QSharedPointer<ConnectivitySettingsView>::create("EX_CONNECTIVITY");

    QObject::connect(pConnSettingsView.data(), &ConnectivitySettingsView::connectivityMetricChanged,
                     pConnectivitySettingsManager.data(), &ConnectivitySettingsManager::onConnectivityMetricChanged);

    QObject::connect(pConnSettingsView.data(), &ConnectivitySettingsView::numberTrialsChanged,
                     [&pConnectivitySettingsManager, &pConnSettingsView, iNumTrials](int a) {
                         pConnSettingsView->setNumberTrials(qMin(iNumTrials,a));
                         pConnectivitySettingsManager.data()->onNumberTrialsChanged(qMin(iNumTrials,a)); });

    QObject::connect(pConnSettingsView.data(), &ConnectivitySettingsView::freqBandChanged,
                     pConnectivitySettingsManager.data(), &ConnectivitySettingsManager::onFreqBandChanged);

    QObject::connect(pConnectivitySettingsManager.data(), &ConnectivitySettingsManager::newConnectivityResultAvailable,
                     [&]([[maybe_unused]] const QString& a, const QString& b, const Network& c) {
                         pBrainView->loadNetwork(c, b);
                         pBrainView->setNetworkThreshold(0.9);
                     });

    //Read and show sensor helmets
    if (!bDoSourceLoc && sChType.contains("meg", Qt::CaseInsensitive)) {
        //Read and show sensor helmets
        pBrainView->loadSensors(QCoreApplication::applicationDirPath() + "/../resources/general/sensorSurfaces/306m_rt.fif");
    } else {
        // Add source estimate for source-level visualization
        int nVertLh = t_clusteredFwd.src[0].nuse;
        InvSourceEstimate stcLh, stcRh;
        stcLh.data = sourceEstimateEvoked.data.topRows(nVertLh);
        stcLh.vertices = sourceEstimateEvoked.vertices.head(nVertLh);
        stcLh.tmin = sourceEstimateEvoked.tmin;
        stcLh.tstep = sourceEstimateEvoked.tstep;
        stcLh.times = sourceEstimateEvoked.times;
        stcRh.data = sourceEstimateEvoked.data.bottomRows(sourceEstimateEvoked.data.rows() - nVertLh);
        stcRh.vertices = sourceEstimateEvoked.vertices.tail(sourceEstimateEvoked.vertices.size() - nVertLh);
        stcRh.tmin = sourceEstimateEvoked.tmin;
        stcRh.tstep = sourceEstimateEvoked.tstep;
        stcRh.times = sourceEstimateEvoked.times;

        QString tmpDir = QDir::tempPath();
        QString lhStcPath = tmpDir + "/mnecpp_connectivity-lh.stc";
        QString rhStcPath = tmpDir + "/mnecpp_connectivity-rh.stc";
        QFile lhStcFile(lhStcPath);
        stcLh.write(lhStcFile);
        QFile rhStcFile(rhStcPath);
        stcRh.write(rhStcFile);

        // Add hemisphere surfaces
        for (auto it = tSurfSetInflated.data().constBegin(); it != tSurfSetInflated.data().constEnd(); ++it) {
            int hIdx = it.key();
            QString hemi = (it.value().hemi() == 0) ? "lh" : "rh";
            QString surfType = it.value().surf().isEmpty() ? "inflated" : it.value().surf();
            pModel->addSurface("sample", hemi, surfType, it.value());
            if (tAnnotSet.size() > hIdx)
                pModel->addAnnotation("sample", hemi, tAnnotSet[hIdx]);
        }

        // Load source estimate and configure
        QObject::connect(pBrainView, &BrainView::sourceEstimateLoaded, [&](int /*nTimePoints*/) {
            pBrainView->setSourceColormap("Jet");
            pBrainView->setSourceThresholds(0.0f, 0.5f, 10.0f);
        });
        pBrainView->loadSourceEstimate(lhStcPath, rhStcPath);
    }

    pConnSettingsView->setNumberTrials(iNumTrials);
    pConnectivitySettingsManager->onNumberTrialsChanged(iNumTrials);
    pConnectivitySettingsManager->onFreqBandChanged(pConnSettingsView->getLowerFreq(), pConnSettingsView->getUpperFreq());

    pConnSettingsView->show();
    pBrainView->show();

    return a.exec();
}