[add]: volumetric cloud rendering

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<float>(m_SagA.m_Rs * 2.2);
         float r2 = static_cast<float>(m_SagA.m_Rs * 5.2);
         float num = 2.0f;
-        float thickness = 1e9f;
-        float diskData[4] = { r1, r2, num, thickness };
+        float thickness = static_cast<float>(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<float>(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<Donut::Object>& 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();