From f469a1d1e8413e07f3e3033e719161ffd17fe463 Mon Sep 17 00:00:00 2001 From: hachem Date: Fri, 22 Aug 2025 04:31:37 +0200 Subject: [add]: volumetric cloud rendering --- Assets/Shaders/Geodesic.glsl | 212 +++++++++++++++++++++++++++++++++++++---- src/Engine/Engine.cpp | 12 +-- src/Engine/Engine.h | 17 +++- src/States/SimulationState.cpp | 26 +++++ 4 files changed, 241 insertions(+), 26 deletions(-) diff --git a/Assets/Shaders/Geodesic.glsl b/Assets/Shaders/Geodesic.glsl index 9f871a7..ace58a5 100644 --- a/Assets/Shaders/Geodesic.glsl +++ b/Assets/Shaders/Geodesic.glsl @@ -20,6 +20,7 @@ layout(std140, binding = 2) uniform Disk float disk_r2; float disk_num; float thickness; + float disk_density; }; layout(std140, binding = 3) uniform Objects @@ -35,7 +36,7 @@ layout(std140, binding = 4) uniform Simulation int maxStepsMoving; int maxStepsStatic; float earlyExitDistance; - int _pad0; + float time; }; const float SagA_rs = 1.269e10; @@ -54,6 +55,112 @@ vec4 objectColor = vec4(0.0); vec3 hitCenter = vec3(0.0); float hitRadius = 0.0; +float hash(float p) +{ + p = fract(p * 0.1031); + p *= p + 33.33; + p *= p + p; + return fract(p); +} + +float hash(vec2 p) +{ + vec3 p3 = fract(vec3(p.xyx) * vec3(0.1031, 0.1030, 0.0973)); + p3 += dot(p3, p3.yzx + 33.33); + return fract((p3.x + p3.y) * p3.z); +} + +float hash(vec3 p) +{ + p = fract(p * vec3(0.1031, 0.1030, 0.0973)); + p += dot(p, p.yxz + 33.33); + return fract((p.x + p.y) * p.z); +} + +float noise(vec3 x) +{ + vec3 i = floor(x); + vec3 frac = fract(x); + + vec3 u = frac * frac * (3.0 - 2.0 * frac); + + float a = hash(i); + float b = hash(i + vec3(1.0, 0.0, 0.0)); + float c = hash(i + vec3(0.0, 1.0, 0.0)); + float d = hash(i + vec3(1.0, 1.0, 0.0)); + float e = hash(i + vec3(0.0, 0.0, 1.0)); + float f = hash(i + vec3(1.0, 0.0, 1.0)); + float g = hash(i + vec3(0.0, 1.0, 1.0)); + float h = hash(i + vec3(1.0, 1.0, 1.0)); + + return mix(mix(mix(a, b, u.x), mix(c, d, u.x), u.y), + mix(mix(e, f, u.x), mix(g, h, u.x), u.y), u.z); +} + +float fbm(vec3 x, int octaves) +{ + float v = 0.0; + float a = 0.5; + float f = 1.0; + vec3 shift = vec3(100, 200, 300); + + for (int i = 0; i < octaves; ++i) + { + v += a * noise(x * f); + x = x * 2.0 + shift; + a *= 0.5; + f *= 2.0; + } + return v; +} + +float GetCloudDensity(vec3 pos) +{ + float r_cyl = length(vec2(pos.x, pos.z)); + float r_norm = (r_cyl - disk_r1) / (disk_r2 - disk_r1); + + if (r_norm < 0.0 || r_norm > 1.0) + return 0.0; + + float h_norm = abs(pos.y) / thickness; + float vertical_falloff = exp(-h_norm * h_norm * 3.0); + float radial_density = 1.0 - r_norm * 0.5; + + vec3 noise_pos = pos * 1e-10; + float keplerian_speed = 1.0 / sqrt(r_norm + 0.1); + + float rotation_angle = time * keplerian_speed * 0.5; + vec3 rotated_pos = vec3( + pos.x * cos(rotation_angle) - pos.z * sin(rotation_angle), + pos.y, + pos.x * sin(rotation_angle) + pos.z * cos(rotation_angle) + ) * 1e-10; + + float large_turbulence = fbm(rotated_pos * 1.2, 5); + + float medium_wisps = fbm(rotated_pos * 2.5, 4); + float small_detail = fbm(rotated_pos * 6.0, 3); + float fine_detail = fbm(rotated_pos * 10.0, 2); + + float noise_mask = large_turbulence * 0.4 + + medium_wisps * 0.3 + + small_detail * 0.2 + + fine_detail * 0.1; + + noise_mask = smoothstep(0.25, 0.75, noise_mask); + + float angle = atan(pos.z, pos.x); + float rotated_angle = angle + time * 0.5; + + float spiral_arms = sin(rotated_angle * 3.0 + r_norm * 15.0) * 0.15 + 0.85; + + float orbital_angle = angle + time * keplerian_speed * 0.8; + float orbital_pattern = sin(orbital_angle * 2.0 + r_norm * 8.0) * 0.2 + 0.8; + + float density = vertical_falloff * radial_density * noise_mask * spiral_arms * orbital_pattern; + return density * disk_density; +} + struct Ray { float x, y, z; @@ -163,11 +270,62 @@ void RK4Step(inout Ray ray, float dL) ray.z = ray.r * cos(ray.theta); } -bool CrossesEquatorialPlane(vec3 oldPos, vec3 newPos) +bool IsInDiskVolume(vec3 pos) +{ + float r_cyl = length(vec2(pos.x, pos.z)); + return (r_cyl >= disk_r1 && r_cyl <= disk_r2 && abs(pos.y) <= thickness); +} + +vec4 SampleDiskColor(vec3 pos) { - bool crossed = (oldPos.y * newPos.y < 0.0); - float r = length(vec2(newPos.x, newPos.z)); - return crossed && (r >= disk_r1 && r <= disk_r2); + float r_cyl = length(vec2(pos.x, pos.z)); + float r_norm = (r_cyl - disk_r1) / (disk_r2 - disk_r1); + + vec3 innerColor = vec3(1.0, 0.9, 0.5); + vec3 midColor = vec3(1.0, 0.6, 0.2); + vec3 outerColor = vec3(0.9, 0.3, 0.1); + + vec3 baseColor; + if (r_norm < 0.5) + baseColor = mix(innerColor, midColor, r_norm * 2.0); + else + baseColor = mix(midColor, outerColor, (r_norm - 0.5) * 2.0); + + float r_norm_rot = (r_cyl - disk_r1) / (disk_r2 - disk_r1); + float keplerian_speed = 1.0 / sqrt(r_norm_rot + 0.1); + + vec3 noise_pos = pos * 1e-10; + + float color_rotation_angle = time * keplerian_speed * 0.3; + vec3 rotated_color_pos = vec3( + pos.x * cos(color_rotation_angle) - pos.z * sin(color_rotation_angle), + pos.y, + pos.x * sin(color_rotation_angle) + pos.z * cos(color_rotation_angle) + ) * 1e-10; + + float large_color = fbm(rotated_color_pos * 1.8, 4); + float medium_color = fbm(rotated_color_pos * 4.0, 3); + float small_color = fbm(rotated_color_pos * 8.0, 2); + float colorVariation = (large_color * 0.5 + medium_color * 0.3 + small_color * 0.2) * 0.6; + baseColor = baseColor * (1.0 + colorVariation); + + float density = GetCloudDensity(pos); + vec3 brightness_noise_pos = pos * 1e-10; + + float brightness_rotation_angle = time * keplerian_speed * 0.7; + vec3 rotated_brightness_pos = vec3( + pos.x * cos(brightness_rotation_angle) - pos.z * sin(brightness_rotation_angle), + pos.y, + pos.x * sin(brightness_rotation_angle) + pos.z * cos(brightness_rotation_angle) + ) * 1e-10; + + float brightness_large = fbm(rotated_brightness_pos * 3.0, 3); + float brightness_medium = fbm(rotated_brightness_pos * 5.0, 2); + float brightness_small = fbm(rotated_brightness_pos * 7.0, 2); + float brightness_noise = (brightness_large * 0.6 + brightness_medium * 0.3 + brightness_small * 0.1); + float brightness = 1.0 + density * 1.0 + brightness_noise * 0.4; + + return vec4(baseColor * brightness, density); } float CalculateAdaptiveStepSize(Ray ray, float baseStepSize) @@ -203,8 +361,10 @@ void main() float lambda = 0.0; bool hitBlackHole = false; - bool hitDisk = false; bool hitObject = false; + + vec4 accumulatedColor = vec4(0.0); + float transmittance = 1.0; int maxSteps = cam.moving ? maxStepsMoving : maxStepsStatic; @@ -246,10 +406,29 @@ void main() vec3 newPos = vec3(ray.x, ray.y, ray.z); - if (CrossesEquatorialPlane(prevPos, newPos)) - { - hitDisk = true; - break; + if (IsInDiskVolume(newPos)) + { + vec4 diskSample = SampleDiskColor(newPos); + float density = diskSample.a; + vec3 diskColor = diskSample.rgb; + + float stepLength = currentStepSize * 1e-8; + float absorption = density * stepLength * 1.2; + float scattering = density * stepLength * 0.8; + float extinction = absorption + scattering; + + float stepTransmittance = exp(-extinction); + + vec3 emission = diskColor * density * stepLength * 2.5 * sqrt(disk_density); + accumulatedColor.rgb += emission * transmittance; + + transmittance *= stepTransmittance; + + if (transmittance < 0.01) + { + accumulatedColor.a = 1.0 - transmittance; + break; + } } if (i % objectCheckInterval == 0 && InterceptObject(ray)) @@ -264,13 +443,9 @@ void main() break; } - if (hitDisk) - { - double r = length(vec3(ray.x, ray.y, ray.z)) / disk_r2; - vec3 diskColor = vec3(1.0, r, 0.2); - color = vec4(diskColor, r); - - } else if (hitBlackHole) + accumulatedColor.a = 1.0 - transmittance; + + if (hitBlackHole) color = vec4(0.0, 0.0, 0.0, 1.0); else if (hitObject) { @@ -284,8 +459,9 @@ void main() vec3 shaded = objectColor.rgb * intensity; color = vec4(shaded, objectColor.a); + color = mix(accumulatedColor, color, color.a); } else - color = vec4(0.0); + color = accumulatedColor; imageStore(outImage, pix, color); } \ No newline at end of file diff --git a/src/Engine/Engine.cpp b/src/Engine/Engine.cpp index 245f6fd..37cafd8 100644 --- a/src/Engine/Engine.cpp +++ b/src/Engine/Engine.cpp @@ -36,14 +36,14 @@ namespace Donut m_ComputeProgram = CreateComputeProgram("Assets/Shaders/Geodesic.glsl"); m_CameraUBO = UniformBuffer::Create(128, 1); - m_DiskUBO = UniformBuffer::Create(sizeof(float) * 4, 2); + m_DiskUBO = UniformBuffer::Create(sizeof(float) * 5, 2); // Added density parameter uint32_t objUBOSize = sizeof(int) + 3 * sizeof(float) + 16 * (sizeof(glm::vec4) + sizeof(glm::vec4)) + 16 * sizeof(float); m_ObjectsUBO = UniformBuffer::Create(objUBOSize, 3); - m_SimulationUBO = UniformBuffer::Create(sizeof(int) * 2 + sizeof(float) * 2, 4); + m_SimulationUBO = UniformBuffer::Create(sizeof(int) * 2 + sizeof(float) * 2, 4); // time added to existing float auto result = QuadVAO(); m_QuadVAO = result.first; @@ -177,8 +177,8 @@ namespace Donut float r1 = static_cast(m_SagA.m_Rs * 2.2); float r2 = static_cast(m_SagA.m_Rs * 5.2); float num = 2.0f; - float thickness = 1e9f; - float diskData[4] = { r1, r2, num, thickness }; + float thickness = static_cast(m_SagA.m_Rs * m_DiskThickness); + float diskData[5] = { r1, r2, num, thickness, m_DiskDensity }; m_DiskUBO->SetData(diskData, sizeof(diskData)); m_DiskUBO->Bind(2); @@ -191,13 +191,13 @@ namespace Donut int maxStepsMoving; int maxStepsStatic; float earlyExitDistance; - int padding; + float time; } data; data.maxStepsMoving = m_MaxStepsMoving; data.maxStepsStatic = m_MaxStepsStatic; data.earlyExitDistance = m_EarlyExitDistance; - data.padding = 0; + data.time = static_cast(glfwGetTime()) * m_RotationSpeed; m_SimulationUBO->SetData(&data, sizeof(data)); m_SimulationUBO->Bind(4); diff --git a/src/Engine/Engine.h b/src/Engine/Engine.h index 328435a..a0b25b5 100644 --- a/src/Engine/Engine.h +++ b/src/Engine/Engine.h @@ -100,6 +100,15 @@ namespace Donut void SetMaxStepsStatic(int steps) { m_MaxStepsStatic = steps; } void SetEarlyExitDistance(float distance) { m_EarlyExitDistance = distance; } + float GetDiskThickness() const { return m_DiskThickness; } + void SetDiskThickness(float thickness) { m_DiskThickness = thickness; } + + float GetDiskDensity() const { return m_DiskDensity; } + void SetDiskDensity(float density) { m_DiskDensity = density; } + + float GetRotationSpeed() const { return m_RotationSpeed; } + void SetRotationSpeed(float speed) { m_RotationSpeed = speed; } + void LoadObjectsFromScene(const std::vector& objects); void ExportHighResFrame(const std::string& filename, int width = 4096, int height = 3072); void PrintObjectInfo() const; @@ -131,8 +140,12 @@ namespace Donut Camera m_Camera; bool m_Gravity = false; - int m_MaxStepsMoving = 60000; - int m_MaxStepsStatic = 30000; + int m_MaxStepsMoving = 60000; + int m_MaxStepsStatic = 30000; float m_EarlyExitDistance = 5e11f; + + float m_DiskThickness = 0.1; + float m_DiskDensity = 1.2f; + float m_RotationSpeed = 1.0f; }; }; diff --git a/src/States/SimulationState.cpp b/src/States/SimulationState.cpp index d961f9e..28504e6 100644 --- a/src/States/SimulationState.cpp +++ b/src/States/SimulationState.cpp @@ -220,6 +220,32 @@ namespace Donut ImGui::Spacing(); + ImGui::TextColored(ImVec4(0.9f, 0.9f, 1.0f, 1.0f), "Accretion Disk"); + ImGui::Separator(); + + float diskThickness = engine.GetDiskThickness(); + if (ImGui::SliderFloat("Cloud Thickness", &diskThickness, 0.1f, 2.0f, "%.2f")) + { + engine.SetDiskThickness(diskThickness); + } + ImGui::TextDisabled("Thickness relative to Schwarzschild radius"); + + float diskDensity = engine.GetDiskDensity(); + if (ImGui::SliderFloat("Cloud Density", &diskDensity, 0.1f, 3.0f, "%.2f")) + { + engine.SetDiskDensity(diskDensity); + } + ImGui::TextDisabled("Overall density multiplier"); + + float rotationSpeed = engine.GetRotationSpeed(); + if (ImGui::SliderFloat("Rotation Speed", &rotationSpeed, 0.0f, 3.0f, "%.2f")) + { + engine.SetRotationSpeed(rotationSpeed); + } + ImGui::TextDisabled("Rotation speed multiplier (0 = no rotation)"); + + ImGui::Spacing(); + ImGui::TextColored(ImVec4(0.9f, 0.9f, 1.0f, 1.0f), "Camera"); ImGui::Separator(); -- cgit v1.3