summaryrefslogtreecommitdiffstats
path: root/harnack_tracing.cginc
blob: bd911faa8a8b233402c8ad12c46b32dfb5503c02 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
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