begin sketching out harnack tracer

1 files changed, 168 insertions, 0 deletions

diff --git a/harnack_tracing.cginc b/harnack_tracing.cginc
new file mode 100644
index 0000000..bd911fa
--- /dev/null
+++ b/harnack_tracing.cginc
@@ -0,0 +1,168 @@
+#ifndef __HARNACK_TRACING_INC

+#define __HARNACK_TRACING_INC

+

+#include "cnlohr.cginc"

+#include "globals.cginc"

+#include "interpolators.cginc"

+

+#if defined(_HARNACK_TRACING)

+

+#define MAX_ITERATIONS 100

+#define UNIT_SHIFT 5.0

+#define WALL_THICKNESS 0.1

+

+#define SPHERE_RADIUS 0.49

+#define OUTER_RADIUS (SPHERE_RADIUS + 0.5)

+

+float gyroid(float3 pos);

+float3 gyroid_gradient(float3 pos);

+bool harnackTrace(float3 ro, float3 rd, out float t, out float3 pos, out float3 normal, float tMax);

+bool intersectSphere(float3 ro, float3 rd, float3 center, float radius, out float t0, out float t1);

+float getRadius(float3 p);

+float getMaxStep4D(float fx, float R, float levelset, float shift);

+bool closeToLevelset(float f, float levelset, float tol, float gradNorm);

+bool betweenLevelsets(float f, float loBound, float hiBound, float tol, float gradNorm);

+

+struct HarnackTracingOutput {

+  float4 color;

+  float3 worldPos;  // intersection point in world space

+  float3 normal;    // normal in world space

+};

+

+// Gyroid function implementation

+float gyroid(float3 pos) {

+    float timeOffset = 2.0 * PI * _Harnack_Tracing_Gyroid_Speed * _Time.y;

+    float3 p = _Harnack_Tracing_Gyroid_Scale * pos + float3(0.0, timeOffset, 0.0);

+    return sin(p.x) * cos(p.y) + sin(p.y) * cos(p.z) + sin(p.z) * cos(p.x);

+}

+

+// Gradient of gyroid function

+float3 gyroid_gradient(float3 pos) {

+    float timeOffset = 2.0 * PI * _Harnack_Tracing_Gyroid_Speed * _Time.y;

+    float3 p = _Harnack_Tracing_Gyroid_Scale * pos + float3(0.0, timeOffset, 0.0);

+    return float3(

+        cos(p.x) * cos(p.y) - sin(p.z) * sin(p.x),

+        cos(p.y) * cos(p.z) - sin(p.x) * sin(p.y),

+        cos(p.z) * cos(p.x) - sin(p.y) * sin(p.z)

+    ) * _Harnack_Tracing_Gyroid_Scale;

+}

+

+// Get radius from point to outer boundary

+float getRadius(float3 p) {

+    return OUTER_RADIUS - length(p);

+}

+

+// Calculate if point is close to levelset

+bool closeToLevelset(float f, float levelset, float tol, float gradNorm) {

+    return abs(f - levelset) < tol;

+}

+

+// Calculate if point is between levelsets (for wall thickness)

+bool betweenLevelsets(float f, float loBound, float hiBound, float tol, float gradNorm) {

+    return max(loBound - f, f - hiBound) < tol;

+}

+

+// Sphere intersection test

+bool intersectSphere(float3 ro, float3 rd, float3 center, float radius, out float t0, out float t1) {

+    float3 oc = ro - center;

+    float b = dot(oc, rd);

+    float c = dot(oc, oc) - radius * radius;

+    float h = b * b - c;

+    

+    if (h < 0.0) return false;

+    

+    h = sqrt(h);

+    t0 = -b - h;

+    t1 = -b + h;

+    

+    return true;

+}

+

+// Calculate maximum step size for Harnack tracing

+float getMaxStep4D(float fx, float R, float levelset, float shift) {

+    float a = (fx + shift) / (levelset + shift);

+    float u = pow(3.0 * sqrt(3.0 * pow(a, 3.0) + 81.0 * pow(a, 2.0)) + 27.0 * a, 1.0 / 3.0);

+    return R * abs(u / 3.0 - a / u - 1.0);

+}

+

+// Main Harnack tracing function

+bool harnackTrace(float3 ro, float3 rd, out float t, out float3 pos, out float3 normal, float tMax) {

+    t = 0.0;

+    float levelset = sin(_Harnack_Tracing_Gyroid_Speed * _Time.y);

+    

+    // Early sphere intersection test

+    float t0, t1;

+    if (!intersectSphere(ro, rd, float3(0,0,0), SPHERE_RADIUS, t0, t1) || tMax < 0.0) 

+        return false;

+    

+    // Optimize bounds

+    t = max(t0, 0.0);

+    tMax = min(t1, tMax);

+    

+    // Check immediate intersection

+    pos = ro + t0 * rd;

+    float val = gyroid(pos);

+    float3 gradF = gyroid_gradient(pos);

+    if (betweenLevelsets(val, levelset - WALL_THICKNESS, levelset + WALL_THICKNESS, 0.025, length(gradF))) {

+        normal = normalize(gradF);

+        return true;

+    }

+

+    float t_overstep = 0.0;

+    

+    [loop]

+    for (int iters = 0; iters < MAX_ITERATIONS && t < tMax; iters++) {

+        pos = ro + t * rd + t_overstep * rd;

+        

+        val = gyroid(pos);

+        gradF = gyroid_gradient(pos);

+        

+        float offset_levelset = clamp(val, levelset - WALL_THICKNESS, levelset + WALL_THICKNESS);

+        

+        float R = getRadius(pos);

+        float shift = exp(sqrt(2.0) * R) * UNIT_SHIFT;

+        float r = getMaxStep4D(val, R, offset_levelset, shift);

+        

+        if (r >= t_overstep && closeToLevelset(val, offset_levelset, 0.025, length(gradF))) {

+            normal = normalize(gradF);

+            return true;

+        }

+        

+        float stepSize = (r >= t_overstep) ? t_overstep + r : 0.0;

+        t_overstep = (r >= t_overstep) ? r * 0.75 : 0.0;

+        t += stepSize;

+    }

+     

+    normal = float3(0,0,0);

+    return false;

+}

+

+HarnackTracingOutput HarnackTracing(v2f i) {

+  HarnackTracingOutput o = (HarnackTracingOutput)0;

+  

+#if defined(_HARNACK_TRACING_GYROID)

+  float3 worldRo = _WorldSpaceCameraPos;

+  float3 ro = mul(unity_WorldToObject, float4(worldRo, 1.0)).xyz;

+  

+  float3 worldRd = normalize(i.worldPos - worldRo);

+  float3 rd = normalize(mul((float3x3)unity_WorldToObject, worldRd));

+  

+  float t;

+  float3 pos;

+  float3 normal;

+  

+  bool hit = harnackTrace(ro, rd, t, pos, normal, 100.0);

+  

+  if (hit) {

+    o.color = 1;

+    o.worldPos = mul(unity_ObjectToWorld, float4(pos, 1.0)).xyz;

+    o.normal = normalize(mul(transpose((float3x3)unity_WorldToObject), normal));

+  }

+#endif

+  

+  return o;

+}

+

+#endif  // _HARNACK_TRACING

+

+#endif  // __HARNACK_TRACING_INC