brainsurface.cpp

Go to the documentation of this file.

//=============================================================================================================
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "brainsurface.h"
 
#include <rhi/qrhi.h>
 
#include <set>
 
namespace DISP3DLIB
{
 
namespace
{
uint32_t withAlpha(uint32_t color, uint32_t alpha)
{
    return (color & 0x00FFFFFFu) | ((alpha & 0xFFu) << 24);
}
 
uint32_t curvatureGray(const QVector<float>& curvature, int index)
{
    return (index >= 0 && index < curvature.size() && curvature[index] > 0.0f)
        ? 0x40u
        : 0xAAu;
}
 
// Overlay alpha is blended into RGB here because the shader keeps the curvature grey in the alpha channel.
uint32_t overCortex(uint32_t overlay, uint32_t gray)
{
    const uint32_t a = overlay >> 24;
    uint32_t rgb = 0;
    for (int shift = 0; shift < 24; shift += 8) {
        const uint32_t c = (overlay >> shift) & 0xFFu;
        rgb |= ((c * a + gray * (255u - a) + 127u) / 255u) << shift;
    }
    return withAlpha(rgb, gray);
}
} // namespace
 
//=============================================================================================================
// 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)
{
    createFromData(vertices, FSLIB::FsSurface::compute_normals(vertices, triangles), 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. Preserve neutral
    // curvature grey in alpha so Surface mode can still render the classic
    // light/dark cortex while RGB is occupied by source-estimate colours.
    for (int i = 0; i < qMin(colors.size(), m_vertexData.size()); ++i) {
        m_vertexData[i].color = overCortex(colors[i], curvatureGray(m_curvature, 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::setColor(const QColor& color)
{
    m_defaultColor = color;
    m_baseColor = color;
    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 = curvatureGray(m_curvature, i);
            m_vertexData[i].color = packABGR(val, 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 = overCortex(m_stcColors[i], curvatureGray(m_curvature, 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)
{
    // Only the geometry is reset; colours and overlays set since loading stay
    for (int i = 0; i < qMin(m_vertexData.size(), m_originalVertexData.size()); ++i) {
        m_vertexData[i].pos = m_originalVertexData[i].pos;
        m_vertexData[i].norm = m_originalVertexData[i].norm;
    }
    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
    // Desktop: VBO is Immutable.  Nothing to do unless the data has
    // actually changed (markVertexDirty()) or the buffer hasn't been
    // created yet.  When dirty we recreate the VBO below to swap in the
    // new vertex data — this is intentional: an Immutable buffer cannot
    // be partially updated, and STC colour animation modifies the colour
    // channel of every vertex anyway.
    if (!m_gpu->dirty && !needsCreate)
        return;
#endif
 
    const quint32 vbufSize = static_cast<quint32>(m_vertexData.size() * sizeof(VertexData));
    const quint32 ibufSize = static_cast<quint32>(m_indexData.size() * sizeof(uint32_t));
 
#ifdef __EMSCRIPTEN__
    if (!m_gpu->vertexBuffer) {
        m_gpu->vertexBuffer.reset(rhi->newBuffer(QRhiBuffer::Immutable, QRhiBuffer::VertexBuffer, vbufSize));
        m_gpu->vertexBuffer->create();
        m_gpu->indexDirty = true; // new VBO implies IBO also needs upload
    }
#else
    // Desktop: recreate the Immutable VBO whenever data is dirty so the
    // new vertex contents take effect.  Using Immutable (re-created on
    // change) instead of QRhiBuffer::Dynamic avoids two issues with
    // Dynamic buffers on Metal/Vulkan for multi-MB vertex data:
    //   1) Dynamic is documented as intended for small, frequently-updated
    //      payloads (UBOs); large Dynamic buffers exhibit unstable
    //      behaviour across in-flight slot rotation, including stale-slot
    //      reads (frozen STC frames) and out-of-range geometry artefacts
    //      ("stretched triangle" glitches).
    //   2) Each in-flight frame slot allocates its own physical buffer,
    //      multiplying memory pressure ~3x for every visible surface.
    if (!m_gpu->vertexBuffer || m_gpu->dirty) {
        m_gpu->vertexBuffer.reset(rhi->newBuffer(QRhiBuffer::Immutable, QRhiBuffer::VertexBuffer, vbufSize));
        m_gpu->vertexBuffer->create();
        m_gpu->indexDirty = true; // new VBO implies IBO also needs upload
    }
#endif
    if (!m_gpu->indexBuffer) {
        m_gpu->indexBuffer.reset(rhi->newBuffer(QRhiBuffer::Immutable, QRhiBuffer::IndexBuffer, ibufSize));
        m_gpu->indexBuffer->create();
        m_gpu->indexDirty = true;
    }
 
    u->uploadStaticBuffer(m_gpu->vertexBuffer.get(), m_vertexData.constData());
    if (m_gpu->indexDirty) {
        u->uploadStaticBuffer(m_gpu->indexBuffer.get(), m_indexData.constData());
        m_gpu->indexDirty = false;
    }
    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();
}
 
} // namespace DISP3DLIB