Connectivity
Namespace: CONNECTIVITYLIB · Library: Connectivity Library
Python equivalent
mne_connectivity.spectral_connectivity_epochs in MNE-Python.
#include <connectivity/connectivity.h>
class CONNECTIVITYLIB::Connectivity
Dispatcher that runs one or more functional-connectivity estimators over a batch of pre-processed trials and returns one Network per estimator.
The selected methods are read from ConnectivitySettings::getConnectivityMethods and matched by string against the supported metric set ("COR", "XCOR", "COH", "IMAGCOH", "PLI", "USPLI", "WPLI", "DSWPLI", "PLV", "GC", "DTF", "PDC"). All spectral metrics share the same DPSS tapered FFT and CSD accumulator inside ConnectivitySettings, so running e.g. PLI + wPLI + dwPLI in one call costs only one FFT pass.
Runs the selected functional-connectivity metrics over a ConnectivitySettings batch.
Public Methods
Connectivity()
Constructs a Connectivity object.
Static Methods
calculate(connectivitySettings)
Computes the network based on the current settings.
Returns:
- QList< Network > — Returns the list with calculated networks for each provided method.
Example
Source: src/examples/ex_connectivity/main.cpp
#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;
//=============================================================================================================
// 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,
[&](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();
}
Authors of this file
- Lorenz Esch <lorenz.esch@tu-ilmenau.de>
- Gabriel Motta <gabrielbenmotta@gmail.com>
- Christoph Dinh <christoph.dinh@mne-cpp.org>