From 481c4ffed1200cf5b65c88f3fade2cb8280914ea Mon Sep 17 00:00:00 2001
From: Niklas Korz <korz.niklask@gmail.com>
Date: Fri, 17 Sep 2021 07:50:48 +0000
Subject: [PATCH] Adaptive integration depending on gradient

---
 src/main_view.rs | 2 +-
 src/main_view.wgsl | 36 ++++++++++++++++++++++++++++++++----
 src/overlay.wgsl | 9 ++++-----
 3 file(s) changed, 37 insertion(s)(+), 10 deletion(s)(-)

diff --git a/src/main_view.rs b/src/main_view.rs
--- a/src/main_view.rs
+++ b/src/main_view.rs
@@ -673,7 +673,7 @@
         let original_shader_src = self.shader_src.clone();
         let high_accuracy_shader_src = self
             .shader_src
-            .replace("let h: f32 = 0.05;", "let h: f32 = 0.001;")
+            .replace("let h: f32 = 0.1;", "let h: f32 = 0.001;")
             .replace("let steps: i32 = 100;", "let steps: i32 = 5000;");
         self.reload_shader(
             device,
diff --git a/src/main_view.wgsl b/src/main_view.wgsl
--- a/src/main_view.wgsl
+++ b/src/main_view.wgsl
@@ -218,8 +218,8 @@
     return vec4<f32>(0.0, 0.0, 0.0, 0.0);
 }
 
-let h: f32 = 0.05;
-let steps: i32 = 100;
+let h: f32 = 0.1;
+let steps: i32 = 5000;
 
 struct NonlinearRayColorResult {
     color: vec4<f32>;
@@ -236,14 +236,23 @@
     var cur_point: vec3<f32> = start_point;
     var cur_dir: vec3<f32> = start_dir;
     var t: f32 = 0.0;
+    var last_v: vec3<f32> = vec3<f32>(0.0, 0.0, 0.0);
+    var last_diff: f32 = -1.0;
+    let h_initial = h;
+    var h: f32 = h_initial;
 
-    for (var i: i32 = 0; i < steps; i = i + 1) {
+    for (; t <= 5.0;) {
         // Runge-Kutta method
         let k1 = field_function(cur_point, cur_point, start_dir, cur_dir, t);
         let k2 = field_function(cur_point, cur_point + 0.5 * h * k1, start_dir, k1, t + 0.5 * h);
         let k3 = field_function(cur_point, cur_point + 0.5 * h * k2, start_dir, k2, t + 0.5 * h);
         let k4 = field_function(cur_point, cur_point + h * k3, start_dir, k3, t + h);
         let v = (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0;
+        let diff = length(v - last_v);
+        if (last_diff >= 0.0 && diff > 10.0 * last_diff && h > 0.002) {
+            h = 0.001;
+            continue;
+        }
         cur_dir = (1.0 - field_weight) * cur_dir + field_weight * v;
 
         let step_dir = cur_dir * h;
@@ -254,6 +263,11 @@
 
         cur_point = cur_point + step_dir;
         t = t + h;
+        last_v = v;
+        if (2.0 * diff < last_diff && h < h_initial) {
+            h = h_initial;
+        }
+        last_diff = diff;
     }
 
     result.mapping_point = vec4<f32>(cur_point, 1.0);
@@ -270,19 +284,28 @@
     var cur_dir: vec3<f32> = start_dir;
     var color: vec4<f32>;
     var t: f32 = 0.0;
+    var last_v: vec3<f32> = vec3<f32>(0.0, 0.0, 0.0);
+    var last_diff: f32 = -1.0;
+    let h_initial = h;
+    var h: f32 = h_initial;
 
     var sample: RaySample;
     sample.color = vec4<f32>(sample_color, 0.5);
     let sample_steps = 100;
     let sample_step_size = steps / sample_steps;
 
-    for (var i: i32 = 0; i < steps; i = i + 1) {
+    for (var i: i32 = 0; i < steps && t <= 5.0; i = i + 1) {
         // Runge-Kutta method
         let k1 = field_function(cur_point, cur_point, start_dir, cur_dir, t);
         let k2 = field_function(cur_point, cur_point + 0.5 * h * k1, start_dir, k1, t + 0.5 * h);
         let k3 = field_function(cur_point, cur_point + 0.5 * h * k2, start_dir, k2, t + 0.5 * h);
         let k4 = field_function(cur_point, cur_point + h * k3, start_dir, k3, t + h);
         let v = (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0;
+        let diff = length(v - last_v);
+        if (last_diff >= 0.0 && diff > 10.0 * last_diff && h > 0.002) {
+            h = 0.001;
+            continue;
+        }
         cur_dir = (1.0 - field_weight) * cur_dir + field_weight * v;
 
         let step_dir = cur_dir * h;
@@ -295,6 +318,11 @@
 
         cur_point = cur_point + step_dir;
         t = t + h;
+        last_v = v;
+        if (2.0 * diff < last_diff && h < h_initial) {
+            h = h_initial;
+        }
+        last_diff = diff;
     }
 }
 
diff --git a/src/overlay.wgsl b/src/overlay.wgsl
--- a/src/overlay.wgsl
+++ b/src/overlay.wgsl
@@ -180,12 +180,11 @@
     if (settings.overlay_mode == 1) {
         let exponent = lyapunov_exponent(coords);
         let scaled_exponent = exp(settings.lyapunov_scaling * exponent - 5.0);
-        return scaled_exponent * vec4<f32>(1.0, 1.0, 1.0, 1.0) + (1.0 - scaled_exponent) * color;
-    } elseif (settings.overlay_mode == 2) {
-        if (ridge_test(coords, size)) {
-            return vec4<f32>(1.0, 1.0, 1.0, 1.0);
+        if (scaled_exponent > 0.8) {
+            return scaled_exponent * vec4<f32>(1.0, 1.0, 1.0, 1.0) + (1.0 - scaled_exponent) * color;
         }
-        return color;
+    } elseif (settings.overlay_mode == 2 && ridge_test(coords, size)) {
+        return vec4<f32>(1.0, 1.0, 1.0, 1.0);
     }
     return color;
 }

--
tangled.sh