mirror/OpenSpace
1
0
Fork 0
mirror of https://github.com/OpenSpace/OpenSpace.git synced 2026-01-04 02:29:49 -06:00
Code Issues Packages Projects Releases Wiki Activity

WIP: Fixed using wrong domain in shader

Browse Source 
Warptable used angles 0 to 2pi in phi while the shader assumed -pi to pi
This commit is contained in:
Emil Wallberg
2025-04-24 21:05:46 +02:00
parent a108f01b75
commit 8116151d0c
4 changed files with 65 additions and 239 deletions
Download Patch File Download Diff File Expand all files Collapse all files

15
modules/blackhole/cuda/kerr.cu
View File

@@ -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;
}
}

2
modules/blackhole/rendering/renderableblackhole.cpp
View File

@@ -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;
}

2
modules/blackhole/rendering/renderableblackhole.h
View File

@@ -46,7 +46,7 @@ namespace openspace {
std::unique_ptr<ghoul::opengl::ProgramObject> _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;

285
modules/blackhole/shaders/kerr_fs.glsl
View File

@@ -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;