ex_connectivity_comparison

Example to compare connectivity methods.

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

Options

Option	Description
--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'.
--snr	The SNR <value> used for computation (for source level usage only).
--aveIdx	The average <index> to choose from the average file.
--coilType	The coil <type> (for sensor level usage only), i.e. 'grad' or 'mag'.
--chType	The channel <type> (for sensor level usage only), i.e. 'eeg' or 'meg'.
--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.
--eve	Path to the event <file>.
--fileIn	Path to the raw <file>.
--subj	Selected <subject> (for source level usage only).
--subjDir	Selected <subjectPath> (for source level usage only).
--fwd	Path to forwad solution <file> (for source level usage only).
--cov	Path to the covariance <file> (for source level usage only).
--ntrials	The number of averaged trials.
--freqmin	The lower frequency bound for connectivity computation.
--freqmax	The upper frequency bound for connectivity computation.

Build and run

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

The default inputs come from the MNE sample dataset.

Used in the API reference

• Connectivity
• ConnectivitySettings
• Network
• VisualizationInfo
• ConnectivitySettingsView
• MNEForwardSolution

Source

src/examples/ex_connectivity_comparison/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>;
 *           Juan GPC <jgarciaprieto@mgh.harvard.edu>;
 *           Lorenz Esch <lorenz.esch@tu-ilmenau.de>
 * @since    0.1.0
 * @date     July, 2016
 * @brief     Example to compare connectivity methods
 */

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

#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 <connectivity/metrics/abstractmetric.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>

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

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

//=============================================================================================================
// 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);

    AbstractMetric::m_bStorageModeIsActive = false;
    //    AbstractMetric::m_iNumberBinStart = 8;
    //    AbstractMetric::m_iNumberBinAmount = 4;

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

    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 snrOption("snr", "The SNR <value> used for computation (for source level usage only).", "value", "1.0");
    QCommandLineOption evokedIndexOption("aveIdx", "The average <index> to choose from the average file.", "index", "3");
    QCommandLineOption coilTypeOption("coilType", "The coil <type> (for sensor level usage only), i.e. 'grad' or 'mag'.", "type", "grad");
    QCommandLineOption chTypeOption("chType", "The channel <type> (for sensor level usage only), i.e. 'eeg' or 'meg'.", "type", "meg");
    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");
    QCommandLineOption eventsFileOption("eve", "Path to the event <file>.", "file", QCoreApplication::applicationDirPath() + "/../resources/data/MNE-sample-data/MEG/sample/sample_audvis_raw-eve.fif");
    QCommandLineOption rawFileOption("fileIn", "Path to the raw <file>.", "file", QCoreApplication::applicationDirPath() + "/../resources/data/MNE-sample-data/MEG/sample/sample_audvis_raw.fif");
    QCommandLineOption subjectOption("subj", "Selected <subject> (for source level usage only).", "subject", "sample");
    QCommandLineOption subjectPathOption("subjDir", "Selected <subjectPath> (for source level usage only).", "subjectPath", QCoreApplication::applicationDirPath() + "/../resources/data/MNE-sample-data/subjects");
    QCommandLineOption fwdOption("fwd", "Path to forwad solution <file> (for source level usage only).", "file", QCoreApplication::applicationDirPath() + "/../resources/data/MNE-sample-data/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif");
    QCommandLineOption covFileOption("cov", "Path to the covariance <file> (for source level usage only).", "file", QCoreApplication::applicationDirPath() + "/../resources/data/MNE-sample-data/MEG/sample/sample_audvis-cov.fif");
    QCommandLineOption trialOption("ntrials", "The number of averaged trials.", "value", "200");
    QCommandLineOption freqMinOption("freqmin", "The lower frequency bound for connectivity computation.", "value", "7.0");
    QCommandLineOption freqMaxOption("freqmax", "The upper frequency bound for connectivity computation.", "value", "13.0");
    parser.addOption(annotOption);
    parser.addOption(subjectOption);
    parser.addOption(subjectPathOption);
    parser.addOption(fwdOption);
    parser.addOption(sourceLocOption);
    parser.addOption(clustOption);
    parser.addOption(covFileOption);
    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.addOption(trialOption);
    parser.addOption(freqMinOption);
    parser.addOption(freqMaxOption);

    parser.process(a);

    // Init from arguments
    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();
    int iNumTrials = parser.value(trialOption).toInt();
    double dFreqMin = parser.value(freqMinOption).toDouble();
    double dFreqMax = parser.value(freqMaxOption).toDouble();

    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 t_fileCov(sCov);
    FiffCov noise_cov(t_fileCov);

    bool keep_comp = false;
    fiff_int_t dest_comp = 0;

    ConnectivitySettings conSettings;

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

    // 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.
    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,
                                                         fTMin,
                                                         fTMax,
                                                         iEvent,
                                                         mapReject);
    data.dropRejected();
    QPair<float, float> pair(fTMin, 0.0f);
    data.applyBaselineCorrection(pair);

    FiffEvoked evoked = data.average(raw.info,
                                     fTMin,
                                     fTMax);
    InvSourceEstimate sourceEstimateEvoked;

    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);
            picks.resize(0);

            for (int i = 0; i < picksTmp.cols() - 1; i += 2) {
                if (!raw.info.bads.contains(raw.info.ch_names.at(picksTmp(i)))) {
                    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);
        }

        // Generate network nodes
        conSettings.setNodePositions(raw.info, picks);

        // 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

        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;
        }
    } 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
        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,
                                                          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);

        // Generate network nodes
        conSettings.setNodePositions(t_clusteredFwd, tSurfSetInflated);
    }

    // Compute the connectivity estimates for the methods to be compared
    conSettings.setConnectivityMethods(QStringList() << "COH" << "COR" << "XCOR" << "PLI" << "IMAGCOH" << "PLV" << "WPLI" << "USPLI" << "DSWPLI");

    iNumTrials = qMin(iNumTrials, matDataList.size());
    for (int i = 0; i < iNumTrials; i++) {
        // Only calculate connectivity for post stim
        int samplesToCutOut = abs(fTMin * raw.info.sfreq);
        conSettings.append(matDataList.at(i).block(0,
                                                   samplesToCutOut,
                                                   matDataList.at(i).rows(),
                                                   matDataList.at(i).cols() - samplesToCutOut));
    }

    conSettings.setSamplingFrequency(raw.info.sfreq);
    conSettings.setWindowType("hanning");

    QList<Network> lNetworks = Connectivity::calculate(conSettings);

    QMap<QString, Vector4i> mColor;
    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));

    for (int j = 0; j < lNetworks.size(); ++j) {
        lNetworks[j].setFrequencyRange(dFreqMin, dFreqMax);
        lNetworks[j].normalize();
        VisualizationInfo visInfo = lNetworks.at(j).getVisualizationInfo();
        visInfo.sMethod = "Color";
        visInfo.colNodes = mColor[lNetworks[j].getConnectivityMethod()];
        visInfo.colEdges = mColor[lNetworks[j].getConnectivityMethod()];
        lNetworks[j].setVisualizationInfo(visInfo);
    }

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

    // Load all networks
    for (int j = 0; j < lNetworks.size(); ++j) {
        pBrainView->loadNetwork(lNetworks[j], lNetworks[j].getConnectivityMethod());
    }
    pBrainView->setNetworkThreshold(0.9);

    // Create ConnectivitySettingsView as standalone widget
    QSharedPointer<ConnectivitySettingsView> pConnSettingsView = QSharedPointer<ConnectivitySettingsView>::create("EX_CONNECTIVITY_COMP");
    pConnSettingsView->setNumberTrials(iNumTrials);
    pConnSettingsView->setFrequencyBand(dFreqMin, dFreqMax);

    //Read and show sensor helmets
    if (!bDoSourceLoc && sChType.contains("meg", Qt::CaseInsensitive)) {
        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_comparison-lh.stc";
        QString rhStcPath = tmpDir + "/mnecpp_connectivity_comparison-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);
    }

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

    QObject::connect(pConnSettingsView.data(), &ConnectivitySettingsView::freqBandChanged,
                     [&pConnSettingsView, dFreqMin, dFreqMax]() { pConnSettingsView->setFrequencyBand(dFreqMin, dFreqMax); });

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

    return a.exec();
}