Added KDTree stars to black hole render

Co-Authored-By: Emil Wallberg <49481622+EmilWallberg@users.noreply.github.com>

modules/blackhole/CMakeLists.txt

@@ -29,6 +29,7 @@ set(HEADER_FILES
   ${CMAKE_CURRENT_SOURCE_DIR}/blackholemodule.h
   ${CMAKE_CURRENT_SOURCE_DIR}/rendering/renderableblackhole.h
   ${CMAKE_CURRENT_SOURCE_DIR}/rendering/viewport.h
+  ${CMAKE_CURRENT_SOURCE_DIR}/rendering/kdtree.h
 )
 
 # Define Source Files
@@ -36,6 +37,7 @@ set(SOURCE_FILES
   ${CMAKE_CURRENT_SOURCE_DIR}/blackholemodule.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/rendering/renderableblackhole.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/rendering/viewport.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/rendering/kdtree.cpp
 )
 
 # Define Shader Files
@@ -63,8 +65,8 @@ project(blackhole_cuda_prj CUDA)
 find_package(CUDA REQUIRED)
 # Define CUDA Files
 set(CUDA_FILES
-  ${CMAKE_CURRENT_SOURCE_DIR}/cuda/blackhole_cuda.cu
   ${CMAKE_CURRENT_SOURCE_DIR}/cuda/blackhole_cuda.h
+  ${CMAKE_CURRENT_SOURCE_DIR}/cuda/blackhole_cuda.cu
 )
 source_group("CUDA Files" FILES ${CUDA_FILES})
 

modules/blackhole/rendering/kdtree.cpp

@@ -0,0 +1,96 @@
+#include "kdtree.h"
+#include <ghoul/filesystem/filesystem.h>
+#include <openspace/data/dataloader.h>
+#include <algorithm>
+
+#include <queue>
+
+namespace {
+    glm::vec3 cartesianToSpherical(const glm::vec3& cartesian) {
+        float radius = glm::length(cartesian);
+        float theta = std::atan2(glm::sqrt(cartesian.x * cartesian.x + cartesian.y * cartesian.y), cartesian.z);
+        float phi = std::atan2(cartesian.y, cartesian.x);
+
+        return glm::vec3(radius, theta, phi);
+    }
+}
+namespace openspace {
+    void KDTree::build(const std::string& filePath, const glm::vec3& localWorldCenter) {
+        const std::filesystem::path file{ absPath(filePath) };
+        dataloader::Dataset dataset = dataloader::data::loadFileWithCache(file);
+        size_t numberOfStars =  dataset.entries.size();
+        tree.resize(numberOfStars);
+        struct NodeInfo {
+            size_t index;
+            size_t depth;
+            size_t start;
+            size_t end;
+        };
+
+        #pragma omp parallel for
+        for (auto& entry : dataset.entries) {
+            entry.position = cartesianToSpherical(entry.position - localWorldCenter);
+        }
+
+        std::queue<NodeInfo> q;
+        q.emplace(0, 0, 0, numberOfStars);
+        while (!q.empty()) {
+            NodeInfo node{ q.front() };
+            q.pop();
+
+            if (node.start >= node.end) continue;
+
+            int axis = node.depth % 2 + 1;
+
+            auto comparetor = [axis](dataloader::Dataset::Entry const& a, dataloader::Dataset::Entry const& b) -> bool {
+                return a.position[axis] < b.position[axis];
+                };
+
+            // Sort to find median
+            std::sort(
+                dataset.entries.begin() + node.start,
+                dataset.entries.begin() + node.end,
+                comparetor
+            );
+
+            size_t medianIndex{ (node.start + node.end) / 2 };
+            dataloader::Dataset::Entry const& entry{ dataset.entries[medianIndex] };
+
+            glm::vec3 const& position{ entry.position };
+            float const color{ entry.data[0]};
+            float const lum{ entry.data[1]};
+            float const absMag{ entry.data[2]};
+
+            if (node.index >= tree.size()) {
+                tree.resize(std::max(tree.size() * 2, node.index + 1));
+            }
+
+            tree.emplace(
+                tree.begin() + node.index, position, color, lum, absMag
+            );
+
+            // Enqueue left and right children
+            q.emplace(2 * node.index + 1, node.depth + 1, node.start, medianIndex);
+            q.emplace(2 * node.index + 2, node.depth + 1, medianIndex + 1, node.end);
+        }
+    }
+
+    std::vector<float> KDTree::flatTree() const {
+        std::vector<float> flatData;
+        flatData.resize(tree.size() * 6);
+
+        for (int i = 0; i < tree.size(); i++) {
+            Node const& node{ tree[i] };
+
+            size_t index = i * 6;
+            flatData[index] = node.position.x;
+            flatData[index + 1] = node.position.y;
+            flatData[index + 2] = node.position.z;
+            flatData[index + 3] = node.color;
+            flatData[index + 4] = node.lum;
+            flatData[index + 5] = node.absMag;
+        }
+
+        return flatData;
+    }
+}

modules/blackhole/rendering/kdtree.h

@@ -0,0 +1,31 @@
+#ifndef __OPENSPACE_MODULE_BLACKHOLE___KDTREE___H__
+#define __OPENSPACE_MODULE_BLACKHOLE___KDTREE___H__
+#include <glm/glm.hpp>
+#include <vector>
+#include <string>
+
+namespace openspace {
+    namespace kdtree{}
+    class KDTree {
+    public:
+        KDTree() {};
+
+        size_t size() { return tree.size() * 6; };
+
+        void build(const std::string& filePath, const glm::vec3& localWorldCenter);
+
+        std::vector<float> flatTree() const;
+
+    private:
+        struct Node {
+            glm::fvec3 position;
+            float color;
+            float lum;
+            float absMag;
+        };
+
+        std::vector<Node> tree{};
+    };
+}
+
+#endif

modules/blackhole/rendering/renderableblackhole.cpp

@@ -66,8 +66,11 @@ namespace openspace {
     RenderableBlackHole::~RenderableBlackHole() {}
 
     void RenderableBlackHole::initialize() {
-        global::navigationHandler->camera()->setRotation(glm::dquat(0,0,0,0));
         _schwarzschildWarpTable = std::vector<float>(_rayCount * 2, 0.f);
+
+        _starKDTree.build("${BASE}/sync/http/stars_du/6/stars.speck", glm::vec3(0));
+
+        flatDataStar = _starKDTree.flatTree();
     }
 
     void RenderableBlackHole::initializeGL() {
@@ -99,6 +102,7 @@ namespace openspace {
     bool RenderableBlackHole::isReady() const {
         return _program != nullptr;
     }
+
     void RenderableBlackHole::update(const UpdateData&) {
         glm::vec3 cameraPosition = global::navigationHandler->camera()->positionVec3();
         glm::vec3 anchorNodePosition = global::navigationHandler->anchorNode()->position();
@@ -110,7 +114,8 @@ namespace openspace {
 
             schwarzchild(_rs, _rEnvmap, _rayCount, _stepsCount, 1.0f / _rCamera, _stepLength, _schwarzschildWarpTable.data());
         }
-        bindSSBOData(_program, "ssbo_warp_table", _ssboDataBinding, _ssboData);
+        bindSSBOData(_program, "ssbo_warp_table", _ssboSchwarzschildDataBinding, _ssboSchwarzschildWarpTable);
+        bindSSBOData(_program, "ssbo_star_map", _ssboStarDataBinding, _ssboStarKDTree);
     }
 
     void RenderableBlackHole::render(const RenderData& renderData, RendererTasks&) {
@@ -153,8 +158,7 @@ namespace openspace {
 
     void RenderableBlackHole::SendSchwarzchildTableToShader()
     {
-        // Update SSBO Index array with accumulated stars in all chunks.
-        glBindBuffer(GL_SHADER_STORAGE_BUFFER, _ssboData);
+        glBindBuffer(GL_SHADER_STORAGE_BUFFER, _ssboSchwarzschildWarpTable);
 
         const size_t indexBufferSize = _schwarzschildWarpTable.size() * sizeof(float);
 
@@ -168,6 +172,22 @@ namespace openspace {
         glBindBuffer(GL_SHADER_STORAGE_BUFFER, 0);
     }
 
+    void RenderableBlackHole::SendStarKDTreeToShader()
+    {
+        glBindBuffer(GL_SHADER_STORAGE_BUFFER, _ssboStarKDTree);
+
+        const size_t indexBufferSize = flatDataStar.size() * sizeof(float);
+
+        glBufferData(
+            GL_SHADER_STORAGE_BUFFER,
+            indexBufferSize,
+            flatDataStar.data(),
+            GL_STREAM_DRAW
+        );
+
+        glBindBuffer(GL_SHADER_STORAGE_BUFFER, 0);
+    }
+
     void RenderableBlackHole::setupShaders() {
         const std::string vertexShaderPath = "${MODULE_BLACKHOLE}/shaders/blackhole_vs.glsl";
         const std::string fragmentShaderPath = "${MODULE_BLACKHOLE}/shaders/blackhole_fs.glsl";

modules/blackhole/rendering/renderableblackhole.h

@@ -7,6 +7,7 @@
 #include <ghoul/opengl/textureunit.h>
 #include <ghoul/opengl/bufferbinding.h>
 #include <openspace/properties/scalar/floatproperty.h>
+#include <modules/blackhole/rendering/kdtree.h>
 
 namespace openspace {
 
@@ -23,13 +24,14 @@ namespace openspace {
         bool isReady() const override;
 
         void render(const RenderData& data, RendererTasks& rendererTask) override;
-        void SendSchwarzchildTableToShader();
         void update(const UpdateData& data) override;
 
 
         static documentation::Documentation Documentation();
 
     private:
+        void SendSchwarzchildTableToShader();
+        void SendStarKDTreeToShader();
         void bindSSBOData(ghoul::opengl::ProgramObject* program,
             const std::string& ssboName,
             std::unique_ptr<ghoul::opengl::BufferBinding<ghoul::opengl::bufferbinding::Buffer::ShaderStorage>>& ssboBinding,
@@ -54,13 +56,21 @@ namespace openspace {
 
         ViewPort _viewport{};
 
+
+
         std::vector<float> _schwarzschildWarpTable;
+        std::vector<float> flatDataStar;
+        KDTree _starKDTree{};
+
         std::unique_ptr<ghoul::opengl::BufferBinding<
-            ghoul::opengl::bufferbinding::Buffer::ShaderStorage>> _ssboDataBinding;
+            ghoul::opengl::bufferbinding::Buffer::ShaderStorage>> _ssboSchwarzschildDataBinding;
+        std::unique_ptr<ghoul::opengl::BufferBinding<
+            ghoul::opengl::bufferbinding::Buffer::ShaderStorage>> _ssboStarDataBinding;
 
         GLuint _quadVao = 0;
         GLuint _quadVbo = 0;
-        GLuint _ssboData = 0;
+        GLuint _ssboSchwarzschildWarpTable = 0;
+        GLuint _ssboStarKDTree = 0;
 
         UniformCache(environmentTexture, viewGrid, worldRotationMatrix, cameraRotationMatrix) _uniformCache;
 

modules/blackhole/shaders/blackhole_fs.glsl

@@ -15,12 +15,18 @@ layout (std430) buffer ssbo_warp_table {
   float schwarzschildWarpTable[];
 };
 
+layout (std430) buffer ssbo_star_map {
+    float starKDTree[];
+};
 
 const float PI = 3.1415926535897932384626433832795f;
 const float VIEWGRIDZ = -1.0f;
 const float INF = 1.0f/0.0f;
 
-// Math
+
+/**********************************************************
+                        Math
+***********************************************************/
 
 float lerp(float P0, float P1, float t) {
     return P0 + t * (P1 - P0);
@@ -34,16 +40,18 @@ float atan2(float a, float b){
     return 0.0f;
 }
 
-// Conversions
+/**********************************************************
+                        Conversions
+***********************************************************/
 
-vec2 cartesianToSpherical(vec3 cartisian) {
-    float theta = atan2(sqrt(cartisian.x * cartisian.x + cartisian.y * cartisian.y) , cartisian.z);
-    float phi = atan2(cartisian.y, cartisian.x);
+vec2 cartesianToSpherical(vec3 cartesian) {
+    float theta = atan2(sqrt(cartesian.x * cartesian.x + cartesian.y * cartesian.y) , cartesian.z);
+    float phi = atan2(cartesian.y, cartesian.x);
 
-    return vec2(phi, theta);
+    return vec2(theta, phi);
 }
 
-vec3 sphericalToCartesian(float phi, float theta){
+vec3 sphericalToCartesian(float theta, float phi){
     float x = sin(theta)*cos(phi);
     float y = sin(theta)*sin(phi);
     float z = cos(theta);
@@ -52,14 +60,15 @@ vec3 sphericalToCartesian(float phi, float theta){
 }
 
 vec2 sphericalToUV(vec2 sphereCoords){
-    float u = sphereCoords.x / (2.0f * PI) + 0.5f;
-    float v = mod(sphereCoords.y, PI) / PI;
+    float u = sphereCoords.y / (2.0f * PI) + 0.5f;
+    float v = mod(sphereCoords.x, PI) / PI;
     
     return vec2(u, v);
 }
 
-//Warp Table
-
+/**********************************************************
+                        Warp Table
+***********************************************************/
 ivec2 bstWarpTable(float phi){
     float midPhi = -1.0f;
     float deltaPhi = -1.0f;
@@ -129,13 +138,51 @@ float getEndAngleFromTable(float phi){
 }
 
 vec2 applyBlackHoleWarp(vec2 cameraOutSphereCoords){
-    float phi = cameraOutSphereCoords.x;
-    float theta = cameraOutSphereCoords.y;
+    float theta = cameraOutSphereCoords.x;
+    float phi = cameraOutSphereCoords.y;
     theta = getEndAngleFromTable(theta);
-    return vec2(phi, theta);
+    return vec2(theta, phi);
 }
 
-// Fragment shader function
+/**********************************************************
+                         Star Map
+***********************************************************/
+
+float angularDist(vec2 a, vec2 b) {
+    float dTheta = a.x - b.x;
+    float dPhi = a.y - b.y;
+    return sqrt(dTheta * dTheta + sin(a.x) * sin(b.x) * dPhi * dPhi);
+}
+
+vec4 searchNearestStar(vec3 sphericalCoords) {
+    const int NODE_SIZE = 6;
+    const int SIZE = starKDTree.length() / NODE_SIZE;
+    int index = 0;
+    int nodeIndex = 0;
+    int depth = 0;
+    int axis = -1;
+
+    while(index < SIZE && starKDTree[nodeIndex] > 0.0f){
+        if (angularDist(sphericalCoords.yz, vec2(starKDTree[nodeIndex + 1], starKDTree[nodeIndex + 2])) < 0.001f){
+                return vec4(0.9f, 0.9f, 0.8f, 0.2f);
+        }
+
+        axis = depth % 2 + 1;
+        if(sphericalCoords[axis] < starKDTree[nodeIndex + axis]){
+            index = 2 * index + 1;
+        } else {
+            index = 2 * index + 2;
+        }
+        nodeIndex = index * NODE_SIZE;
+        depth += 1;
+    }
+    return vec4(0.0f);
+}
+
+
+/**********************************************************
+                        Fragment shader
+***********************************************************/
 
 Fragment getFragment() {
     Fragment frag;
@@ -152,7 +199,7 @@ Fragment getFragment() {
     #if SHOW_BLACK_HOLE == 1
     // Apply black hole warping to spherical coordinates
     envMapSphericalCoords = applyBlackHoleWarp(sphericalCoords);
-    if (isnan(envMapSphericalCoords.y)) {
+    if (isnan(envMapSphericalCoords.x)) {
         // If inside the event horizon
         frag.color = vec4(0.0f);
         return frag;
@@ -181,6 +228,8 @@ Fragment getFragment() {
     vec2 uv = sphericalToUV(sphericalCoords);
     vec4 texColor = texture(environmentTexture, uv);
     
+    texColor = clamp(texColor + searchNearestStar(vec3(0.0f, sphericalCoords.x, sphericalCoords.y)), 0.f, 1.f);
+    
     frag.color = texColor;
     return frag;
 }