rt_averaging.cpp

Go to the documentation of this file.

//=============================================================================================================
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "rt_averaging.h"
 
#include <mne/mne_epoch_data_list.h>
 
#include <utils/ioutils.h>
#include "rt_detect_trigger.h"
#include <math/numerics.h>
 
//=============================================================================================================
// QT INCLUDES
//=============================================================================================================
 
#include <QDebug>
 
//=============================================================================================================
// EIGEN INCLUDES
//=============================================================================================================
 
//=============================================================================================================
// USED NAMESPACES
//=============================================================================================================
 
using namespace RTPROCESSINGLIB;
using namespace FIFFLIB;
using namespace UTILSLIB;
using namespace Eigen;
using namespace MNELIB;
 
//=============================================================================================================
// DEFINE MEMBER METHODS RtAveragingWorker
//=============================================================================================================
 
RtAveragingWorker::RtAveragingWorker(quint32 numAverages,
                                     quint32 iPreStimSamples,
                                     quint32 iPostStimSamples,
                                     qint32 iBaselineFromMSecs,
                                     qint32 iBaselineToMSecs,
                                     quint32 iTriggerIndex,
                                     FiffInfo::SPtr pFiffInfo)
: QObject()
, m_iNumAverages(numAverages)
, m_iPreStimSamples(iPreStimSamples)
, m_iPostStimSamples(iPostStimSamples)
, m_iTriggerChIndex(-1)
, m_iNewTriggerIndex(iTriggerIndex)
, m_fTriggerThreshold(0.5f)
, m_bActivateThreshold(false)
, m_bDoBaselineCorrection(false)
, m_pairBaselineSec(qMakePair(float(iBaselineFromMSecs) / 1000.0f, float(iBaselineToMSecs) / 1000.0f))
, m_pFiffInfo(pFiffInfo)
{
    m_mapThresholds["eog"] = 300e-6;
 
    m_stimEvokedSet.info = *m_pFiffInfo.data();
 
    m_iNewPreStimSamples = m_iPreStimSamples;
    m_iNewPostStimSamples = m_iPostStimSamples;
 
    if (m_iNumAverages <= 0) {
        qDebug() << "RtAveragingWorker::RtAveragingWorker - Number of averages <= 0. Setting to 1 as default.";
        m_iNumAverages = 1;
    }
}
 
//=============================================================================================================
 
void RtAveragingWorker::doWork(const MatrixXd& rawSegment)
{
    if (this->thread()->isInterruptionRequested()) {
        return;
    }
 
    if (controlValuesChanged()) {
        reset();
    }
 
    doAveraging(rawSegment);
}
 
//=============================================================================================================
 
void RtAveragingWorker::setAverageNumber(qint32 numAve)
{
    if (numAve <= 0) {
        qDebug() << "[RtAveragingWorker::setAverageNumber] Number of averages <= 0 are not allowed. Returning.";
        return;
    }
 
    if (numAve < m_iNumAverages) {
        int iDiff = m_iNumAverages - numAve;
 
        //Do averaging for each trigger type
        QMutableMapIterator<double, QList<Eigen::MatrixXd>> idx(m_mapStimAve);
 
        while (idx.hasNext()) {
            idx.next();
 
            if (idx.value().size() > iDiff) {
                //Pop data from buffer
                for (int i = 0; i < iDiff; ++i) {
                    idx.value().pop_front();
                }
            }
        }
    }
 
    m_iNumAverages = numAve;
}
 
//=============================================================================================================
 
void RtAveragingWorker::setPreStim(qint32 samples, qint32 secs)
{
    Q_UNUSED(secs);
 
    m_iNewPreStimSamples = samples;
}
 
//=============================================================================================================
 
void RtAveragingWorker::setPostStim(qint32 samples, qint32 secs)
{
    Q_UNUSED(secs);
 
    m_iNewPostStimSamples = samples;
}
 
//=============================================================================================================
 
void RtAveragingWorker::setTriggerChIndx(qint32 idx)
{
    m_iNewTriggerIndex = idx;
}
 
//=============================================================================================================
 
void RtAveragingWorker::setArtifactReduction(const QMap<QString, double>& mapThresholds)
{
    if (mapThresholds["Active"] == 0.0) {
        m_bActivateThreshold = false;
    } else {
        m_bActivateThreshold = true;
    }
 
    m_mapThresholds = mapThresholds;
}
 
//=============================================================================================================
 
void RtAveragingWorker::setBaselineActive(bool activate)
{
    m_bDoBaselineCorrection = activate;
 
    if (!m_bDoBaselineCorrection) {
        for (int i = 0; i < m_stimEvokedSet.evoked.size(); ++i) {
            m_stimEvokedSet.evoked[i].baseline = qMakePair(-1.0f, -1.0f);
        }
    } else {
        for (int i = 0; i < m_stimEvokedSet.evoked.size(); ++i) {
            m_stimEvokedSet.evoked[i].baseline = m_pairBaselineSec;
        }
    }
}
 
//=============================================================================================================
 
void RtAveragingWorker::setBaselineFrom(int fromSamp,
                                        int fromMSec)
{
    m_pairBaselineSec.first = float(fromMSec) / 1000.0f;
    m_pairBaselineSamp.first = float(fromSamp);
 
    for (int i = 0; i < m_stimEvokedSet.evoked.size(); ++i) {
        m_stimEvokedSet.evoked[i].baseline.first = float(fromMSec) / 1000.0f;
    }
}
 
//=============================================================================================================
 
void RtAveragingWorker::setBaselineTo(int toSamp,
                                      int toMSec)
{
    m_pairBaselineSec.second = float(toMSec) / 1000.0f;
    m_pairBaselineSamp.second = float(toSamp);
 
    for (int i = 0; i < m_stimEvokedSet.evoked.size(); ++i) {
        m_stimEvokedSet.evoked[i].baseline.second = float(toMSec) / 1000.0f;
    }
}
 
//=============================================================================================================
 
void RtAveragingWorker::doAveraging(const MatrixXd& rawSegment)
{
    //Detect trigger
    QList<QPair<int, double>> lDetectedTriggers = RTPROCESSINGLIB::detectTriggerFlanksMax(rawSegment, m_iTriggerChIndex, 0, m_fTriggerThreshold, true);
 
    //TODO: This does not permit the same trigger type twice in one data block
    for (int i = 0; i < lDetectedTriggers.size(); ++i) {
        if (!m_mapFillingBackBuffer.contains(lDetectedTriggers.at(i).second)) {
            double dTriggerType = lDetectedTriggers.at(i).second;
 
            //qDebug()<<"Adding dTriggerType"<<dTriggerType;
 
            m_mapFillingBackBuffer[dTriggerType] = false;
            m_mapMatDataPostIdx[dTriggerType] = 0;
            m_mapDataPost[dTriggerType].resize(m_pFiffInfo->chs.size(), m_iPostStimSamples);
            m_mapDataPost[dTriggerType].setZero();
        }
    }
 
    // Pre-stimulus samples of a trigger may reach back into the previous blocks
    MatrixXd matHistoryAndBlock(rawSegment.rows(), m_matPreHistory.cols() + rawSegment.cols());
    matHistoryAndBlock << m_matPreHistory, rawSegment;
 
    //Do averaging for each trigger type
    QMutableMapIterator<double, bool> idx(m_mapFillingBackBuffer);
    QStringList lResponsibleTriggerTypes;
 
    while (idx.hasNext()) {
        idx.next();
 
        double dTriggerType = idx.key();
 
        //Fill back buffer and decide when to do the data packing of the different buffers
        if (m_mapFillingBackBuffer[dTriggerType]) {
            if (m_mapMatDataPostIdx[dTriggerType] != m_iPostStimSamples) {
                fillBackBuffer(rawSegment, dTriggerType);
            }
 
            if (m_mapMatDataPostIdx[dTriggerType] == m_iPostStimSamples) {
                m_mapMatDataPostIdx[dTriggerType] = 0;
                m_mapFillingBackBuffer[dTriggerType] = false;
                emitEvoked(dTriggerType, lResponsibleTriggerTypes);
            }
        } else {
            for (int i = 0; i < lDetectedTriggers.size(); ++i) {
                if (dTriggerType == lDetectedTriggers.at(i).second) {
                    int iTriggerPos = lDetectedTriggers.at(i).first;
 
                    m_mapDataPre[dTriggerType] = matHistoryAndBlock.middleCols(iTriggerPos, m_iPreStimSamples);
 
                    //Do back buffer stuff
                    if (rawSegment.cols() - iTriggerPos >= m_mapDataPost[dTriggerType].cols()) {
                        m_mapDataPost[dTriggerType] = rawSegment.block(0,
                                                                       iTriggerPos,
                                                                       m_mapDataPost[dTriggerType].rows(),
                                                                       m_mapDataPost[dTriggerType].cols());
                        emitEvoked(dTriggerType, lResponsibleTriggerTypes);
                        m_mapMatDataPostIdx[dTriggerType] = 0;
                        m_mapFillingBackBuffer[dTriggerType] = false;
                    } else {
                        m_mapDataPost[dTriggerType].block(0,
                                                          0,
                                                          m_mapDataPost[dTriggerType].rows(),
                                                          rawSegment.cols() - iTriggerPos) = rawSegment.block(0,
                                                                                                              iTriggerPos,
                                                                                                              rawSegment.rows(),
                                                                                                              rawSegment.cols() - iTriggerPos);
                        m_mapMatDataPostIdx[dTriggerType] = rawSegment.cols() - iTriggerPos;
                        m_mapFillingBackBuffer[dTriggerType] = true;
                    }
                }
            }
        }
    }
 
    m_matPreHistory = matHistoryAndBlock.rightCols(m_iPreStimSamples);
}
 
//=============================================================================================================
 
void RtAveragingWorker::emitEvoked(double dTriggerType, QStringList& lResponsibleTriggerTypes)
{
    //Merge the different buffers
    mergeData(dTriggerType);
 
    //Calculate the final average/evoked data
    generateEvoked(dTriggerType);
 
    //List of all trigger types which lead to the recent emit of a new evoked set. */
    if (!lResponsibleTriggerTypes.contains(QString::number(dTriggerType))) {
        lResponsibleTriggerTypes << QString::number(dTriggerType);
    }
 
    if (m_stimEvokedSet.evoked.size() > 0) {
        emit resultReady(m_stimEvokedSet, lResponsibleTriggerTypes);
    }
 
    //    qDebug()<<"RtAveragingWorker::emitEvoked() - dTriggerType:" << dTriggerType;
    //    qDebug()<<"RtAveragingWorker::emitEvoked() - m_mapStimAve[dTriggerType].size():" << m_mapStimAve[dTriggerType].size();
}
 
//=============================================================================================================
 
void RtAveragingWorker::fillBackBuffer(const MatrixXd& data, double dTriggerType)
{
    int iResidualCols = data.cols();
    if (m_mapMatDataPostIdx[dTriggerType] + data.cols() > m_iPostStimSamples) {
        iResidualCols = m_iPostStimSamples - m_mapMatDataPostIdx[dTriggerType];
        m_mapDataPost[dTriggerType].block(0, m_mapMatDataPostIdx[dTriggerType], m_mapDataPost[dTriggerType].rows(), iResidualCols) = data.block(0, 0, data.rows(), iResidualCols);
    } else {
        m_mapDataPost[dTriggerType].block(0, m_mapMatDataPostIdx[dTriggerType], m_mapDataPost[dTriggerType].rows(), iResidualCols) = data;
    }
 
    m_mapMatDataPostIdx[dTriggerType] += iResidualCols;
}
 
//=============================================================================================================
 
void RtAveragingWorker::mergeData(double dTriggerType)
{
    if (m_mapDataPre[dTriggerType].rows() != m_mapDataPost[dTriggerType].rows()) {
        qDebug() << "[RtAveragingWorker::mergeData] Rows of m_mapDataPre (" << m_mapDataPre[dTriggerType].rows() << ") and m_mapDataPost (" << m_mapDataPost[dTriggerType].rows() << ") are not the same. Returning.";
        return;
    }
 
    MatrixXd mergedData(m_mapDataPre[dTriggerType].rows(), m_mapDataPre[dTriggerType].cols() + m_mapDataPost[dTriggerType].cols());
 
    mergedData << m_mapDataPre[dTriggerType], m_mapDataPost[dTriggerType];
 
    //Perform artifact threshold
    bool bArtifactDetected = false;
 
    if (m_bActivateThreshold && m_pFiffInfo) {
        qDebug() << "[RtAveragingWorker::mergeData] Doing artifact reduction for" << m_mapThresholds;
 
        bArtifactDetected = MNEEpochDataList::checkForArtifact(mergedData,
                                                               *m_pFiffInfo,
                                                               m_mapThresholds);
    }
 
    if (!bArtifactDetected) {
        //Add cut data to average buffer
        m_mapStimAve[dTriggerType].append(mergedData);
 
        //Pop data from buffer
        int iDiff = m_mapStimAve[dTriggerType].size() - m_iNumAverages;
        if (iDiff > 0) {
            for (int i = 0; i < iDiff; ++i) {
                m_mapStimAve[dTriggerType].pop_front();
            }
        }
    }
}
 
//=============================================================================================================
 
void RtAveragingWorker::generateEvoked(double dTriggerType)
{
    if (m_mapStimAve[dTriggerType].isEmpty()) {
        qDebug() << "[RtAveragingWorker::generateEvoked] m_mapStimAve is empty for type" << dTriggerType << "Returning.";
        return;
    }
 
    //Init evoked
    m_stimEvokedSet.info = *m_pFiffInfo.data();
    FiffEvoked evoked;
    evoked.setInfo(*m_pFiffInfo.data());
    int iEvokedIdx = -1;
 
    for (int i = 0; i < m_stimEvokedSet.evoked.size(); ++i) {
        if (m_stimEvokedSet.evoked.at(i).comment == QString::number(dTriggerType)) {
            evoked = m_stimEvokedSet.evoked.at(i);
            iEvokedIdx = i;
            break;
        }
    }
 
    //If the evoked is not yet present add it here
    if (iEvokedIdx == -1) {
        evoked.baseline = m_pairBaselineSec;
        evoked.times = RowVectorXf::LinSpaced(m_iPreStimSamples + m_iPostStimSamples,
                                              -m_iPreStimSamples,
                                              m_iPostStimSamples - 1) /
            m_pFiffInfo->sfreq;
        evoked.first = -m_iPreStimSamples;
        evoked.last = m_iPostStimSamples - 1;
        evoked.comment = QString::number(dTriggerType);
    }
 
    // Generate final evoked
    MatrixXd finalAverage = MatrixXd::Zero(m_mapStimAve[dTriggerType].first().rows(), m_iPreStimSamples + m_iPostStimSamples);
 
    for (int i = 0; i < m_mapStimAve[dTriggerType].size(); ++i) {
        finalAverage += m_mapStimAve[dTriggerType].at(i);
    }
 
    if (!m_mapStimAve[dTriggerType].isEmpty()) {
        finalAverage = finalAverage / m_mapStimAve[dTriggerType].size();
    }
 
    if (m_bDoBaselineCorrection) {
        finalAverage = Numerics::rescale(finalAverage, evoked.times, m_pairBaselineSec, QString("mean"));
    }
 
    evoked.data = finalAverage;
 
    evoked.nave = m_mapStimAve[dTriggerType].size();
 
    //Add new data to evoked data set
    if (iEvokedIdx != -1) {
        //Evoked data is already present
        m_stimEvokedSet.evoked[iEvokedIdx] = evoked;
    } else {
        //Evoked data is not present yet
        m_stimEvokedSet.evoked.append(evoked);
    }
}
 
//=============================================================================================================
 
void RtAveragingWorker::reset()
{
    //Reset
    m_iPreStimSamples = m_iNewPreStimSamples;
    m_iPostStimSamples = m_iNewPostStimSamples;
    m_iTriggerChIndex = m_iNewTriggerIndex;
 
    //Clear all evoked data information
    m_stimEvokedSet.evoked.clear();
 
    //Clear all maps
    m_mapStimAve.clear();
    m_mapDataPre.clear();
    m_matPreHistory = MatrixXd::Zero(m_pFiffInfo->chs.size(), m_iPreStimSamples);
    m_mapDataPost.clear();
    m_mapMatDataPostIdx.clear();
    m_mapFillingBackBuffer.clear();
}
 
//=============================================================================================================
// DEFINE MEMBER METHODS RtAveraging
//=============================================================================================================
 
RtAveraging::RtAveraging(quint32 numAverages,
                         quint32 iPreStimSamples,
                         quint32 iPostStimSamples,
                         qint32 iBaselineFromMSecs,
                         qint32 iBaselineToMSecs,
                         quint32 iTriggerIndex,
                         FiffInfo::SPtr pFiffInfo,
                         QObject* parent)
: QObject(parent)
{
    qRegisterMetaType<Eigen::MatrixXd>("Eigen::MatrixXd");
 
    RtAveragingWorker* worker = new RtAveragingWorker(numAverages,
                                                      iPreStimSamples,
                                                      iPostStimSamples,
                                                      iBaselineFromMSecs,
                                                      iBaselineToMSecs,
                                                      iTriggerIndex,
                                                      pFiffInfo);
    worker->moveToThread(&m_workerThread);
 
    connect(&m_workerThread, &QThread::finished,
            worker, &QObject::deleteLater);
 
    connect(this, &RtAveraging::operate,
            worker, &RtAveragingWorker::doWork);
 
    connect(worker, &RtAveragingWorker::resultReady,
            this, &RtAveraging::handleResults, Qt::DirectConnection);
 
    connect(this, &RtAveraging::averageNumberChanged,
            worker, &RtAveragingWorker::setAverageNumber);
    connect(this, &RtAveraging::averagePreStimChanged,
            worker, &RtAveragingWorker::setPreStim);
    connect(this, &RtAveraging::averagePostStimChanged,
            worker, &RtAveragingWorker::setPostStim);
    connect(this, &RtAveraging::averageTriggerChIdxChanged,
            worker, &RtAveragingWorker::setTriggerChIndx);
    connect(this, &RtAveraging::averageArtifactReductionChanged,
            worker, &RtAveragingWorker::setArtifactReduction);
    connect(this, &RtAveraging::averageBaselineActiveChanged,
            worker, &RtAveragingWorker::setBaselineActive);
    connect(this, &RtAveraging::averageBaselineFromChanged,
            worker, &RtAveragingWorker::setBaselineFrom);
    connect(this, &RtAveraging::averageBaselineToChanged,
            worker, &RtAveragingWorker::setBaselineTo);
    connect(this, &RtAveraging::averageResetRequested,
            worker, &RtAveragingWorker::reset);
 
    m_workerThread.start();
}
 
//=============================================================================================================
 
RtAveraging::~RtAveraging()
{
    stop();
}
 
//=============================================================================================================
 
void RtAveraging::append(const MatrixXd& data)
{
    emit operate(data);
}
 
//=============================================================================================================
 
void RtAveraging::handleResults(const FiffEvokedSet& evokedStimSet,
                                const QStringList& lResponsibleTriggerTypes)
{
    emit evokedStim(evokedStimSet,
                    lResponsibleTriggerTypes);
}
 
//=============================================================================================================
 
void RtAveraging::restart(quint32 numAverages,
                          quint32 iPreStimSamples,
                          quint32 iPostStimSamples,
                          qint32 iBaselineFromMSecs,
                          qint32 iBaselineToMSecs,
                          quint32 iTriggerIndex,
                          FiffInfo::SPtr pFiffInfo)
{
    stop();
 
    RtAveragingWorker* worker = new RtAveragingWorker(numAverages,
                                                      iPreStimSamples,
                                                      iPostStimSamples,
                                                      iBaselineFromMSecs,
                                                      iBaselineToMSecs,
                                                      iTriggerIndex,
                                                      pFiffInfo);
    worker->moveToThread(&m_workerThread);
 
    connect(&m_workerThread, &QThread::finished,
            worker, &QObject::deleteLater);
 
    connect(this, &RtAveraging::operate,
            worker, &RtAveragingWorker::doWork);
 
    connect(worker, &RtAveragingWorker::resultReady,
            this, &RtAveraging::handleResults, Qt::DirectConnection);
 
    connect(this, &RtAveraging::averageNumberChanged,
            worker, &RtAveragingWorker::setAverageNumber);
    connect(this, &RtAveraging::averagePreStimChanged,
            worker, &RtAveragingWorker::setPreStim);
    connect(this, &RtAveraging::averagePostStimChanged,
            worker, &RtAveragingWorker::setPostStim);
    connect(this, &RtAveraging::averageTriggerChIdxChanged,
            worker, &RtAveragingWorker::setTriggerChIndx);
    connect(this, &RtAveraging::averageArtifactReductionChanged,
            worker, &RtAveragingWorker::setArtifactReduction);
    connect(this, &RtAveraging::averageBaselineActiveChanged,
            worker, &RtAveragingWorker::setBaselineActive);
    connect(this, &RtAveraging::averageBaselineFromChanged,
            worker, &RtAveragingWorker::setBaselineFrom);
    connect(this, &RtAveraging::averageBaselineToChanged,
            worker, &RtAveragingWorker::setBaselineTo);
    connect(this, &RtAveraging::averageResetRequested,
            worker, &RtAveragingWorker::reset);
 
    m_workerThread.start();
}
 
//=============================================================================================================
 
void RtAveraging::stop()
{
    m_workerThread.requestInterruption();
    m_workerThread.quit();
    m_workerThread.wait();
}
 
//=============================================================================================================
 
void RtAveraging::setAverageNumber(qint32 numAve)
{
    emit averageNumberChanged(numAve);
}
 
//=============================================================================================================
 
void RtAveraging::setPreStim(qint32 samples,
                             qint32 secs)
{
    Q_UNUSED(secs);
 
    emit averagePreStimChanged(samples,
                               secs);
}
 
//=============================================================================================================
 
void RtAveraging::setPostStim(qint32 samples,
                              qint32 secs)
{
    Q_UNUSED(secs);
 
    emit averagePostStimChanged(samples,
                                secs);
}
 
//=============================================================================================================
 
void RtAveraging::setTriggerChIndx(qint32 idx)
{
    emit averageTriggerChIdxChanged(idx);
}
 
//=============================================================================================================
 
void RtAveraging::setArtifactReduction(const QMap<QString, double>& mapThresholds)
{
    emit averageArtifactReductionChanged(mapThresholds);
}
 
//=============================================================================================================
 
void RtAveraging::setBaselineActive(bool activate)
{
    emit averageBaselineActiveChanged(activate);
}
 
//=============================================================================================================
 
void RtAveraging::setBaselineFrom(int fromSamp,
                                  int fromMSec)
{
    emit averageBaselineFromChanged(fromSamp,
                                    fromMSec);
}
 
//=============================================================================================================
 
void RtAveraging::setBaselineTo(int toSamp,
                                int toMSec)
{
    emit averageBaselineToChanged(toSamp,
                                  toMSec);
}
 
//=============================================================================================================
 
void RtAveraging::reset()
{
    emit averageResetRequested();
}