dipoleobject.cpp

Go to the documentation of this file.

//=============================================================================================================
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "dipoleobject.h"
 
#include <cmath>
#include <QQuaternion>
#include <QRandomGenerator>
 
namespace DISP3DLIB
{
 
//=============================================================================================================
// PIMPL
//=============================================================================================================
 
DipoleObject::DipoleObject() = default;
 
DipoleObject::~DipoleObject() = default;
 
InstancedStreams DipoleObject::streams() const
{
    return {{m_vertexData.constData(), quint32(m_vertexData.size()), m_geometryGeneration},
            {m_indexData.constData(), quint32(m_indexData.size()), m_geometryGeneration},
            {m_instanceData.constData(), quint32(m_instanceCount * sizeof(InstanceData)), m_instanceGeneration},
            m_indexCount,
            m_instanceCount};
}
 
void DipoleObject::load(const INVLIB::InvEcdSet& ecdSet)
{
    createGeometry();
 
    m_instanceCount = ecdSet.size();
    m_instanceData.resize(m_instanceCount * sizeof(InstanceData));
    m_loadedModels.resize(m_instanceCount);
    InstanceData* data = reinterpret_cast<InstanceData*>(m_instanceData.data());
 
    QVector3D from(0.0f, 1.0f, 0.0f); // Cone points up Y axis
 
    float maxMag = 0.0f;
    for (int i = 0; i < ecdSet.size(); ++i) {
        float mag = std::sqrt(std::pow(ecdSet[i].Q(0), 2) + std::pow(ecdSet[i].Q(1), 2) + std::pow(ecdSet[i].Q(2), 2));
        maxMag = std::max(maxMag, mag);
    }
 
    for (int i = 0; i < ecdSet.size(); ++i) {
        const auto& dip = ecdSet[i];
 
        QVector3D pos(dip.rd(0), dip.rd(1), dip.rd(2));
        QVector3D Q(dip.Q(0), dip.Q(1), dip.Q(2));
        float mag = Q.length();
 
        QVector3D to = Q.normalized();
        // Handle case where 'to' is parallel to 'from'
        QQuaternion rot;
        if (QVector3D::dotProduct(from, to) > 0.99f) {
            rot = QQuaternion();
        } else if (QVector3D::dotProduct(from, to) < -0.99f) {
            rot = QQuaternion::fromAxisAndAngle(1.0f, 0.0f, 0.0f, 180.0f);
        } else {
            rot = QQuaternion::rotationTo(from, to);
        }
 
        // Scale based on magnitude relative to max, with a minimum size of 20%
        float scaleFactor = (maxMag > 0.0f) ? (0.2f + 0.8f * (mag / maxMag)) : 1.0f;
 
        QMatrix4x4 m;
        m.translate(pos);
        m.rotate(rot);
        m.scale(scaleFactor);
        m_loadedModels[i] = m;
 
        const float* mPtr = m.constData();
        for (int j = 0; j < 16; ++j) {
            data[i].model[j] = mPtr[j];
        }
 
        // Random Color
        data[i].color[0] = QRandomGenerator::global()->generateDouble();
        data[i].color[1] = QRandomGenerator::global()->generateDouble();
        data[i].color[2] = QRandomGenerator::global()->generateDouble();
        data[i].color[3] = 1.0f;
 
        // Selection state
        data[i].isSelected = 0.0f;
    }
 
    ++m_instanceGeneration;
}
 
void DipoleObject::applyTransform(const QMatrix4x4& trans)
{
    if (m_instanceCount == 0)
        return;
 
    InstanceData* data = reinterpret_cast<InstanceData*>(m_instanceData.data());
 
    for (int i = 0; i < m_instanceCount; ++i) {
        const QMatrix4x4 newModel = trans * m_loadedModels[i];
        const float* newPtr = newModel.constData();
        for (int j = 0; j < 16; ++j) {
            data[i].model[j] = newPtr[j];
        }
    }
 
    ++m_instanceGeneration;
}
 
bool DipoleObject::boundingBox(QVector3D& min, QVector3D& max) const
{
    if (m_instanceCount == 0)
        return false;
    const InstanceData* data = reinterpret_cast<const InstanceData*>(m_instanceData.constData());
    for (int i = 0; i < m_instanceCount; ++i) {
        // Column-major model matrix: the translation is elements 12..14
        const QVector3D pos(data[i].model[12], data[i].model[13], data[i].model[14]);
        min = i == 0 ? pos : QVector3D(std::min(min.x(), pos.x()), std::min(min.y(), pos.y()), std::min(min.z(), pos.z()));
        max = i == 0 ? pos : QVector3D(std::max(max.x(), pos.x()), std::max(max.y(), pos.y()), std::max(max.z(), pos.z()));
    }
    return true;
}
 
void DipoleObject::createGeometry()
{
    if (!m_vertexData.isEmpty())
        return;
 
    // Create a simple cone
    // Radius 0.001, Height 0.003
    // 32 segments
 
    float radius = 0.005f; // Slightly larger for visibility
    float height = 0.01f;
    int segments = 16;
 
    std::vector<VertexData> vertices;
    std::vector<uint32_t> indices;
 
    // Tip
    VertexData tip = {0, height, 0, 0, 1, 0};
    vertices.push_back(tip);
 
    // Base center
    VertexData baseCenter = {0, 0, 0, 0, -1, 0};
    vertices.push_back(baseCenter); // Index 1
 
    // Rim vertices
    for (int i = 0; i < segments; ++i) {
        float angle = 2.0f * M_PI * i / segments;
        float x = radius * cos(angle);
        float z = radius * sin(angle);
 
        // Side normal
        // Slope vector: (cos, -h/r, sin) -> normalized
        QVector3D n(x, radius / height, z);
        n.normalize();
 
        VertexData vSide = {x, 0, z, n.x(), n.y(), n.z()};
        vertices.push_back(vSide);
 
        VertexData vBase = {x, 0, z, 0, -1, 0};
        vertices.push_back(vBase);
    }
 
    // Indices
    for (int i = 0; i < segments; ++i) {
        // Cone sides
        // Tip (0), current(i), next(next)
        // Side vertices start at firstRimIdx.
        // Layout: Tip, BaseCenter, Side0, Base0, Side1, Base1...
        // Wait, easier to keep separate lists or just flat
 
        // Re-do layout for ease:
        // 0: Tip
        // 1: Base Center
        // 2..2+segments-1: Side rim vertices
        // 2+segments..2+2*segments-1: Base rim vertices
    }
 
    vertices.clear();
    vertices.push_back(tip);        // 0
    vertices.push_back(baseCenter); // 1
 
    for (int i = 0; i < segments; ++i) {
        float angle = 2.0f * M_PI * i / segments;
        float x = radius * cos(angle);
        float z = radius * sin(angle);
 
        QVector3D n(x, 0, z); // Approximate normal for side (flat shading effectively if we don't slant)
        // Better normal: vector perpendicular to slope
        QVector3D slope(x, -height, z);
        QVector3D tangent(-sin(angle), 0, cos(angle));
        QVector3D sideNormal = QVector3D::crossProduct(tangent, slope).normalized(); // Wait, slope is vector down side.
        // Actually simple: normal y component is radius/height ratio related.
 
        vertices.push_back({x, 0, z, sideNormal.x(), sideNormal.y(), sideNormal.z()}); // Side rim
    }
 
    for (int i = 0; i < segments; ++i) {
        float angle = 2.0f * M_PI * i / segments;
        float x = radius * cos(angle);
        float z = radius * sin(angle);
        vertices.push_back({x, 0, z, 0, -1, 0}); // Base rim
    }
 
    int sideStart = 2;
    int baseStart = 2 + segments;
 
    for (int i = 0; i < segments; ++i) {
        int next = (i + 1) % segments;
 
        // Cone face
        indices.push_back(0);
        indices.push_back(sideStart + next);
        indices.push_back(sideStart + i);
 
        // Base face
        indices.push_back(1);
        indices.push_back(baseStart + i);
        indices.push_back(baseStart + next);
    }
 
    m_indexCount = static_cast<int>(indices.size());
 
    m_vertexData.resize(vertices.size() * sizeof(VertexData));
    memcpy(m_vertexData.data(), vertices.data(), m_vertexData.size());
 
    m_indexData.resize(indices.size() * sizeof(uint32_t));
    memcpy(m_indexData.data(), indices.data(), m_indexData.size());
 
    ++m_geometryGeneration;
}
 
int DipoleObject::intersect(const QVector3D& rayOrigin, const QVector3D& rayDir, float& dist) const
{
    if (m_instanceCount == 0)
        return -1;
 
    int closestIdx = -1;
    float closestDist = std::numeric_limits<float>::max();
 
    const InstanceData* data = reinterpret_cast<const InstanceData*>(m_instanceData.constData());
 
    // Geometry radius ~ 0.005 (base) to 0.01 (height).
    // Let's use a slightly larger hit radius to make selection easier and more accurate.
    const float baseRadius = 0.02f;
 
    for (int i = 0; i < m_instanceCount; ++i) {
        // Extract translation (last column)
        QVector3D center(data[i].model[12], data[i].model[13], data[i].model[14]);
 
        // Extract scale (length of first column) - assumes uniform scale roughly
        QVector3D col0(data[i].model[0], data[i].model[1], data[i].model[2]);
        float scale = col0.length();
 
        float radius = baseRadius * scale;
 
        // Ray-Sphere Intersection
        QVector3D L = center - rayOrigin;
        float tca = QVector3D::dotProduct(L, rayDir);
 
        if (tca < 0)
            continue; // Behind ray
 
        float d2 = QVector3D::dotProduct(L, L) - tca * tca;
        if (d2 > radius * radius)
            continue; // Miss
 
        float thc = std::sqrt(radius * radius - d2);
        float t0 = tca - thc;
        float t1 = tca + thc;
 
        float t = t0;
        if (t < 0)
            t = t1;
        if (t < 0)
            continue;
 
        if (t < closestDist) {
            closestDist = t;
            closestIdx = i;
        }
    }
 
    if (closestIdx != -1) {
        dist = closestDist;
        return closestIdx;
    }
 
    return -1;
}
 
void DipoleObject::setSelected(int index, bool selected)
{
    if (index < 0 || index >= m_instanceCount)
        return;
 
    InstanceData* data = reinterpret_cast<InstanceData*>(m_instanceData.data());
 
    data[index].isSelected = selected ? 1.0f : 0.0f;
 
    ++m_instanceGeneration;
}
 
} // namespace DISP3DLIB