channeldatamodel.cpp

Go to the documentation of this file.

//=============================================================================================================
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "channeldatamodel.h"
#include <fiff/fiff_info.h>
#include <fiff/fiff_constants.h>
 
//=============================================================================================================
// QT INCLUDES
//=============================================================================================================
 
#include <QtMath>
#include <QWriteLocker>
#include <QReadLocker>
 
//=============================================================================================================
// USED NAMESPACES
//=============================================================================================================
 
using namespace DISPLIB;
using namespace FIFFLIB;
using namespace Eigen;
 
//=============================================================================================================
// Default amplitude scales (physical units) per channel kind.
// Keys match FIFF channel kind constants.
//=============================================================================================================
 
namespace
{
constexpr float kScaleMEGGrad = 400e-13f; // T/m
constexpr float kScaleMEGMag = 1.2e-12f;  // T
constexpr float kScaleEEG = 30e-6f;       // V
constexpr float kScaleEOG = 150e-6f;      // V
constexpr float kScaleEMG = 1e-3f;        // V
constexpr float kScaleECG = 1e-3f;        // V
constexpr float kScaleSTIM = 5.0f;        // AU
constexpr float kScaleMISC = 1.0f;        // AU
constexpr float kScaleFallback = 1.0f;
 
// Internal pseudo-kind keys for separate MEG grad/mag scales
constexpr qint32 kMEGGradKind = -1;
constexpr qint32 kMEGMagKind = -2;
}
 
//=============================================================================================================
// DEFINE MEMBER METHODS
//=============================================================================================================
 
ChannelDataModel::ChannelDataModel(QObject* parent)
: QObject(parent)
{
    // Default scales
    m_scaleMap[FIFFV_MEG_CH] = kScaleMEGGrad;
    m_scaleMap[FIFFV_EEG_CH] = kScaleEEG;
    m_scaleMap[FIFFV_EOG_CH] = kScaleEOG;
    m_scaleMap[FIFFV_EMG_CH] = kScaleEMG;
    m_scaleMap[FIFFV_ECG_CH] = kScaleECG;
    m_scaleMap[FIFFV_STIM_CH] = kScaleSTIM;
    m_scaleMap[FIFFV_MISC_CH] = kScaleMISC;
}
 
//=============================================================================================================
 
void ChannelDataModel::init(QSharedPointer<FiffInfo> pFiffInfo)
{
    QWriteLocker lk(&m_lock);
    m_pFiffInfo = pFiffInfo;
 
    int nCh = (pFiffInfo ? pFiffInfo->nchan : 0) + m_virtualDisplayInfo.size();
    m_channelData.resize(nCh);
    for (auto& ch : m_channelData)
        ch.clear();
    m_firstSample = 0;
 
    lk.unlock();
    rebuildDisplayInfo();
    emit metaChanged();
}
 
//=============================================================================================================
 
void ChannelDataModel::setData(const MatrixXd& data, int firstSample)
{
    if (data.rows() == 0 || data.cols() == 0)
        return;
 
    QWriteLocker lk(&m_lock);
    int nCh = static_cast<int>(data.rows());
    m_channelData.resize(nCh);
    for (int ch = 0; ch < nCh; ++ch) {
        m_channelData[ch].resize(static_cast<int>(data.cols()));
        for (int s = 0; s < static_cast<int>(data.cols()); ++s)
            m_channelData[ch][s] = static_cast<float>(data(ch, s));
    }
    m_firstSample = firstSample;
    lk.unlock();
 
    emit dataChanged();
}
 
//=============================================================================================================
 
void ChannelDataModel::appendData(const MatrixXd& data)
{
    if (data.rows() == 0 || data.cols() == 0)
        return;
 
    QWriteLocker lk(&m_lock);
    int nCh = static_cast<int>(data.rows());
    int nNew = static_cast<int>(data.cols());
 
    if (m_channelData.size() != nCh)
        m_channelData.resize(nCh);
 
    for (int ch = 0; ch < nCh; ++ch) {
        auto& buf = m_channelData[ch];
        int oldSize = buf.size();
        buf.resize(oldSize + nNew);
        for (int s = 0; s < nNew; ++s)
            buf[oldSize + s] = static_cast<float>(data(ch, s));
 
        // Drop oldest samples when capacity exceeded
        if (m_maxStoredSamples > 0 && buf.size() > m_maxStoredSamples) {
            int drop = buf.size() - m_maxStoredSamples;
            buf.remove(0, drop);
            if (ch == 0)
                m_firstSample += drop;
        }
    }
    lk.unlock();
 
    emit dataChanged();
}
 
//=============================================================================================================
 
void ChannelDataModel::clearData()
{
    {
        QWriteLocker lk(&m_lock);
        for (auto& ch : m_channelData)
            ch.clear();
        m_firstSample = 0;
    }
    emit dataChanged();
}
 
//=============================================================================================================
 
void ChannelDataModel::setScaleMap(const QMap<qint32, float>& scaleMap)
{
    {
        QWriteLocker lk(&m_lock);
        m_scaleMap = scaleMap;
    }
    rebuildDisplayInfo();
    emit metaChanged();
}
 
//=============================================================================================================
 
void ChannelDataModel::setScaleMapFromStrings(const QMap<QString, double>& scaleMap)
{
    QMap<qint32, float> intMap;
    if (scaleMap.contains(QStringLiteral("MEG_grad")))
        intMap[kMEGGradKind] = static_cast<float>(scaleMap.value(QStringLiteral("MEG_grad")));
    if (scaleMap.contains(QStringLiteral("MEG_mag")))
        intMap[kMEGMagKind] = static_cast<float>(scaleMap.value(QStringLiteral("MEG_mag")));
    if (scaleMap.contains(QStringLiteral("MEG_EEG")))
        intMap[FIFFV_EEG_CH] = static_cast<float>(scaleMap.value(QStringLiteral("MEG_EEG")));
    if (scaleMap.contains(QStringLiteral("MEG_EOG")))
        intMap[FIFFV_EOG_CH] = static_cast<float>(scaleMap.value(QStringLiteral("MEG_EOG")));
    if (scaleMap.contains(QStringLiteral("MEG_EMG")))
        intMap[FIFFV_EMG_CH] = static_cast<float>(scaleMap.value(QStringLiteral("MEG_EMG")));
    if (scaleMap.contains(QStringLiteral("MEG_ECG")))
        intMap[FIFFV_ECG_CH] = static_cast<float>(scaleMap.value(QStringLiteral("MEG_ECG")));
    if (scaleMap.contains(QStringLiteral("MEG_MISC")))
        intMap[FIFFV_MISC_CH] = static_cast<float>(scaleMap.value(QStringLiteral("MEG_MISC")));
    if (scaleMap.contains(QStringLiteral("MEG_STIM")))
        intMap[FIFFV_STIM_CH] = static_cast<float>(scaleMap.value(QStringLiteral("MEG_STIM")));
    setScaleMap(intMap);
}
 
//=============================================================================================================
 
void ChannelDataModel::setVirtualChannels(const QVector<ChannelDisplayInfo>& virtualChannels)
{
    {
        QWriteLocker lk(&m_lock);
        m_virtualDisplayInfo = virtualChannels;
 
        const int realChannelCount = m_pFiffInfo ? m_pFiffInfo->nchan : 0;
        const int totalChannelCount = realChannelCount + m_virtualDisplayInfo.size();
        m_channelData.resize(totalChannelCount);
    }
 
    rebuildDisplayInfo();
    emit metaChanged();
    emit dataChanged();
}
 
//=============================================================================================================
 
void ChannelDataModel::setSignalColor(const QColor& color)
{
    {
        QWriteLocker lk(&m_lock);
        m_signalColor = color;
    }
    rebuildDisplayInfo();
    emit metaChanged();
}
 
//=============================================================================================================
 
void ChannelDataModel::setMaxStoredSamples(int n)
{
    QWriteLocker lk(&m_lock);
    m_maxStoredSamples = (n > 0) ? n : 0;
}
 
//=============================================================================================================
 
void ChannelDataModel::setRemoveDC(bool remove)
{
    setDetrendMode(remove ? DetrendMode::Mean : DetrendMode::None);
}
 
//=============================================================================================================
 
void ChannelDataModel::setDetrendMode(DetrendMode mode)
{
    {
        QWriteLocker lk(&m_lock);
        m_detrendMode = mode;
    }
    emit dataChanged();
}
 
//=============================================================================================================
 
void ChannelDataModel::setChannelBad(int channelIdx, bool bad)
{
    QWriteLocker lk(&m_lock);
    if (channelIdx >= 0 && channelIdx < m_displayInfo.size())
        m_displayInfo[channelIdx].bad = bad;
 
    if (m_pFiffInfo && channelIdx >= 0 && channelIdx < m_pFiffInfo->nchan) {
        const QString channelName = m_pFiffInfo->ch_names.value(channelIdx);
        const int badIndex = m_pFiffInfo->bads.indexOf(channelName);
 
        if (bad) {
            if (badIndex < 0)
                m_pFiffInfo->bads.append(channelName);
        } else if (badIndex >= 0) {
            m_pFiffInfo->bads.removeAt(badIndex);
        }
    }
 
    lk.unlock();
    emit metaChanged();
}
 
//=============================================================================================================
 
int ChannelDataModel::channelCount() const
{
    QReadLocker lk(&m_lock);
    const int realChannelCount = m_pFiffInfo ? m_pFiffInfo->nchan : 0;
    return qMax(m_channelData.size(), realChannelCount + m_virtualDisplayInfo.size());
}
 
//=============================================================================================================
 
int ChannelDataModel::firstSample() const
{
    QReadLocker lk(&m_lock);
    return m_firstSample;
}
 
//=============================================================================================================
 
int ChannelDataModel::totalSamples() const
{
    QReadLocker lk(&m_lock);
    return m_channelData.isEmpty() ? 0 : m_channelData[0].size();
}
 
//=============================================================================================================
 
ChannelDisplayInfo ChannelDataModel::channelInfo(int channelIdx) const
{
    QReadLocker lk(&m_lock);
    if (channelIdx >= 0 && channelIdx < m_displayInfo.size())
        return m_displayInfo[channelIdx];
    return {};
}
 
//=============================================================================================================
 
float ChannelDataModel::channelRms(int channelIdx, int first, int last) const
{
    QReadLocker lk(&m_lock);
    if (channelIdx < 0 || channelIdx >= m_channelData.size())
        return 0.f;
    const QVector<float>& src = m_channelData[channelIdx];
    int bufFirst = qBound(0, first - m_firstSample, src.size());
    int bufLast = qBound(0, last - m_firstSample, src.size());
    if (bufLast <= bufFirst)
        return 0.f;
    // Cap at 1000 samples: use the last kMax samples of the window for speed
    constexpr int kMax = 1000;
    if (bufLast - bufFirst > kMax)
        bufFirst = bufLast - kMax;
    double sum = 0.0;
    for (int i = bufFirst; i < bufLast; ++i)
        sum += static_cast<double>(src[i]) * src[i];
    return static_cast<float>(qSqrt(sum / (bufLast - bufFirst)));
}
 
//=============================================================================================================
 
float ChannelDataModel::sampleValueAt(int channelIdx, int sample) const
{
    QReadLocker lk(&m_lock);
    if (channelIdx < 0 || channelIdx >= m_channelData.size())
        return 0.f;
    const QVector<float>& src = m_channelData[channelIdx];
    int bufIdx = sample - m_firstSample;
    if (bufIdx < 0 || bufIdx >= src.size())
        return 0.f;
    return src[bufIdx];
}
 
//=============================================================================================================
 
QVector<float> ChannelDataModel::decimatedVertices(int channelIdx,
                                                   int firstSample,
                                                   int lastSample,
                                                   int pixelWidth,
                                                   int& vboFirstSample) const
{
    QReadLocker lk(&m_lock);
 
    vboFirstSample = firstSample; // may be updated below after clamping
 
    if (channelIdx < 0 || channelIdx >= m_channelData.size() || pixelWidth <= 0 || firstSample >= lastSample) {
        return {};
    }
 
    const QVector<float>& src = m_channelData[channelIdx];
    // Map absolute sample indices to buffer indices
    int bufFirst = firstSample - m_firstSample;
    int bufLast = lastSample - m_firstSample;
 
    // Clamp to available data
    bufFirst = qBound(0, bufFirst, src.size());
    bufLast = qBound(0, bufLast, src.size());
 
    // Update to the actual absolute sample index at vertex 0 after clamping.
    // When firstSample < m_firstSample (tile extends before buffer start),
    // bufFirst is clamped to 0 so vboFirstSample must reflect m_firstSample,
    // not the original (possibly negative) firstSample.
    vboFirstSample = m_firstSample + bufFirst;
 
    if (bufLast <= bufFirst)
        return {};
 
    int nSamples = bufLast - bufFirst;
 
    // ── Detrending: compute offset / trend over the window ────────────
    float dcOffset = 0.f;
    float linearSlope = 0.f;
    float linearIntercept = 0.f;
    const bool useLinear = (m_detrendMode == DetrendMode::Linear);
    const bool useMean = (m_detrendMode == DetrendMode::Mean);
 
    if (useMean) {
        double sum = 0.0;
        for (int i = bufFirst; i < bufLast; ++i)
            sum += src[i];
        dcOffset = static_cast<float>(sum / nSamples);
    } else if (useLinear) {
        // Least-squares fit: y = slope * t + intercept, where t = 0..nSamples-1
        double sumX = 0.0, sumY = 0.0, sumXX = 0.0, sumXY = 0.0;
        for (int i = 0; i < nSamples; ++i) {
            double x = static_cast<double>(i);
            double y = static_cast<double>(src[bufFirst + i]);
            sumX += x;
            sumY += y;
            sumXX += x * x;
            sumXY += x * y;
        }
        double denom = nSamples * sumXX - sumX * sumX;
        if (qAbs(denom) > 1e-30) {
            linearSlope = static_cast<float>((nSamples * sumXY - sumX * sumY) / denom);
            linearIntercept = static_cast<float>((sumY - linearSlope * sumX) / nSamples);
        }
    }
 
    if (nSamples <= pixelWidth * 2) {
        // ── Raw path: one vertex per sample ────────────────────────────
        QVector<float> result;
        result.reserve(nSamples * 2);
        for (int i = 0; i < nSamples; ++i) {
            float trend = useMean ? dcOffset
                : useLinear       ? (linearSlope * i + linearIntercept)
                                  : 0.f;
            result.append(static_cast<float>(i));
            result.append(src[bufFirst + i] - trend);
        }
        return result;
    }
 
    // ── Decimation path: min/max envelope, 2 vertices per pixel ─────────
    // Interleaved as (x_at_col, max), (x_at_col, min) per column.
    // Connected as LINE_STRIP this draws: vertical spike at each column,
    // diagonals between columns — the classic oscilloscope envelope look.
 
    QVector<float> result;
    result.reserve(pixelWidth * 4); // 2 vertices × 2 floats
 
    float spp = static_cast<float>(nSamples) / pixelWidth; // samples per pixel
 
    for (int px = 0; px < pixelWidth; ++px) {
        int sBegin = bufFirst + static_cast<int>(px * spp);
        int sEnd = bufFirst + static_cast<int>((px + 1) * spp);
        sEnd = qMin(sEnd, bufLast);
        if (sBegin >= sEnd)
            sBegin = qMax(sEnd - 1, bufFirst);
 
        float minV = src[sBegin];
        float maxV = src[sBegin];
        for (int s = sBegin + 1; s < sEnd; ++s) {
            if (src[s] < minV)
                minV = src[s];
            if (src[s] > maxV)
                maxV = src[s];
        }
 
        // Subtract trend at the center of this pixel bin
        float tCenter = static_cast<float>(sBegin - bufFirst) + (sEnd - sBegin) * 0.5f;
        float trend = useMean ? dcOffset
            : useLinear       ? (linearSlope * tCenter + linearIntercept)
                              : 0.f;
        minV -= trend;
        maxV -= trend;
 
        float xOffset = px * spp;
 
        // Emit max first, then min: draws ascending spike first
        result.append(xOffset);
        result.append(maxV);
        result.append(xOffset);
        result.append(minV);
    }
 
    return result;
}
 
//=============================================================================================================
// Private
//=============================================================================================================
 
void ChannelDataModel::rebuildDisplayInfo()
{
    QWriteLocker lk(&m_lock);
    const int realChannelCount = m_pFiffInfo ? m_pFiffInfo->nchan : 0;
    const int nCh = qMax(m_channelData.size(), realChannelCount + m_virtualDisplayInfo.size());
    m_displayInfo.resize(nCh);
    for (int ch = 0; ch < nCh; ++ch) {
        if (ch >= realChannelCount && ch - realChannelCount < m_virtualDisplayInfo.size()) {
            ChannelDisplayInfo info = m_virtualDisplayInfo.at(ch - realChannelCount);
            if (info.name.isEmpty())
                info.name = QString("Virtual %1").arg(ch - realChannelCount + 1);
            if (info.typeLabel.isEmpty())
                info.typeLabel = QStringLiteral("MISC");
            if (!info.color.isValid())
                info.color = m_signalColor;
            if (info.amplitudeMax <= 0.f)
                info.amplitudeMax = kScaleFallback;
            m_displayInfo[ch] = info;
            m_displayInfo[ch].isVirtualChannel = true;
            continue;
        }
 
        m_displayInfo[ch].amplitudeMax = amplitudeMaxForChannel(ch);
        m_displayInfo[ch].color = colorForChannel(ch);
        if (m_pFiffInfo && ch < realChannelCount)
            m_displayInfo[ch].name = m_pFiffInfo->ch_names[ch];
        else
            m_displayInfo[ch].name = QString("CH %1").arg(ch + 1);
        m_displayInfo[ch].typeLabel = typeLabelForChannel(ch);
        m_displayInfo[ch].bad = (m_pFiffInfo && ch < realChannelCount)
            ? m_pFiffInfo->bads.contains(m_displayInfo[ch].name)
            : false;
        m_displayInfo[ch].isVirtualChannel = false;
    }
}
 
//=============================================================================================================
 
float ChannelDataModel::amplitudeMaxForChannel(int ch) const
{
    if (!m_pFiffInfo || ch >= m_pFiffInfo->nchan)
        return kScaleFallback;
 
    const auto& info = m_pFiffInfo->chs[ch];
    qint32 kind = info.kind;
 
    // MEG: distinguish gradiometer vs. magnetometer by unit
    if (kind == FIFFV_MEG_CH) {
        if (info.unit == FIFF_UNIT_T_M)
            return m_scaleMap.value(kMEGGradKind,
                                    m_scaleMap.value(FIFFV_MEG_CH, kScaleMEGGrad));
        else
            return m_scaleMap.value(kMEGMagKind,
                                    m_scaleMap.value(FIFFV_MEG_CH, kScaleMEGMag));
    }
    if (m_scaleMap.contains(kind))
        return m_scaleMap.value(kind);
    return kScaleFallback;
}
 
//=============================================================================================================
 
QColor ChannelDataModel::colorForChannel(int ch) const
{
    if (!m_pFiffInfo || ch >= m_pFiffInfo->nchan)
        return m_signalColor;
 
    // Dark colours chosen to be readable on a light (near-white) background
    switch (m_pFiffInfo->chs[ch].kind) {
        case FIFFV_MEG_CH:
            return QColor(20, 90, 180); // dark blue for MEG
        case FIFFV_EEG_CH:
            return QColor(170, 55, 10); // dark orange for EEG
        case FIFFV_EOG_CH:
            return QColor(130, 0, 130); // dark purple for EOG
        case FIFFV_ECG_CH:
            return QColor(190, 15, 45); // dark crimson for ECG
        case FIFFV_EMG_CH:
            return QColor(20, 110, 20); // dark green for EMG
        case FIFFV_STIM_CH:
            return QColor(180, 100, 0); // dark amber for STIM
        default:
            return m_signalColor;
    }
}
 
//=============================================================================================================
 
QString ChannelDataModel::typeLabelForChannel(int ch) const
{
    if (!m_pFiffInfo || ch >= m_pFiffInfo->nchan)
        return QStringLiteral("MISC");
    switch (m_pFiffInfo->chs[ch].kind) {
        case FIFFV_MEG_CH:
            if (m_pFiffInfo->chs[ch].unit == FIFF_UNIT_T_M)
                return QStringLiteral("MEG grad");
            return QStringLiteral("MEG mag");
        case FIFFV_EEG_CH:
            return QStringLiteral("EEG");
        case FIFFV_EOG_CH:
            return QStringLiteral("EOG");
        case FIFFV_ECG_CH:
            return QStringLiteral("ECG");
        case FIFFV_EMG_CH:
            return QStringLiteral("EMG");
        case FIFFV_STIM_CH:
            return QStringLiteral("STIM");
        default:
            return QStringLiteral("MISC");
    }
}
 
//=============================================================================================================
 
float ChannelDataModel::sfreq() const
{
    QReadLocker lk(&m_lock);
    return m_pFiffInfo ? static_cast<float>(m_pFiffInfo->sfreq) : 0.f;
}