rtsensorinterpolationmatworker.cpp

Go to the documentation of this file.

//=============================================================================================================
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "rtsensorinterpolationmatworker.h"
#include "../scene/sensorfieldmapper.h"
#include <fwd/fwd_field_map.h>
 
#include <fiff/fiff_ch_info.h>
#include <fiff/fiff_constants.h>
#include <fiff/fiff_coord_trans.h>
#include <fwd/fwd_coil_set.h>
#include <fs/fs_surface.h>
 
#include <QCoreApplication>
#include <QDebug>
#include <QVector3D>
#include <QMatrix4x4>
#include <cmath>
 
using namespace FIFFLIB;
using namespace DISP3DLIB;
 
//=============================================================================================================
// ANONYMOUS HELPERS
//=============================================================================================================
 
namespace
{

Eigen::Vector3f applyTransform(const Eigen::Vector3f& point,
                               const FiffCoordTrans& trans)
{
    if (trans.isEmpty())
        return point;
    float r[3] = {point.x(), point.y(), point.z()};
    FiffCoordTrans::apply_trans(r, trans, FIFFV_MOVE);
    return Eigen::Vector3f(r[0], r[1], r[2]);
}

QMatrix4x4 toQMatrix4x4(const Eigen::Matrix4f& m)
{
    QMatrix4x4 q;
    for (int r = 0; r < 4; ++r)
        for (int c = 0; c < 4; ++c)
            q(r, c) = m(r, c);
    return q;
}
 
} // anonymous namespace
 
//=============================================================================================================
// MEMBER METHODS
//=============================================================================================================
 
RtSensorInterpolationMatWorker::RtSensorInterpolationMatWorker(QObject* parent)
: QObject(parent)
{
}
 
//=============================================================================================================
 
void RtSensorInterpolationMatWorker::setEvoked(const FiffEvoked& evoked)
{
    QMutexLocker locker(&m_mutex);
    m_evoked = evoked;
    m_hasEvoked = (evoked.data.rows() > 0 && evoked.data.cols() > 0);
}
 
//=============================================================================================================
 
void RtSensorInterpolationMatWorker::setTransform(const FiffCoordTrans& trans,
                                                  bool applySensorTrans)
{
    QMutexLocker locker(&m_mutex);
    m_headToMriTrans = trans;
    m_applySensorTrans = applySensorTrans;
}
 
//=============================================================================================================
 
void RtSensorInterpolationMatWorker::setMegFieldMapOnHead(bool onHead)
{
    QMutexLocker locker(&m_mutex);
    m_megOnHead = onHead;
}
 
//=============================================================================================================
 
void RtSensorInterpolationMatWorker::setMegSurface(const QString& surfaceKey,
                                                   const Eigen::MatrixX3f& vertices,
                                                   const Eigen::MatrixX3f& normals,
                                                   const Eigen::MatrixX3i& triangles)
{
    QMutexLocker locker(&m_mutex);
    m_megSurfaceKey = surfaceKey;
    m_megVertices = vertices;
    m_megNormals = normals;
    m_megTriangles = triangles;
    m_hasMegSurface = (vertices.rows() > 0);
}
 
//=============================================================================================================
 
void RtSensorInterpolationMatWorker::setEegSurface(const QString& surfaceKey,
                                                   const Eigen::MatrixX3f& vertices)
{
    QMutexLocker locker(&m_mutex);
    m_eegSurfaceKey = surfaceKey;
    m_eegVertices = vertices;
    m_hasEegSurface = (vertices.rows() > 0);
}
 
//=============================================================================================================
 
void RtSensorInterpolationMatWorker::setBadChannels(const QStringList& bads)
{
    QMutexLocker locker(&m_mutex);
    m_bads = bads;
}
 
//=============================================================================================================
 
void RtSensorInterpolationMatWorker::computeMapping()
{
    // ── Snapshot parameters under lock ─────────────────────────────────
    FiffEvoked evoked;
    FiffCoordTrans headToMriTrans;
    bool applySensorTrans;
    bool megOnHead;
    QString megSurfaceKey, eegSurfaceKey;
    Eigen::MatrixX3f megVerts, megNorms, eegVerts;
    Eigen::MatrixX3i megTris;
    bool hasMegSurface, hasEegSurface, hasEvoked;
    QStringList bads;
 
    {
        QMutexLocker locker(&m_mutex);
        evoked = m_evoked;
        hasEvoked = m_hasEvoked;
        headToMriTrans = m_headToMriTrans;
        applySensorTrans = m_applySensorTrans;
        megOnHead = m_megOnHead;
        megSurfaceKey = m_megSurfaceKey;
        eegSurfaceKey = m_eegSurfaceKey;
        megVerts = m_megVertices;
        megNorms = m_megNormals;
        megTris = m_megTriangles;
        hasMegSurface = m_hasMegSurface;
        eegVerts = m_eegVertices;
        hasEegSurface = m_hasEegSurface;
        bads = m_bads;
    }
 
    if (!hasEvoked) {
        qDebug() << "RtSensorInterpolationMatWorker: No evoked data set, skipping.";
        return;
    }
 
    qDebug() << "RtSensorInterpolationMatWorker: Computing mapping matrices...";
 
    // ── Build coordinate transforms ────────────────────────────────────
    bool hasDevHead = false;
    QMatrix4x4 devHeadQt;
    if (!evoked.info.dev_head_t.isEmpty() &&
        evoked.info.dev_head_t.from == FIFFV_COORD_DEVICE &&
        evoked.info.dev_head_t.to == FIFFV_COORD_HEAD &&
        !evoked.info.dev_head_t.trans.isIdentity()) {
        hasDevHead = true;
        devHeadQt = toQMatrix4x4(evoked.info.dev_head_t.trans);
    }
 
    QMatrix4x4 headToMri;
    if (applySensorTrans && !headToMriTrans.isEmpty()) {
        headToMri = toQMatrix4x4(headToMriTrans.trans);
    }
 
    // ── Classify channels ──────────────────────────────────────────────
    QList<FiffChInfo> megChs, eegChs;
    QVector<int> megPick, eegPick;
 
    for (int k = 0; k < evoked.info.chs.size(); ++k) {
        const auto& ch = evoked.info.chs[k];
        if (bads.contains(ch.ch_name) || evoked.info.bads.contains(ch.ch_name))
            continue;
 
        QVector3D pos(ch.chpos.r0(0), ch.chpos.r0(1), ch.chpos.r0(2));
 
        if (ch.kind == FIFFV_MEG_CH) {
            if (hasDevHead)
                pos = devHeadQt.map(pos);
            if (applySensorTrans && !headToMriTrans.isEmpty())
                pos = headToMri.map(pos);
            megPick.append(k);
            megChs.append(ch);
        } else if (ch.kind == FIFFV_EEG_CH) {
            if (applySensorTrans && !headToMriTrans.isEmpty())
                pos = headToMri.map(pos);
            eegPick.append(k);
            eegChs.append(ch);
        }
    }
 
    // ── Constants (matching MNE-Python) ────────────────────────────────
    constexpr float kIntrad = 0.06f;
    constexpr float kMegMiss = 1e-4f;
    constexpr float kEegMiss = 1e-3f;
    // origin="auto" of MNE-Python's make_field_map, as in SensorFieldMapper::buildMapping()
    const Eigen::Vector3f fittedOrigin = SensorFieldMapper::fitSphereOrigin(evoked.info);
 
    FiffCoordTrans headMri = (applySensorTrans && !headToMriTrans.isEmpty())
        ? headToMriTrans
        : FiffCoordTrans();
    FiffCoordTrans devHead = (!evoked.info.dev_head_t.isEmpty() &&
                              evoked.info.dev_head_t.from == FIFFV_COORD_DEVICE &&
                              evoked.info.dev_head_t.to == FIFFV_COORD_HEAD)
        ? evoked.info.dev_head_t
        : FiffCoordTrans();
 
    // ── MEG mapping ────────────────────────────────────────────────────
    if (hasMegSurface && !megChs.isEmpty()) {
        Eigen::MatrixX3f norms = megNorms;
 
        // Recompute normals if missing
        if (norms.rows() != megVerts.rows()) {
            if (megTris.rows() > 0) {
                norms = FSLIB::FsSurface::compute_normals(megVerts, megTris);
            }
        }
 
        if (megVerts.rows() > 0 && norms.rows() == megVerts.rows()) {
            const QString coilPath = QCoreApplication::applicationDirPath() + "/../resources/general/coilDefinitions/coil_def.dat";
            auto templates =
                FWDLIB::FwdCoilSet::read_coil_defs(coilPath);
 
            if (templates) {
                FiffCoordTrans devToTarget;
                if (megOnHead && !headMri.isEmpty()) {
                    if (!devHead.isEmpty()) {
                        devToTarget = FiffCoordTrans::combine(
                            FIFFV_COORD_DEVICE, FIFFV_COORD_MRI,
                            devHead, headMri);
                    }
                } else if (!devHead.isEmpty()) {
                    devToTarget = devHead;
                }
 
                Eigen::Vector3f origin = fittedOrigin;
                if (megOnHead && !headMri.isEmpty()) {
                    origin = applyTransform(origin, headMri);
                }
 
                auto coils = templates->create_meg_coils(
                    megChs, megChs.size(), FWDLIB::FWD_COIL_ACCURACY_NORMAL, devToTarget);
 
                if (coils && coils->ncoil() > 0) {
                    auto mat = FWDLIB::FwdFieldMap::computeMegMapping(
                        *coils, megVerts, norms, origin, kIntrad, kMegMiss);
 
                    if (mat && mat->rows() > 0) {
                        qDebug() << "RtSensorInterpolationMatWorker: MEG mapping computed:"
                                 << mat->rows() << "x" << mat->cols();
                        emit newMegMappingAvailable(megSurfaceKey,
                                                    std::shared_ptr<Eigen::MatrixXf>(std::move(mat)), megPick);
                    }
                }
            }
        }
    }
 
    // ── EEG mapping ────────────────────────────────────────────────────
    if (hasEegSurface && !eegChs.isEmpty()) {
        if (eegVerts.rows() > 0) {
            Eigen::Vector3f origin = fittedOrigin;
            if (!headMri.isEmpty())
                origin = applyTransform(origin, headMri);
 
            auto eegCoils =
                FWDLIB::FwdCoilSet::create_eeg_els(
                    eegChs, eegChs.size(), headMri);
 
            if (eegCoils && eegCoils->ncoil() > 0) {
                auto mat = FWDLIB::FwdFieldMap::computeEegMapping(
                    *eegCoils, eegVerts, origin, kIntrad, kEegMiss);
 
                if (mat && mat->rows() > 0) {
                    qDebug() << "RtSensorInterpolationMatWorker: EEG mapping computed:"
                             << mat->rows() << "x" << mat->cols();
                    emit newEegMappingAvailable(eegSurfaceKey,
                                                std::shared_ptr<Eigen::MatrixXf>(std::move(mat)), eegPick);
                }
            }
        }
    }
 
    qDebug() << "RtSensorInterpolationMatWorker: Mapping computation complete.";
}