diff --git a/modules/blackhole/cuda/kerr.cu b/modules/blackhole/cuda/kerr.cu index 36f9dadb0d..8c981cf0c0 100644 --- a/modules/blackhole/cuda/kerr.cu +++ b/modules/blackhole/cuda/kerr.cu @@ -22,7 +22,7 @@ __constant__ unsigned int c_num_steps = 15000; __constant__ unsigned int c_layers = 1; __constant__ float c_M = 1.0f; // Mass parameter __constant__ float c_epsilon = 1e-10; // Numerical tolerance -__constant__ float3 worldUp = { 0.0f, 1.0f, 0.0f }; +__constant__ float3 worldUp = { 0.0f, 0.0f, 1.0f }; @@ -310,29 +310,34 @@ __global__ void simulateRayKernel(float3 pos, size_t num_rays_per_dim, float* lo int idx_entry = 0; entry[idx_entry] = theta; - entry[idx_entry++ + 1] = phi; + entry[idx_entry + 1] = phi; + idx_entry += 2; for (int step = 0; step < c_num_steps; step++) { // Terminate integration if ray is inside the horizon or outside the environment. if (y[0] < 2.f) { while (idx_entry <= c_layers) { entry[idx_entry] = nanf(""); - entry[idx_entry++ + 1] = nanf(""); + entry[idx_entry + 1] = nanf(""); + idx_entry += 2; } break; } else if (y[0] > 100.f) { //TODO Check collision with the correct env map and save to entry entry[idx_entry] = y[1]; - entry[idx_entry++ + 1] = y[2]; + entry[idx_entry + 1] = y[2]; + idx_entry += 2; if (idx_entry > c_layers) { break; } } // Advance one RK4 step. rk4(y, c_h, E, L, k_val); + } if (idx_entry <= c_layers) { entry[idx_entry] = y[1]; - entry[idx_entry++ + 1] = y[2]; + entry[idx_entry + 1] = y[2]; + idx_entry += 2; } } diff --git a/modules/blackhole/rendering/renderableblackhole.cpp b/modules/blackhole/rendering/renderableblackhole.cpp index feeee8dc25..e3f55b2841 100644 --- a/modules/blackhole/rendering/renderableblackhole.cpp +++ b/modules/blackhole/rendering/renderableblackhole.cpp @@ -138,7 +138,7 @@ namespace openspace { if (glm::distance(cameraPosition, _chacedCameraPos) > 0.01f && rerender) { //kerr(2e11, 0, 0, _rs, 0.99f, _rayCount, _stepsCount, _schwarzschildWarpTable); - kerr(cameraPosition.x, cameraPosition.y, cameraPosition.z, _rs, 0.5f, _rayCount, _stepsCount, _schwarzschildWarpTable); + kerr(cameraPosition.x, cameraPosition.y, cameraPosition.z, _rs, 0.99f, _rayCount, _stepsCount, _schwarzschildWarpTable); _chacedCameraPos = cameraPosition; //rerender = false; } diff --git a/modules/blackhole/rendering/renderableblackhole.h b/modules/blackhole/rendering/renderableblackhole.h index 925fe92201..42962fce79 100644 --- a/modules/blackhole/rendering/renderableblackhole.h +++ b/modules/blackhole/rendering/renderableblackhole.h @@ -46,7 +46,7 @@ namespace openspace { std::unique_ptr _cullProgram = nullptr; glm::dvec3 _chachedTranslation{}; glm::dvec3 _chacedCameraPos{}; - static constexpr size_t _rayCount{ 100 }; + static constexpr size_t _rayCount{ 250 }; static constexpr size_t _stepsCount = 100000; static constexpr float _stepLength = 0.0001f; static constexpr size_t _mapCount = 3; diff --git a/modules/blackhole/shaders/kerr_fs.glsl b/modules/blackhole/shaders/kerr_fs.glsl index 36ba3852af..8184d524ff 100644 --- a/modules/blackhole/shaders/kerr_fs.glsl +++ b/modules/blackhole/shaders/kerr_fs.glsl @@ -7,40 +7,29 @@ in vec2 TexCoord; uniform sampler2D environmentTexture; uniform sampler2D viewGrid; uniform sampler1D colorBVMap; -uniform mat4 cameraRotationMatrix;;; +uniform mat4 cameraRotationMatrix; uniform mat4 worldRotationMatrix; layout (std430) buffer ssbo_warp_table { float schwarzschildWarpTable[]; }; -// layout (std430) buffer ssbo_star_map { -// float starMapKDTrees[]; -// }; - -// layout (std430) buffer ssbo_star_map_start_indices { -// int starMapKDTreesIndices[]; -// }; - const float PI = 3.1415926535897932384626433832795f; const float VIEWGRIDZ = -1.0f; const float INF = 1.0f/0.0f; const float LUM_LOWER_CAP = 0.01; const int NODE_SIZE = 6; -// const int StarArrySize = starMapKDTrees.length() / NODE_SIZE; -// const int layerCount = starMapKDTreesIndices.length(); const int layerCount = 1; -// const int tableNodeSize = starMapKDTreesIndices.length() + 1; const float starRadius = 0.002f; const int STACKSIZE = 32; const int tableSize = schwarzschildWarpTable.length() / 2; const int num_rays = schwarzschildWarpTable.length() / (layerCount + 1); - const float max_theta = PI; -const float max_phi = 2*PI; -const int num_rays_per_dim = 100; +const float max_phi = 2.0 * PI; +const int num_rays_per_dim = 250; + /********************************************************** Math ***********************************************************/ @@ -62,9 +51,14 @@ float atan2(float a, float b){ ***********************************************************/ vec2 cartesianToSpherical(vec3 cartesian) { - float theta = atan2(sqrt(cartesian.x * cartesian.x + cartesian.y * cartesian.y) , cartesian.z); + float theta = atan2(sqrt(cartesian.x * cartesian.x + cartesian.y * cartesian.y), cartesian.z); float phi = atan2(cartesian.y, cartesian.x); + // Remap phi from [-PI, +PI] to [0, 2*PI) + if (phi < 0.0) { + phi += 2.0 * PI; + } + return vec2(theta, phi); } @@ -77,84 +71,58 @@ vec3 sphericalToCartesian(float theta, float phi){ } vec2 sphericalToUV(vec2 sphereCoords){ - float u = sphereCoords.y / (2.0f * PI) + 0.5f; - float v = mod(sphereCoords.x, PI) / PI; - + float u = sphereCoords.y / (2.0f * PI); // phi ∈ [0, 2π] → u ∈ [0, 1] + float v = sphereCoords.x / PI; // theta ∈ [0, π] → v ∈ [0, 1] + return vec2(u, v); } - /********************************************************** Warp Table ***********************************************************/ vec2 getInterpolatedWarpAngles(float input_theta, float input_phi, int layer) { - // 1) Convert to float grid indices - float theta_f = (input_theta) - * float(num_rays_per_dim - 1) - / (max_theta); - float phi_f = input_phi - * float(num_rays_per_dim -1) - / (2.0 * PI); + float phi = input_phi; // already in [0, 2PI) + float theta = input_theta; - // 2) Integer bounds - int theta_low = int(clamp(floor(theta_f), 0.0, float(num_rays_per_dim-1))); - int theta_high = int(clamp(theta_low + 1, 0, num_rays_per_dim-1)); + float theta_f = clamp(theta / max_theta * float(num_rays_per_dim - 1), + 0.0, float(num_rays_per_dim - 1)); + float phi_f = clamp(phi / max_phi * float(num_rays_per_dim - 1), + 0.0, float(num_rays_per_dim - 1)); - int phi_low = int(floor(phi_f)) % num_rays_per_dim; - int phi_high = (phi_low + 1) % num_rays_per_dim; + int theta0 = int(floor(theta_f)); + int phi0 = int(floor(phi_f)); + int theta1 = min(theta0 + 1, num_rays_per_dim - 1); + int phi1 = (phi0 + 1) % num_rays_per_dim; - // 3) Flat indices - int N = num_rays_per_dim; - int idx00 = theta_low * N + phi_low; - int idx10 = theta_high * N + phi_low; - int idx01 = theta_low * N + phi_high; - int idx11 = theta_high * N + phi_high; + float t = fract(theta_f); + float u = fract(phi_f); - // 4) Record stride & per-layer offset - int record_stride = 2 + layerCount * 2; - int offset = 2 + layer * 2; + int N = num_rays_per_dim; + int stride = 2 + layerCount * 2; + int off = 2 + layer * 2; - // 5) Base pointers into the table - int base00 = idx00 * record_stride; - int base10 = idx10 * record_stride; - int base01 = idx01 * record_stride; - int base11 = idx11 * record_stride; + #define FETCH(t_idx, p_idx) \ + vec2( \ + schwarzschildWarpTable[((t_idx)*N + (p_idx))*stride + off], \ + schwarzschildWarpTable[((t_idx)*N + (p_idx))*stride + off + 1] \ + ) - // 6) Extract the original theta/phi stored in the table - float theta00 = schwarzschildWarpTable[base00 ]; - float phi00 = schwarzschildWarpTable[base00 + 1]; - float theta10 = schwarzschildWarpTable[base10 ]; - float phi01 = schwarzschildWarpTable[base01 + 1]; + vec2 v00 = FETCH(theta0, phi0); + vec2 v10 = FETCH(theta1, phi0); + vec2 v01 = FETCH(theta0, phi1); + vec2 v11 = FETCH(theta1, phi1); - // 7) Interpolation weights - float t = (input_theta - theta00) / (theta10 - theta00); - float u = (input_phi - phi00) / (phi01 - phi00); + bool beyoundTheHorizion = isnan(v00.x) || isnan(v10.x) || isnan(v01.x) || isnan(v11.x); + if (beyoundTheHorizion) { + return vec2(0/0); + } - // 8) Fetch the four target warp‐angle pairs - vec2 v00 = vec2( - schwarzschildWarpTable[base00 + offset ], - schwarzschildWarpTable[base00 + offset + 1] - ); - vec2 v10 = vec2( - schwarzschildWarpTable[base10 + offset ], - schwarzschildWarpTable[base10 + offset + 1] - ); - vec2 v01 = vec2( - schwarzschildWarpTable[base01 + offset ], - schwarzschildWarpTable[base01 + offset + 1] - ); - vec2 v11 = vec2( - schwarzschildWarpTable[base11 + offset ], - schwarzschildWarpTable[base11 + offset + 1] - ); - return v00; - // 9) Bilinear interpolate - vec2 interp_low = mix(v00, v10, t); - vec2 interp_high = mix(v01, v11, t); - return mix(interp_low, interp_high, u); + vec2 low = mix(v00, v10, t); + vec2 high = mix(v01, v11, t); + return mix(low, high, u); } - +#undef FETCH vec2 applyBlackHoleWarp(vec2 cameraOutAngles, int layer){ float theta = cameraOutAngles.x; @@ -163,179 +131,32 @@ vec2 applyBlackHoleWarp(vec2 cameraOutAngles, int layer){ } /********************************************************** - Star Map -***********************************************************/ - -vec3 BVIndex2rgb(float color) { - // BV is [-0.4, 2.0] - float st = (color + 0.4) / (2.0 + 0.4); - - return texture(colorBVMap, st).rgb; -} - -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); -} - -// Search a single KD-tree given its start and node count. -// Returns a vec4 where rgb is the accumulated star color and a is the total alpha. -// vec4 searchTree(int treeStart, int treeNodeCount, vec3 sphericalCoords) { -// // Local struct for tree traversal. -// struct TreeIndex { -// int index; // local index (starting at 0 for this tree) -// int depth; -// }; - -// // Initialize candidate tracking for this tree. -// float bestDist = INF; -// int bestLocalIndex = -1; - -// TreeIndex stack[STACKSIZE]; -// int stackIndex = 0; -// TreeIndex cur; -// cur.index = 0; // root of the tree (local index 0) -// cur.depth = 0; - -// // Perform iterative search on the single KD-tree. -// while (true) { -// // If current index is out-of-bounds, backtrack if possible. -// if (cur.index < 0 || cur.index >= treeNodeCount) { -// if (stackIndex > 0) { -// cur = stack[--stackIndex]; -// continue; -// } else { -// break; -// } -// } - -// // Convert local index to absolute index in the flat array. -// int absIndex = treeStart + cur.index; -// int base = absIndex * NODE_SIZE; - -// float d = angularDist(sphericalCoords.yz, -// vec2(starMapKDTrees[base + 1], starMapKDTrees[base + 2])); -// if (d < bestDist) { -// bestDist = d; -// bestLocalIndex = cur.index; -// } - -// // Determine axis (1 or 2) using depth (sphericalCoords: index 1=theta, index 2=phi) -// int axis = cur.depth % 2 + 1; -// float diff = sphericalCoords[axis] - starMapKDTrees[base + axis]; - -// int offset = int(diff >= 0.0); -// int closerLocal = 2 * cur.index + 1 + offset; -// int fartherLocal = 2 * cur.index + 2 - offset; - -// int nextDepth = cur.depth + 1; - -// // If the farther branch might hold a closer point, push it onto the stack. -// if (abs(diff) < bestDist && (fartherLocal < treeNodeCount)) { -// stack[stackIndex].index = fartherLocal; -// stack[stackIndex].depth = nextDepth; -// stackIndex++; -// } - -// // Continue down the closer branch. -// cur.index = closerLocal; -// cur.depth = nextDepth; -// } // End of tree traversal loop - -// if (bestLocalIndex == -1 || bestDist >= starRadius) { -// return vec4(0.0); -// } - -// int bestAbsIndex = treeStart + bestLocalIndex; -// int base = bestAbsIndex * NODE_SIZE; - -// float observedDistance = starMapKDTrees[base]; -// float luminosity = pow(10.0, 1.89 - 0.4 * starMapKDTrees[base + 4]); -// luminosity /= pow(observedDistance, 1.1); -// luminosity = max(luminosity, LUM_LOWER_CAP); - -// float alpha = 1.0 - pow(bestDist / starRadius, 2.2); -// vec3 starColor = BVIndex2rgb(starMapKDTrees[base + 3]); - -// // Return the star contribution (with pre-multiplied alpha). -// return vec4(starColor * alpha * luminosity, alpha * luminosity); -// } - -// vec4 searchNearestStar(vec3 sphericalCoords, int layer) { -// vec4 accumulatedColor = vec4(0.0); -// float totalAlpha = 0.0; - -// int kdTreeCount = starMapKDTreesIndices.length(); -// int totalNodes = starMapKDTrees.length() / NODE_SIZE; - -// // Determine the start offset for the current tree -// int treeStartFloat = starMapKDTreesIndices[layer]; -// int treeStart = treeStartFloat / NODE_SIZE; - -// // Determine the end of this tree. -// int treeEnd; -// if (layer < kdTreeCount - 1) { -// int nextTreeStartFloat = int(starMapKDTreesIndices[layer+1]); -// treeEnd = nextTreeStartFloat / NODE_SIZE; -// } else { -// treeEnd = totalNodes; -// } -// int treeNodeCount = treeEnd - treeStart; - -// // Search the current tree. -// return searchTree(treeStart, treeNodeCount, sphericalCoords); -// } - - -/********************************************************** - Fragment shader + Fragment shader ***********************************************************/ Fragment getFragment() { Fragment frag; vec4 viewCoords = normalize(vec4(texture(viewGrid, TexCoord).xy, VIEWGRIDZ, 0.0f)); - - // User local input rotation of the black hole vec4 rotatedViewCoords = cameraRotationMatrix * viewCoords; - - vec2 sphericalCoords; - - vec4 accumulatedColor = vec4(0.0f); - float accumulatedWeight = 0.0f; // Track total weight of blending - // Apply black hole warping to spherical coordinates + vec2 sphericalCoords; + vec4 accumulatedColor = vec4(0.0f); + float accumulatedWeight = 0.0f; + for (int l = 0; l < layerCount; ++l) { sphericalCoords = cartesianToSpherical(rotatedViewCoords.xyz); sphericalCoords = applyBlackHoleWarp(sphericalCoords, l); - if (isnan(sphericalCoords.x)) { - // If inside the event horizon - frag.color = vec4(.5f); + + if (sphericalCoords == vec2(0.0f/0.0f)) { + frag.color = vec4(0.0f); return frag; } - - // vec4 envMapCoords = vec4(sphericalToCartesian(sphericalCoords.x, sphericalCoords.y), 0.0f); - - // User world input rotation of the black hole - // envMapCoords = worldRotationMatrix * envMapCoords; - - // sphericalCoords = cartesianToSpherical(envMapCoords.xyz); - // vec4 starColor = searchNearestStar(vec3(0.0f, sphericalCoords.x, sphericalCoords.y), l); - - // if (starColor.a > 0.0) { - // float layerWeight = exp(-0.5 * l); // Earlier layers have more weight - - // // Blend using weighted alpha blending - // accumulatedColor.rgb = (accumulatedColor.rgb * accumulatedWeight + starColor.rgb * starColor.a * layerWeight) / (accumulatedWeight + starColor.a * layerWeight); - // accumulatedWeight += starColor.a * layerWeight; - // } } vec2 uv = sphericalToUV(sphericalCoords); vec4 texColor = texture(environmentTexture, uv); - // frag.color.rgb = accumulatedColor.rgb * accumulatedWeight + texColor.rgb * (1.0 - accumulatedWeight); frag.color.rgb = texColor.rgb; frag.color.a = 1.0; return frag;