From cfe824bddc89dfbffc82f05f2bd156f673571d8a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Wilhelm=20Bj=C3=B6rkstr=C3=B6m?=
 <143391787+Grantallkotten@users.noreply.github.com>
Date: Mon, 10 Mar 2025 14:49:53 +0100
Subject: [PATCH] Interpolation between angels for camera to env-map

Co-Authored-By: Emil Wallberg <49481622+EmilWallberg@users.noreply.github.com>
---
 modules/blackhole/cuda/blackhole_cuda.cu      |  4 ++
 .../blackhole/rendering/renderableblackhole.h |  2 +-
 modules/blackhole/shaders/blackhole_fs.glsl   | 52 +++++++++++++++----
 3 files changed, 47 insertions(+), 11 deletions(-)

diff --git a/modules/blackhole/cuda/blackhole_cuda.cu b/modules/blackhole/cuda/blackhole_cuda.cu
index eec6517af4..a75b767309 100644
--- a/modules/blackhole/cuda/blackhole_cuda.cu
+++ b/modules/blackhole/cuda/blackhole_cuda.cu
@@ -85,6 +85,10 @@ void schwarzchild(
 
     cudaMemcpy(angle_out, d_angle_values, num_rays * 2 * sizeof(float), cudaMemcpyDeviceToHost);
 
+    // Add handeling of special case straight backwards
+    angle_out[(num_rays - 1) * 2] = 0.0f;
+    angle_out[(num_rays - 1) * 2 + 1] = 0.0f;
+
     cudaFree(d_dudphi_0_values);
     cudaFree(d_angle_values);
 }
diff --git a/modules/blackhole/rendering/renderableblackhole.h b/modules/blackhole/rendering/renderableblackhole.h
index 9e1f56c7da..3255514197 100644
--- a/modules/blackhole/rendering/renderableblackhole.h
+++ b/modules/blackhole/rendering/renderableblackhole.h
@@ -41,7 +41,7 @@ namespace openspace {
 
         ghoul::opengl::ProgramObject* _program = nullptr;
 
-        size_t _rayCount = 10000;
+        size_t _rayCount = 1000;
         size_t _stepsCount = 50000;
         float _stepLength = 0.001f;
 
diff --git a/modules/blackhole/shaders/blackhole_fs.glsl b/modules/blackhole/shaders/blackhole_fs.glsl
index 5f49feef10..1b65ae4eb4 100644
--- a/modules/blackhole/shaders/blackhole_fs.glsl
+++ b/modules/blackhole/shaders/blackhole_fs.glsl
@@ -19,25 +19,43 @@ const float PI = 3.1415926535897932384626433832795f;
 const float VIEWGRIDZ = -1.0f;
 const float INF = 1.0f/0.0f;
 
-float bstWarpTable(float phi){
+float lerp(float start, float end, float t) {
+    return start + t * (end - start);
+}
+
+float interpelateWarpTable(int indexStart, int indexEnd, float localPhi){
+    float envMapPhiStart = schwarzschildWarpTable[indexStart * 2 + 1];
+    float envMapPhiEnd = schwarzschildWarpTable[indexEnd * 2 + 1];
+
+    float localPhiStart = schwarzschildWarpTable[indexStart * 2];
+    float localPhiEnd = schwarzschildWarpTable[indexEnd * 2];
+
+    float t = (localPhi - localPhiStart) / (localPhiEnd - localPhiStart);
+    
+    return lerp(envMapPhiStart, envMapPhiEnd, t);
+}
+
+ivec2 bstWarpTable(float phi){
     float midPhi = -1.0f;
     float deltaPhi = -1.0f;
 
     float minDeltaPhi = INF;
-    int index = -1;
+    int closestIndex = -1;
 
     int left = 0;
     int mid = -1;
-    int right = schwarzschildWarpTable.length() / 2 - 1;
+    const int tableSize = schwarzschildWarpTable.length() / 2;
+    int right = tableSize - 1;
 
     while(left <= right){
         mid = (left + right) / 2;
         midPhi = schwarzschildWarpTable[mid * 2];
 
         deltaPhi = abs(midPhi - phi);
+        
         if(deltaPhi < minDeltaPhi){
             minDeltaPhi = deltaPhi;
-            index = mid;
+            closestIndex = mid;
         }
 
         if (phi > midPhi) {
@@ -46,12 +64,30 @@ float bstWarpTable(float phi){
             left = mid + 1;
         }
     }
+    
+    int leftIndex = closestIndex - 1;
+    int rightIndex = closestIndex + 1;
+    float leftDist = INF; 
+    float rightDist = INF;
 
-    return (index != -1) ? schwarzschildWarpTable[index * 2 + 1] : 0.0f;
+    if (rightIndex < tableSize - 1){
+        rightDist = abs(schwarzschildWarpTable[rightIndex * 2] - phi);
+    }
+    if (leftIndex > 0){
+        leftDist = abs(schwarzschildWarpTable[leftIndex * 2] - phi);
+    }
+
+    int nextClosestIndex = (rightDist < leftDist) ? rightIndex : leftIndex;
+
+    float v1 = schwarzschildWarpTable[closestIndex * 2 + 1];
+    float v2 = schwarzschildWarpTable[nextClosestIndex * 2 + 1];
+
+    return v1 < v2 ? ivec2(closestIndex, nextClosestIndex) : ivec2(nextClosestIndex, closestIndex);
 }
 
 float getEndAngleFromTable(float phi){
-    return bstWarpTable(phi);
+    ivec2 indices = bstWarpTable(phi);
+    return interpelateWarpTable(indices.x, indices.y, phi);
 }
 
 vec2 applyBlackHoleWarp(vec2 cameraOutSphereCoords){
@@ -92,10 +128,6 @@ vec2 sphericalToUV(vec2 sphereCoords){
     return vec2(u, v);
 }
 
-float lerp(float start, float end, float t) {
-    return start + t * (end - start);
-}
-
 Fragment getFragment() {
     Fragment frag;