sourceestimateoverlay.cpp

Go to the documentation of this file.

//=============================================================================================================
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "sourceestimateoverlay.h"
#include "renderable/brainsurface.h"
#include "core/rendertypes.h"
 
#include <disp/plots/helpers/colormap.h>
#include "helpers/interpolation.h"
#include "helpers/geometryinfo.h"
 
#include <QFile>
#include <QDebug>
#include <cmath>
 
//=============================================================================================================
// DEFINE MEMBER METHODS
//=============================================================================================================
 
SourceEstimateOverlay::SourceEstimateOverlay()
{
}
 
//=============================================================================================================
 
SourceEstimateOverlay::~SourceEstimateOverlay()
{
}
 
//=============================================================================================================
 
bool SourceEstimateOverlay::loadStc(const QString &path, int hemi)
{
    QFile file(path);
    // Note: InvSourceEstimate::read() opens the file internally, don't open it here
 
    INVLIB::InvSourceEstimate stc;
    if (!INVLIB::InvSourceEstimate::read(file, stc)) {
        qWarning() << "SourceEstimateOverlay::loadStc - Failed to read STC file:" << path;
        return false;
    }
 
    if (hemi == 0) {
        m_stcLh = stc;
        m_hasLh = true;
        qDebug() << "SourceEstimateOverlay: Loaded LH with" << stc.data.rows() << "vertices," 
                 << stc.data.cols() << "time points";
    } else {
        m_stcRh = stc;
        m_hasRh = true;
        qDebug() << "SourceEstimateOverlay: Loaded RH with" << stc.data.rows() << "vertices," 
                 << stc.data.cols() << "time points";
    }
    invalidateColorCache();
 
    // Auto-set thresholds based on data range
    if (m_hasLh || m_hasRh) {
        double minVal, maxVal;
        getDataRange(minVal, maxVal);
        m_threshMin = minVal;
        m_threshMax = maxVal;
        m_threshMid = (minVal + maxVal) / 2.0;
        qDebug() << "SourceEstimateOverlay: Auto thresholds set to" << m_threshMin << m_threshMid << m_threshMax;
    }
 
    return true;
}
 
//=============================================================================================================
 
bool SourceEstimateOverlay::isLoaded() const
{
    return m_hasLh || m_hasRh;
}
 
//=============================================================================================================
 
void SourceEstimateOverlay::applyToSurface(BrainSurface *surface, int timeIndex)
{
    if (!surface) return;
 
    int hemi = surface->hemi();
    const INVLIB::InvSourceEstimate *stc = nullptr;
    QSharedPointer<Eigen::SparseMatrix<float>> interpMat;
 
    if (hemi == 0 && m_hasLh) {
        stc = &m_stcLh;
        interpMat = m_interpolationMatLh;
    } else if (hemi == 1 && m_hasRh) {
        stc = &m_stcRh;
        interpMat = m_interpolationMatRh;
    } else {
        return; // No data for this hemisphere
    }
 
    if (stc->isEmpty()) return;
 
    // Clamp time index
    int tIdx = qBound(0, timeIndex, static_cast<int>(stc->data.cols()) - 1);
    const uint32_t vertexCount = surface->vertexCount();
 
    // ── Color cache fast-path ──────────────────────────────────────
    // Hot loop while scrubbing or auto-looping: avoid the sparse
    // interpolation matvec and per-vertex colormap evaluation by
    // returning the pre-computed buffer for this (hemi, vertexCount,
    // timeIndex). The cache is invalidated whenever colormap,
    // thresholds or source data change.
    const ColorCacheKey cacheKey(hemi, static_cast<int>(vertexCount));
    auto bucketIt = m_colorCache.find(cacheKey);
    if (bucketIt != m_colorCache.end()) {
        auto frameIt = bucketIt.value().constFind(tIdx);
        if (frameIt != bucketIt.value().constEnd()) {
            surface->applySourceEstimateColors(frameIt.value());
            return;
        }
    }
 
    // Get source data for this time point
    Eigen::VectorXf sourceData = stc->data.col(tIdx).cwiseAbs().cast<float>();
 
    // Create color array for all surface vertices
    QVector<uint32_t> colors(vertexCount, 0xFF808080); // Default gray
 
    // Determine if we have an interpolation matrix
    Eigen::VectorXf interpolatedData;
    
    if (interpMat && interpMat->rows() == static_cast<int>(vertexCount) && 
        interpMat->cols() == sourceData.size()) {
        // Use interpolation to spread values to all vertices
        // Note: interpolateSignal returns by value, not QSharedPointer, so assignment matches
        interpolatedData = DISP3DLIB::Interpolation::interpolateSignal(interpMat, QSharedPointer<Eigen::VectorXf>::create(sourceData));
    } else {
        // Fall back to sparse visualization (direct mapping)
        interpolatedData = Eigen::VectorXf::Zero(vertexCount);
        const Eigen::VectorXi &srcVertices = stc->vertices;
        for (int i = 0; i < srcVertices.size() && i < sourceData.size(); ++i) {
            int vertIdx = srcVertices(i);
            if (vertIdx >= 0 && vertIdx < static_cast<int>(vertexCount)) {
                interpolatedData(vertIdx) = sourceData(i);
            }
        }
    }
 
    // Convert interpolated values to colors
    for (int i = 0; i < static_cast<int>(vertexCount); ++i) {
        float value = interpolatedData(i);
 
        // Normalize based on thresholds
        double normalized = 0.0;
        if (m_threshMax > m_threshMin) {
            normalized = (value - m_threshMin) / (m_threshMax - m_threshMin);
            normalized = qBound(0.0, normalized, 1.0);
        }
 
        // Calculate alpha based on threshold
        uint8_t alpha = 255;
        if (value < m_threshMin) {
            alpha = 0; // Fully transparent below minimum
        } else if (value < m_threshMid) {
            // Fade in from min to mid
            float range = m_threshMid - m_threshMin;
            if (range > 0) {
                alpha = static_cast<uint8_t>(255.0f * (value - m_threshMin) / range);
            }
        }
 
        colors[i] = valueToColor(normalized, alpha);
    }
 
    // Store in cache before handing off to the surface so back-scrubs
    // and auto-loop iterations are free.
    m_colorCache[cacheKey].insert(tIdx, colors);
 
    surface->applySourceEstimateColors(colors);
}
 
//=============================================================================================================
 
void SourceEstimateOverlay::setColormap(const QString &name)
{
    if (m_colormap == name) return;
    m_colormap = name;
    invalidateColorCache();
}
 
//=============================================================================================================
 
void SourceEstimateOverlay::setThresholds(float min, float mid, float max)
{
    if (m_threshMin == min && m_threshMid == mid && m_threshMax == max) return;
    m_threshMin = min;
    m_threshMid = mid;
    m_threshMax = max;
    invalidateColorCache();
}
 
//=============================================================================================================
 
int SourceEstimateOverlay::numTimePoints() const
{
    if (m_hasLh) return m_stcLh.data.cols();
    if (m_hasRh) return m_stcRh.data.cols();
    return 0;
}
 
//=============================================================================================================
 
float SourceEstimateOverlay::timeAtIndex(int idx) const
{
    if (m_hasLh && idx < m_stcLh.times.size()) {
        return m_stcLh.times(idx);
    }
    if (m_hasRh && idx < m_stcRh.times.size()) {
        return m_stcRh.times(idx);
    }
    return 0.0f;
}
 
//=============================================================================================================
 
float SourceEstimateOverlay::tmin() const
{
    if (m_hasLh) return m_stcLh.tmin;
    if (m_hasRh) return m_stcRh.tmin;
    return 0.0f;
}
 
//=============================================================================================================
 
float SourceEstimateOverlay::tstep() const
{
    if (m_hasLh) return m_stcLh.tstep;
    if (m_hasRh) return m_stcRh.tstep;
    return 0.0f;
}
 
//=============================================================================================================
 
void SourceEstimateOverlay::getDataRange(double &minVal, double &maxVal) const
{
    minVal = std::numeric_limits<double>::max();
    maxVal = std::numeric_limits<double>::lowest();
 
    if (m_hasLh) {
        double lhMin = m_stcLh.data.minCoeff();
        double lhMax = m_stcLh.data.maxCoeff();
        minVal = qMin(minVal, std::abs(lhMin));
        maxVal = qMax(maxVal, std::abs(lhMax));
    }
 
    if (m_hasRh) {
        double rhMin = m_stcRh.data.minCoeff();
        double rhMax = m_stcRh.data.maxCoeff();
        minVal = qMin(minVal, std::abs(rhMin));
        maxVal = qMax(maxVal, std::abs(rhMax));
    }
 
    // If no data, set defaults
    if (minVal > maxVal) {
        minVal = 0.0;
        maxVal = 1.0;
    }
}
 
//=============================================================================================================
 
uint32_t SourceEstimateOverlay::valueToColor(double value, uint8_t alpha) const
{
    QRgb rgb = DISPLIB::ColorMap::valueToColor(value, m_colormap);
 
    uint32_t r = qRed(rgb);
    uint32_t g = qGreen(rgb);
    uint32_t b = qBlue(rgb);
 
    // Pack as ABGR (same format as BrainSurface uses)
    return packABGR(r, g, b, static_cast<uint32_t>(alpha));
}
 
//=============================================================================================================
 
void SourceEstimateOverlay::computeInterpolationMatrix(BrainSurface *surface, int hemi, double cancelDist)
{
    if (!surface) return;
 
    const INVLIB::InvSourceEstimate *stc = nullptr;
    QSharedPointer<Eigen::SparseMatrix<float>> *pMatPtr = nullptr;
 
    if (hemi == 0 && m_hasLh) {
        stc = &m_stcLh;
        pMatPtr = &m_interpolationMatLh;
    } else if (hemi == 1 && m_hasRh) {
        stc = &m_stcRh;
        pMatPtr = &m_interpolationMatRh;
    } else {
        return;
    }
 
    if (stc->isEmpty()) return;
 
    qDebug() << "SourceEstimateOverlay: Computing interpolation matrix for hemi" << hemi;
 
    // Get vertices and neighbor information needed for Dijkstra (SCDC)
    Eigen::MatrixX3f matVertices = surface->verticesAsMatrix();
    std::vector<Eigen::VectorXi> vecNeighbors = surface->computeNeighbors();
 
    // Source vertex subset from STC (already a VectorXi)
    Eigen::VectorXi vecSourceVertices = stc->vertices;
 
    qDebug() << "SourceEstimateOverlay: FsSurface has" << matVertices.rows() << "vertices,"
             << vecSourceVertices.size() << "sources";
    
    if (vecSourceVertices.size() == 0) {
        qWarning() << "SourceEstimateOverlay: No source vertices found";
        return;
    }
 
    // 1. Calculate Distance Table (Geodesic distance on surface)
    // This uses Dijkstra's algorithm via GeometryInfo::scdc
    // Note: This can be slow for many sources!
    qDebug() << "SourceEstimateOverlay: Computing distance table (SCDC)...";
    QSharedPointer<Eigen::MatrixXd> distTable = DISP3DLIB::GeometryInfo::scdc(
        matVertices,
        vecNeighbors,
        vecSourceVertices,
        cancelDist
    );
 
    if (!distTable || distTable->rows() == 0) {
        qWarning() << "SourceEstimateOverlay: Failed to compute distance table";
        return;
    }
 
    // 2. Create Interpolation Matrix
    qDebug() << "SourceEstimateOverlay: Creating interpolation matrix...";
    *pMatPtr = DISP3DLIB::Interpolation::createInterpolationMat(
        vecSourceVertices,
        distTable,
        DISP3DLIB::Interpolation::cubic,  // Use cubic interpolation function
        cancelDist
    );
 
    if (*pMatPtr && (*pMatPtr)->rows() > 0) {
        qDebug() << "SourceEstimateOverlay: Interpolation matrix created:" 
                 << (*pMatPtr)->rows() << "x" << (*pMatPtr)->cols();
    } else {
        qWarning() << "SourceEstimateOverlay: Failed to compute interpolation matrix";
    }
}
 
//=============================================================================================================
 
void SourceEstimateOverlay::setStcData(const INVLIB::InvSourceEstimate &stc, int hemi)
{
    if (hemi == 0) {
        m_stcLh = stc;
        m_hasLh = true;
    } else {
        m_stcRh = stc;
        m_hasRh = true;
    }
    invalidateColorCache();
}
 
//=============================================================================================================
 
void SourceEstimateOverlay::setInterpolationMatrix(QSharedPointer<Eigen::SparseMatrix<float>> mat, int hemi)
{
    if (hemi == 0) {
        m_interpolationMatLh = mat;
    } else {
        m_interpolationMatRh = mat;
    }
    invalidateColorCache();
}
 
//=============================================================================================================
 
void SourceEstimateOverlay::updateThresholdsFromData()
{
    if (m_hasLh || m_hasRh) {
        double minVal, maxVal;
        getDataRange(minVal, maxVal);
        m_threshMin = minVal;
        m_threshMax = maxVal;
        m_threshMid = (minVal + maxVal) / 2.0;
        invalidateColorCache();
        qDebug() << "SourceEstimateOverlay: Auto thresholds set to" << m_threshMin << m_threshMid << m_threshMax;
    }
}
 
//=============================================================================================================
 
Eigen::VectorXd SourceEstimateOverlay::sourceDataColumn(int timeIndex) const
{
    int nLh = m_hasLh ? m_stcLh.data.rows() : 0;
    int nRh = m_hasRh ? m_stcRh.data.rows() : 0;
 
    if (nLh == 0 && nRh == 0) {
        return Eigen::VectorXd();
    }
 
    Eigen::VectorXd result(nLh + nRh);
 
    if (m_hasLh) {
        int tIdx = qBound(0, timeIndex, static_cast<int>(m_stcLh.data.cols()) - 1);
        result.head(nLh) = m_stcLh.data.col(tIdx);
    }
 
    if (m_hasRh) {
        int tIdx = qBound(0, timeIndex, static_cast<int>(m_stcRh.data.cols()) - 1);
        result.segment(nLh, nRh) = m_stcRh.data.col(tIdx);
    }
 
    return result;
}