brainsurface.cpp

Go to the documentation of this file.

//=============================================================================================================
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "brainsurface.h"
 
#include <rhi/qrhi.h>
 
#include <set>
 
//=============================================================================================================
// PIMPL
//=============================================================================================================
 
struct BrainSurface::GpuBuffers
{
    std::unique_ptr<QRhiBuffer> vertexBuffer;
    std::unique_ptr<QRhiBuffer> indexBuffer;
    bool dirty = true;
    bool indexDirty = true;   // IBO needs (re-)upload (topology change)
};
 
//=============================================================================================================
 
void BrainSurface::markVertexDirty()
{
    m_gpu->dirty = true;
    ++m_vertexGeneration;
}
 
 
//=============================================================================================================
// DEFINE MEMBER METHODS
//=============================================================================================================
 
 
//=============================================================================================================
 
BrainSurface::BrainSurface()
    : m_gpu(std::make_unique<GpuBuffers>())
{
}
 
//=============================================================================================================
 
BrainSurface::~BrainSurface() = default;
 
//=============================================================================================================
 
QRhiBuffer* BrainSurface::vertexBuffer() const { return m_gpu->vertexBuffer.get(); }
QRhiBuffer* BrainSurface::indexBuffer()  const { return m_gpu->indexBuffer.get(); }
 
//=============================================================================================================
 
void BrainSurface::fromSurface(const FSLIB::FsSurface &surf)
{
    m_vertexData.clear();
    m_indexData.clear();
 
    const Eigen::MatrixXf &rr = surf.rr();
    const Eigen::MatrixXf &nn = surf.nn();
    const Eigen::MatrixXi &tris = surf.tris();
    const Eigen::VectorXf &curv = surf.curv();
    m_curvature.resize(curv.size());
    for(int i=0; i<curv.size(); ++i) m_curvature[i] = curv[i];
 
    // Populate vertex data
    m_vertexData.reserve(rr.rows());
    
    for (int i = 0; i < rr.rows(); ++i) {
        VertexData v;
        v.pos = QVector3D(rr(i, 0), rr(i, 1), rr(i, 2));
        v.norm = QVector3D(nn(i, 0), nn(i, 1), nn(i, 2));
        v.color = 0xFFFFFFFF; // Default white (overwritten by updateVertexColors)
        v.colorAnnotation = 0x00000000; // No annotation yet
        m_vertexData.append(v);
    }
 
    m_indexData.reserve(tris.rows() * 3);
    for (int i = 0; i < tris.rows(); ++i) {
        m_indexData.append(tris(i, 0));
        m_indexData.append(tris(i, 1));
        m_indexData.append(tris(i, 2));
    }
    m_indexCount = m_indexData.size();
    
    markVertexDirty();
    m_bAABBDirty = true;
    
    // Initial coloring based on current visualization mode
    updateVertexColors();
 
    m_originalVertexData = m_vertexData;
}
 
//=============================================================================================================
 
//=============================================================================================================
 
void BrainSurface::fromBemSurface(const MNELIB::MNEBemSurface &surf, const QColor &color)
{
    m_vertexData.clear();
    m_indexData.clear();
    m_curvature.clear(); // BEM has no curvature info usually
 
    int nVerts = surf.rr.rows();
    m_vertexData.reserve(nVerts);
    
    // Compute normals if missing
    Eigen::MatrixX3f nn = surf.nn;
    if (nn.rows() != nVerts) {
        nn = FSLIB::FsSurface::compute_normals(Eigen::MatrixX3f(surf.rr), Eigen::MatrixX3i(surf.itris));
    }
    
    m_defaultColor = color;
    m_baseColor = color;
    uint32_t colorVal = packABGR(color.red(), color.green(), color.blue(), color.alpha());
 
    for (int i = 0; i < nVerts; ++i) {
        VertexData v;
        v.pos = QVector3D(surf.rr(i, 0), surf.rr(i, 1), surf.rr(i, 2));
        v.norm = QVector3D(nn(i, 0), nn(i, 1), nn(i, 2));
        v.color = colorVal;
        v.colorAnnotation = 0x00000000;
        m_vertexData.append(v);
    }
 
    int nTris = surf.itris.rows();
    m_indexData.reserve(nTris * 3);
    for (int i = 0; i < nTris; ++i) {
        m_indexData.append(surf.itris(i, 0));
        m_indexData.append(surf.itris(i, 1));
        m_indexData.append(surf.itris(i, 2));
    }
    m_indexCount = m_indexData.size();
    
    m_originalVertexData = m_vertexData;
    markVertexDirty();
}
 
void BrainSurface::createFromData(const Eigen::MatrixX3f &vertices, const Eigen::MatrixX3i &triangles, const QColor &color)
{
    // Compute spherical normals (legacy behavior / fallback)
    // Note: detailed normals should be passed via the overload for specific shapes
    Eigen::MatrixX3f normals(vertices.rows(), 3);
    for(int i=0; i<vertices.rows(); ++i) {
        QVector3D p(vertices(i, 0), vertices(i, 1), vertices(i, 2));
        QVector3D n = p.normalized();
        normals(i, 0) = n.x();
        normals(i, 1) = n.y();
        normals(i, 2) = n.z();
    }
    createFromData(vertices, normals, triangles, color);
}
 
void BrainSurface::createFromData(const Eigen::MatrixX3f &vertices, const Eigen::MatrixX3f &normals, const Eigen::MatrixX3i &triangles, const QColor &color)
{
    m_vertexData.clear();
    m_indexData.clear();
    m_curvature.clear(); 
 
    int nVerts = vertices.rows();
    m_vertexData.reserve(nVerts);
    
    m_defaultColor = color;
    m_baseColor = color;
    uint32_t colorVal = packABGR(color.red(), color.green(), color.blue(), color.alpha());
 
    for (int i = 0; i < nVerts; ++i) {
        VertexData v;
        v.pos = QVector3D(vertices(i, 0), vertices(i, 1), vertices(i, 2));
        v.norm = QVector3D(normals(i, 0), normals(i, 1), normals(i, 2));
        v.color = colorVal;
        v.colorAnnotation = 0x00000000;
        m_vertexData.append(v);
    }
 
    int nTris = triangles.rows();
    m_indexData.reserve(nTris * 3);
    for (int i = 0; i < nTris; ++i) {
        m_indexData.append(triangles(i, 0));
        m_indexData.append(triangles(i, 1));
        m_indexData.append(triangles(i, 2));
    }
    m_indexCount = m_indexData.size();
    
    m_originalVertexData = m_vertexData;
    markVertexDirty();
}
 
//=============================================================================================================
 
Eigen::MatrixX3f BrainSurface::vertexPositions() const
{
    Eigen::MatrixX3f rr(m_vertexData.size(), 3);
    for (int i = 0; i < m_vertexData.size(); ++i) {
        rr(i, 0) = m_vertexData[i].pos.x();
        rr(i, 1) = m_vertexData[i].pos.y();
        rr(i, 2) = m_vertexData[i].pos.z();
    }
    return rr;
}
 
//=============================================================================================================
 
Eigen::MatrixX3f BrainSurface::vertexNormals() const
{
    Eigen::MatrixX3f nn(m_vertexData.size(), 3);
    for (int i = 0; i < m_vertexData.size(); ++i) {
        nn(i, 0) = m_vertexData[i].norm.x();
        nn(i, 1) = m_vertexData[i].norm.y();
        nn(i, 2) = m_vertexData[i].norm.z();
    }
    return nn;
}
 
//=============================================================================================================
 
bool BrainSurface::loadAnnotation(const QString &path)
{
    if (!FSLIB::FsAnnotation::read(path, m_annotation)) {
        qWarning() << "BrainSurface: Failed to load annotation from" << path;
        return false;
    }
    m_hasAnnotation = true;
    updateVertexColors();
    markVertexDirty();
    return true;
}
 
void BrainSurface::addAnnotation(const FSLIB::FsAnnotation &annotation)
{
    m_annotation = annotation;
    m_hasAnnotation = true;
    updateVertexColors();
    markVertexDirty();
}
 
 
//=============================================================================================================
 
void BrainSurface::setVisible(bool visible)
{
    m_visible = visible;
}
 
//=============================================================================================================
 
void BrainSurface::setVisualizationMode(VisualizationMode mode)
{
    if (m_visMode == mode)
        return;
 
    m_visMode = mode;
 
    // Scientific and SourceEstimate both read from the primary colour
    // channel (v_color).  When switching between them the vertex buffer
    // must be refreshed so the channel holds curvature grays (Scientific)
    // or STC colours (SourceEstimate).
    updateVertexColors();
    markVertexDirty();
}
 
//=============================================================================================================
 
void BrainSurface::applySourceEstimateColors(const QVector<uint32_t> &colors)
{
    m_visMode = ModeSourceEstimate;
    m_stcColors = colors;
    
    // Write STC colours into the primary color channel
    for (int i = 0; i < qMin(colors.size(), m_vertexData.size()); ++i) {
        m_vertexData[i].color = colors[i];
    }
    
    markVertexDirty();
}
 
//=============================================================================================================
 
void BrainSurface::clearSourceEstimateColors()
{
    m_stcColors.clear();
    m_visMode = ModeSurface;
    updateVertexColors();
    markVertexDirty();
}
 
//=============================================================================================================
 
void BrainSurface::setUseDefaultColor(bool useDefault)
{
    m_baseColor = useDefault ? m_defaultColor : Qt::white;
    updateVertexColors();
}
 
void BrainSurface::updateVertexColors()
{
    // ── 1. Populate the primary "color" channel.
    //       Brain surfaces (with curvature): curvature-based gray.
    //       Non-brain surfaces (BEM, sensors, etc.): use m_baseColor.
    //       The shader uses this for Scientific mode and as a lighting base;
    //       in FsSurface mode the shader overrides with white for brain tissue.
    if (!m_curvature.isEmpty() && m_curvature.size() == m_vertexData.size()) {
        for (int i = 0; i < m_vertexData.size(); ++i) {
            const uint32_t val = (m_curvature[i] > 0.0f) ? 0x40u : 0xAAu;
            m_vertexData[i].color = packABGR(val, val, val);
        }
    } else {
        // Use the base colour (set via setUseDefaultColor / fromBemSurface).
        // This preserves BEM Red/Green/Blue colours and sensor colours.
        const uint32_t baseVal = packABGR(m_baseColor.red(), m_baseColor.green(),
                                          m_baseColor.blue(), m_baseColor.alpha());
        for (int i = 0; i < m_vertexData.size(); ++i) {
            m_vertexData[i].color = baseVal;
        }
    }
 
    // Always keep STC colours in the primary colour channel when available.
    // Each viewport selects the overlay mode via a per-draw shader uniform
    // (overlayMode), so the vertex buffer must hold STC data for any
    // viewport that shows Source Estimate, regardless of this surface's
    // m_visMode.  FsAnnotation data lives in a separate vertex attribute
    // (colorAnnotation) and is unaffected.
    if (!m_stcColors.isEmpty()) {
        for (int i = 0; i < qMin(m_stcColors.size(), m_vertexData.size()); ++i) {
            m_vertexData[i].color = m_stcColors[i];
        }
    }
 
    // ── 2. Populate colorAnnotation from loaded annotation data.
    for (auto &v : m_vertexData) {
        v.colorAnnotation = 0x00000000;
    }
 
    if (m_hasAnnotation && !m_vertexData.isEmpty()) {
        const Eigen::VectorXi &vertices = m_annotation.getVertices();
        const Eigen::VectorXi &labelIds = m_annotation.getLabelIds();
        const FSLIB::FsColortable &ct = m_annotation.getColortable();
 
        for (int i = 0; i < labelIds.rows(); ++i) {
            int vertexIdx = vertices(i);
            if (vertexIdx >= 0 && vertexIdx < m_vertexData.size()) {
                int colorIdx = -1;
                for (int c = 0; c < ct.numEntries; ++c) {
                    if (ct.table(c, 4) == labelIds(i)) {
                        colorIdx = c;
                        break;
                    }
                }
                if (colorIdx >= 0) {
                    uint32_t r = ct.table(colorIdx, 0);
                    uint32_t g = ct.table(colorIdx, 1);
                    uint32_t b = ct.table(colorIdx, 2);
                    m_vertexData[vertexIdx].colorAnnotation =
                        packABGR(r, g, b);
                }
            }
        }
    }
 
    // ── 3. Selection highlighting.
    //       Tint the selected region / vertex range gold in the vertex
    //       buffer so the user gets precise per-region feedback.
    //       On WASM the merged-draw path re-reads vertexDataRef() every
    //       frame, so this is safe without separate buffer re-uploads.
    if (m_selected) {
        const uint32_t gold = packABGR(255, 200, 60);
        if (m_selectedRegionId != -1 && m_hasAnnotation) {
            // Highlight a specific annotation region
            const Eigen::VectorXi &vertices = m_annotation.getVertices();
            const Eigen::VectorXi &labelIds = m_annotation.getLabelIds();
            for (int i = 0; i < labelIds.rows(); ++i) {
                if (labelIds(i) == m_selectedRegionId) {
                    int idx = vertices(i);
                    if (idx >= 0 && idx < m_vertexData.size()) {
                        m_vertexData[idx].color = gold;
                        m_vertexData[idx].colorAnnotation = gold;
                    }
                }
            }
        } else if (m_selectedVertexStart >= 0 && m_selectedVertexCount > 0) {
            // Highlight a vertex range (e.g. single digitizer sphere)
            const int end = qMin(m_selectedVertexStart + m_selectedVertexCount,
                                 m_vertexData.size());
            for (int i = m_selectedVertexStart; i < end; ++i) {
                m_vertexData[i].color = gold;
            }
        }
        // Whole-surface selection (no region / vertex range) is handled
        // by the shader gold glow via the isSelected uniform.
    }
 
    markVertexDirty();
}
 
float BrainSurface::minX() const
{
    float minVal = std::numeric_limits<float>::max();
    for (const auto &v : m_vertexData) {
        if (v.pos.x() < minVal) minVal = v.pos.x();
    }
    return minVal;
}
 
float BrainSurface::maxX() const
{
    float maxVal = std::numeric_limits<float>::lowest();
    for (const auto &v : m_vertexData) {
        if (v.pos.x() > maxVal) maxVal = v.pos.x();
    }
    return maxVal;
}
 
void BrainSurface::translateX(float offset)
{
    for (auto &v : m_vertexData) {
        v.pos.setX(v.pos.x() + offset);
    }
    markVertexDirty();
    m_bAABBDirty = true;
}
 
//=============================================================================================================
 
void BrainSurface::transform(const QMatrix4x4 &m)
{
    // Extract 3x3 normal matrix (inverse transpose of upper-left 3x3)
    // QMatrix4x4::normalMatrix() returns QMatrix3x3.
    QMatrix3x3 normalMat = m.normalMatrix();
 
    for (auto &v : m_vertexData) {
        // Transform position
        v.pos = m.map(v.pos);
        
        // Transform normal
        // Note: QMatrix3x3 * QVector3D isn't directly supported by some Qt versions conveniently,
        // but let's assume standard multiplication works or do manually.
        // Actually QVector3D operator*(QMatrix4x4) exists but is row-vector mul.
        // QMatrix4x4 operator*(QVector3D) is standard column-vector mul.
        
        // Use generic map method or just manual multiply if needed.
        // Easier: mapVector for vectors (ignores translation) but needs to be normal matrix for non-uniform scales.
        // If scale is uniform, mapVector is fine.
        
        // Let's do it manually to be safe with QMatrix3x3
        const float *d = normalMat.constData();
        float nx = d[0]*v.norm.x() + d[3]*v.norm.y() + d[6]*v.norm.z();
        float ny = d[1]*v.norm.x() + d[4]*v.norm.y() + d[7]*v.norm.z();
        float nz = d[2]*v.norm.x() + d[5]*v.norm.y() + d[8]*v.norm.z();
        v.norm = QVector3D(nx, ny, nz).normalized();
    }
    markVertexDirty();
    m_bAABBDirty = true;
}
 
//=============================================================================================================
 
void BrainSurface::applyTransform(const QMatrix4x4 &m)
{
    m_vertexData = m_originalVertexData;
    if (!m.isIdentity()) {
        transform(m);
    } else {
        markVertexDirty();
        m_bAABBDirty = true;
    }
}
 
//=============================================================================================================
 
void BrainSurface::updateBuffers(QRhi *rhi, QRhiResourceUpdateBatch *u)
{
    const bool needsCreate = !m_gpu->vertexBuffer || !m_gpu->indexBuffer;
 
#ifdef __EMSCRIPTEN__
    // On WebGL/WASM, always re-upload vertex and index data.
    // The QRhi GLES2 backend's VAO cache can lose its element-buffer
    // binding between frames, causing draws to produce no output.
    // Re-uploading via uploadStaticBuffer triggers the necessary
    // glBindBuffer calls that refresh the GL state.
    if (!needsCreate && !m_gpu->dirty) {
        // Buffers exist and data hasn't changed — still re-upload
        u->uploadStaticBuffer(m_gpu->vertexBuffer.get(), m_vertexData.constData());
        u->uploadStaticBuffer(m_gpu->indexBuffer.get(), m_indexData.constData());
        return;
    }
#else
    if (!m_gpu->dirty && !needsCreate) return;
#endif
 
    const quint32 vbufSize = m_vertexData.size() * sizeof(VertexData);
    const quint32 ibufSize = m_indexData.size() * sizeof(uint32_t);
 
    if (!m_gpu->vertexBuffer) {
        // Use Dynamic VBO on desktop: allows in-place updateDynamicBuffer()
        // for STC animation colour updates without buffer destroy/recreate.
        // WASM uses merged per-category buffers, so this path is only
        // relevant for desktop (where it's Immutable as a fallback).
#ifdef __EMSCRIPTEN__
        m_gpu->vertexBuffer.reset(rhi->newBuffer(QRhiBuffer::Immutable, QRhiBuffer::VertexBuffer, vbufSize));
#else
        m_gpu->vertexBuffer.reset(rhi->newBuffer(QRhiBuffer::Dynamic, QRhiBuffer::VertexBuffer, vbufSize));
#endif
        m_gpu->vertexBuffer->create();
        m_gpu->indexDirty = true; // new VBO implies IBO also needs upload
    }
    if (!m_gpu->indexBuffer) {
        m_gpu->indexBuffer.reset(rhi->newBuffer(QRhiBuffer::Immutable, QRhiBuffer::IndexBuffer, ibufSize));
        m_gpu->indexBuffer->create();
        m_gpu->indexDirty = true;
    }
 
#ifdef __EMSCRIPTEN__
    u->uploadStaticBuffer(m_gpu->vertexBuffer.get(), m_vertexData.constData());
    u->uploadStaticBuffer(m_gpu->indexBuffer.get(), m_indexData.constData());
#else
    // Desktop: Dynamic VBO — in-place update, no buffer recreation.
    // Only re-upload IBO when topology actually changed (initial load,
    // geometry swap).  STC colour animation only touches the VBO.
    u->updateDynamicBuffer(m_gpu->vertexBuffer.get(), 0, vbufSize, m_vertexData.constData());
    if (m_gpu->indexDirty) {
        u->uploadStaticBuffer(m_gpu->indexBuffer.get(), m_indexData.constData());
        m_gpu->indexDirty = false;
    }
#endif
    m_gpu->dirty = false;
}
 
//=============================================================================================================
 
std::vector<Eigen::VectorXi> BrainSurface::computeNeighbors() const
{
    // Use temporary std::vector<std::set<int>> for dedup during construction
    std::vector<std::set<int>> tempNeighbors(m_vertexData.size());
    
    // Triangles are stored in m_indexData as triplets
    for (int i = 0; i + 2 < m_indexData.size(); i += 3) {
        int v0 = m_indexData[i];
        int v1 = m_indexData[i + 1];
        int v2 = m_indexData[i + 2];
        
        // Add bidirectional edges (set handles dedup)
        tempNeighbors[v0].insert(v1);
        tempNeighbors[v0].insert(v2);
        tempNeighbors[v1].insert(v0);
        tempNeighbors[v1].insert(v2);
        tempNeighbors[v2].insert(v0);
        tempNeighbors[v2].insert(v1);
    }
    
    // Convert to std::vector<VectorXi>
    std::vector<Eigen::VectorXi> neighbors(tempNeighbors.size());
    for (size_t k = 0; k < tempNeighbors.size(); ++k) {
        const auto& s = tempNeighbors[k];
        neighbors[k].resize(static_cast<Eigen::Index>(s.size()));
        Eigen::Index idx = 0;
        for (int val : s) {
            neighbors[k][idx++] = val;
        }
    }
    
    return neighbors;
}
 
//=============================================================================================================
 
Eigen::MatrixX3f BrainSurface::verticesAsMatrix() const
{
    Eigen::MatrixX3f mat(m_vertexData.size(), 3);
    for (int i = 0; i < m_vertexData.size(); ++i) {
        mat(i, 0) = m_vertexData[i].pos.x();
        mat(i, 1) = m_vertexData[i].pos.y();
        mat(i, 2) = m_vertexData[i].pos.z();
    }
    return mat;
}
 
 
//=============================================================================================================
 
void BrainSurface::boundingBox(QVector3D &min, QVector3D &max) const
{
    if (!m_bAABBDirty) {
        min = m_aabbMin;
        max = m_aabbMax;
        return;
    }
 
    if (m_vertexData.isEmpty()) {
        min = QVector3D(0,0,0);
        max = QVector3D(0,0,0);
        return;
    }
 
    min = m_vertexData[0].pos;
    max = m_vertexData[0].pos;
 
    for (const auto &v : m_vertexData) {
        min.setX(std::min(min.x(), v.pos.x()));
        min.setY(std::min(min.y(), v.pos.y()));
        min.setZ(std::min(min.z(), v.pos.z()));
        
        max.setX(std::max(max.x(), v.pos.x()));
        max.setY(std::max(max.y(), v.pos.y()));
        max.setZ(std::max(max.z(), v.pos.z()));
    }
 
    m_aabbMin = min;
    m_aabbMax = max;
    m_bAABBDirty = false;
}
 
bool BrainSurface::intersects(const QVector3D &rayOrigin, const QVector3D &rayDir, float &dist, int &vertexIdx) const
{
    vertexIdx = -1;
    if (m_vertexData.isEmpty()) return false;
 
    // 1. AABB Check (Cached)
    QVector3D min, max;
    boundingBox(min, max);
    
    // Ray-AABB slab method (Double Precision for stability)
    double eps = 1e-4;
    double origin[3] = {rayOrigin.x(), rayOrigin.y(), rayOrigin.z()};
    double dir[3] = {rayDir.x(), rayDir.y(), rayDir.z()};
    double minB[3] = {min.x() - eps, min.y() - eps, min.z() - eps};
    double maxB[3] = {max.x() + eps, max.y() + eps, max.z() + eps};
 
    // Extract individual components for triangle intersection (Möller–Trumbore)
    double originX = origin[0], originY = origin[1], originZ = origin[2];
    double dirX = dir[0], dirY = dir[1], dirZ = dir[2];
 
    double tmin = -std::numeric_limits<double>::max();
    double tmax = std::numeric_limits<double>::max();
 
    for (int i = 0; i < 3; ++i) {
        if (std::abs(dir[i]) < 1e-15) {
            if (origin[i] < minB[i] || origin[i] > maxB[i]) return false;
        } else {
            double t1 = (minB[i] - origin[i]) / dir[i];
            double t2 = (maxB[i] - origin[i]) / dir[i];
            if (t1 > t2) std::swap(t1, t2);
            if (t1 > tmin) tmin = t1;
            if (t2 < tmax) tmax = t2;
            if (tmin > tmax) return false;
        }
    }
    
    if (tmax < 1e-7) return false;
 
    // 2. Triangle intersection
    double closestDist = std::numeric_limits<double>::max();
    bool hit = false;
    int closestVert = -1;
    
    // Brute-force triangle intersection using Double Precision Möller–Trumbore
    // Note: For 100k+ vertices this is slow, ideally we'd use an Octree/BVH.
    // However, since we only do this on mouse-over, results are usually acceptable if not too many surfaces are active.
    for (int i = 0; i < m_indexData.size(); i += 3) {
        int i0 = m_indexData[i];
        int i1 = m_indexData[i+1];
        int i2 = m_indexData[i+2];
        const QVector3D &v0q = m_vertexData[i0].pos;
        const QVector3D &v1q = m_vertexData[i1].pos;
        const QVector3D &v2q = m_vertexData[i2].pos;
 
        double v0x = v0q.x(), v0y = v0q.y(), v0z = v0q.z();
        double v1x = v1q.x(), v1y = v1q.y(), v1z = v1q.z();
        double v2x = v2q.x(), v2y = v2q.y(), v2z = v2q.z();
        
        double edge1x = v1x - v0x, edge1y = v1y - v0y, edge1z = v1z - v0z;
        double edge2x = v2x - v0x, edge2y = v2y - v0y, edge2z = v2z - v0z;
        
        double hx = dirY * edge2z - dirZ * edge2y;
        double hy = dirZ * edge2x - dirX * edge2z;
        double hz = dirX * edge2y - dirY * edge2x;
        
        double a = edge1x * hx + edge1y * hy + edge1z * hz;
        if (std::abs(a) < 1e-18) continue; // Purely parallel
        
        double f = 1.0 / a;
        double sx = originX - v0x, sy = originY - v0y, sz = originZ - v0z;
        double u = f * (sx * hx + sy * hy + sz * hz);
        if (u < -1e-7 || u > 1.0000001) continue;
        
        double qx = sy * edge1z - sz * edge1y;
        double qy = sz * edge1x - sx * edge1z;
        double qz = sx * edge1y - sy * edge1x;
        
        double v = f * (dirX * qx + dirY * qy + dirZ * qz);
        if (v < -1e-7 || u + v > 1.0000001) continue;
        
        double t = f * (edge2x * qx + edge2y * qy + edge2z * qz);
        if (t > 1e-7 && t < closestDist) {
            // Check barycentric coordinates with Fixed Relative Tolerance (25%)
            // Relative tolerance scales with triangle size:
            // - Large triangles (Helmet): Tolerates large gaps (~cm scale)
            // - Small triangles (Brain): Tolerates small errors (~mm scale)
            // This prevents clicking 'through' sparse meshes while staying precise on dense ones.
            constexpr double tol = 0.25; 
            
            if (u >= -tol && v >= -tol && u + v <= 1.0 + tol) {
                closestDist = t;
                hit = true;
                
                // Find closest vertex of the hit triangle to the hit point
                // (Used for region lookup)
                double hitX = originX + t * dirX;
                double hitY = originY + t * dirY;
                double hitZ = originZ + t * dirZ;
                
                double d0 = (v0x - hitX)*(v0x - hitX) + (v0y - hitY)*(v0y - hitY) + (v0z - hitZ)*(v0z - hitZ);
                double d1 = (v1x - hitX)*(v1x - hitX) + (v1y - hitY)*(v1y - hitY) + (v1z - hitZ)*(v1z - hitZ);
                double d2 = (v2x - hitX)*(v2x - hitX) + (v2y - hitY)*(v2y - hitY) + (v2z - hitZ)*(v2z - hitZ);
                
                if (d0 < d1 && d0 < d2) closestVert = i0;
                else if (d1 < d2) closestVert = i1;
                else closestVert = i2;
            }
        }
    }
    
    if (hit) {
        dist = static_cast<float>(closestDist);
        vertexIdx = closestVert;
        return true;
    }
    
    return false;
}
 
//=============================================================================================================
 
QString BrainSurface::getAnnotationLabel(int vertexIdx) const
{
    if (!m_hasAnnotation || vertexIdx < 0 || vertexIdx >= m_vertexData.size()) {
        return "";
    }
 
    const Eigen::VectorXi &vertices = m_annotation.getVertices();
    const Eigen::VectorXi &labelIds = m_annotation.getLabelIds();
    const FSLIB::FsColortable &ct = m_annotation.getColortable();
 
    // The .annot file might not contain all vertices if it's sparse, 
    // but usually it contains a mapping for all.
    // Let's find the labelId for this vertex.
    int labelId = -1;
    for (int i = 0; i < vertices.rows(); ++i) {
        if (vertices(i) == vertexIdx) {
            labelId = labelIds(i);
            break;
        }
    }
 
    if (labelId == -1) return "Unknown";
 
    // Find the name in colortable
    for (int i = 0; i < ct.numEntries; ++i) {
        if (ct.table(i, 4) == labelId) {
            QString name = ct.struct_names[i];
            // Remove null characters and trailing whitespace that might cause "strange signs"
            while (!name.isEmpty() && (name.endsWith('\0') || name.endsWith(' '))) {
                name.chop(1);
            }
            return name;
        }
    }
 
    return "Unknown";
}
 
int BrainSurface::getAnnotationLabelId(int vertexIdx) const
{
    if (!m_hasAnnotation || vertexIdx < 0) return -1;
 
    const Eigen::VectorXi &vertices = m_annotation.getVertices();
    const Eigen::VectorXi &labelIds = m_annotation.getLabelIds();
 
    for (int i = 0; i < vertices.rows(); ++i) {
        if (vertices(i) == vertexIdx) {
            return labelIds(i);
        }
    }
 
    return -1;
}
 
void BrainSurface::setSelectedRegion(int regionId)
{
    m_selectedRegionId = regionId;
    updateVertexColors();
}
 
void BrainSurface::setSelected(bool selected)
{
    m_selected = selected;
    updateVertexColors();
}
 
void BrainSurface::setSelectedVertexRange(int start, int count)
{
    m_selectedVertexStart = start;
    m_selectedVertexCount = count;
    updateVertexColors();
}