Sphere tracing [ag-0010]
✍️source

    Definition 1. Point-to-set distance [hart1996sphere, def. 1] [ag-000R]

      The point-to-set distance defines the distance from a point \(\boldsymbol {x} \in \mathbb {R}^3\) to a set \(A \subset \mathbb {R}^3\) as the distance from \(x\) to the closest point in \(A\), \[ d(\boldsymbol {x}, A)=\min _{y \in A}\lVert \boldsymbol {x}-\boldsymbol {y}\rVert \]

      Definition 2. Signed distance bound [hart1996sphere, def. 2] [ag-000S]

        A function \(f: \mathbb {R}^3 \rightarrow \mathbb {R}\) is a signed distance bound of its implicit surface \(f^{-1}(0)\) if and only if \[ |f(\boldsymbol {x})| \leq d\left (\boldsymbol {x}, f^{-1}(0)\right ) \] If equality holds, then \(f\) is a signed distance function.

        Definition 3. Lipschitz bound [hart1996sphere, def. 3] [ag-000P]

          A function \(f: \mathbb {R}^3 \rightarrow \mathbb {R}\) is Lipschitz over a domain \(D\) if and only if for all \(\boldsymbol {x}, \boldsymbol {y} \in D\), there exists a positive finite constant \(\lambda \), called Lipschitz bound, such that \[ |f(\boldsymbol {x})-f(\boldsymbol {y})| \leq \lambda \lVert \boldsymbol {x}-\boldsymbol {y}\rVert . \] The Lipschitz constant, denoted \(\operatorname {Lip} f\), is the minimum \(\lambda \) satisfying the condition above.

          Theorem 4.  [gillespie2024ray, thm. 1] [ag-000T]

            Let \(f\) be Lipschitz with Lipschitz bound \(\lambda \geq \) Lip \(f\). Then the function \(f / \lambda \) is a signed distance bound of its implicit surface.

            Definition 5. Sphere tracing [hart1996sphere, sec. 2.3, eq. 12] [ag-001I]

              The sphere tracing algorithm is to march a ray by an adaptive safe step size, which is the absolute value of the signed distance bound calculated by \(F/\lambda \) per Theorem 4, the rest is the same as ray marching (naïve).

              Formally, the root found by sphere tracing is the limit point of the sequence defined by the recurrence equation \[t_{i+1} = t_i + \frac {|F(t_i)|}{\lambda } = t_i + \frac {f(\boldsymbol {r}(t_i))}{\lambda } \] where \(F\) is the ray intersection, and \(f\) is the SDF.

              Usually \(\lambda = 1\), and \(F(t_i)\) stays positive until the stopping condition \(F(t_i) < \epsilon \) is met, so the above can be simplified to \[t_{i+1} = t_i + F(t_i) = t_i + f(\boldsymbol {r}(t_i)) \]

              Figure 6. Sphere tracing [ag-001H]

                Inspired by [coulon2022ray, fig. 2.2]

                Figure 7. Raymarching in raymarching (sphere tracing) [ag-001J]

                  #define MAT_SCREEN 0.
#define MAT_SPHERE1 1.
#define MAT_FLOOR  2.
#define MAT_CORNER 3.
#define MAT_SCREENZ 4.
#define MAT_MARCHSPHERE 5.
#define MAT_MARCHROUTE 6.
#define MAT_MARCHRADIUS 7.
#define MAT_SPHERE2 8.
#define MAT_SPHERE3 9.
#define MAT_RAYDIRECTION 10.
#define MAT_HITPOINT 11.

float uvToP = 0.0;
float colToUv = 1.0;
float screenZ = 2.5;

float sphereAnim = 1.0;
float radiusAnim = 1.0;
float routeAnim = 1.0;

float raySphereAlpha = 0.0;

float cornersAlpha = 0.0;
float cornersAnim = 0.0;

float screenZAlpha = 0.0;

float radiusAlpha = 1.0;

float rayDirectionAnim = 0.0;
float rayDirectionAnim2 = 0.0;

float screenAlpha = 0.0;

float hitSphereID = 0.0;

// =====================Camera========================
vec3 cameraPos, cameraTarget;

//================================

float sum = 0.0;

float tl(float val, float offset, float range)
{
    float im = sum + offset;
    float ix = im + range;
    sum += offset + range;
    return clamp((val - im) / (ix - im), 0.0, 1.0);
}

float cio(float t) {
	return t < 0.5
	? 0.5 * (1.0 - sqrt(1.0 - 4.0 * t * t))
	: 0.5 * (sqrt((3.0 - 2.0 * t) * (2.0 * t - 1.0)) + 1.0);
}

float eio(float t) {
	return t == 0.0 || t == 1.0
		? t
		: t < 0.5
		? +0.5 * pow(2.0, (20.0 * t) - 10.0)
		: -0.5 * pow(2.0, 10.0 - (t * 20.0)) + 1.0;
}

float fio(float t) {
	return t == 0.0 || t == 1.0
		? t
		: t <= 0.5
		? -0.5 * pow(2.0, 5.0 - ((t + 0.5) * 10.0)) + 1.0 - 0.5
		: +0.5 * pow(2.0, (10.0 * (t - 0.5)) - 5.0) + 0.5;
}

void timeLine(float time)
{
    cameraPos = vec3(5., 5., -3.); // mix(vec3(0.0), , eio(t));
    cameraTarget = vec3(0.0, 0.0, 5.);
    
    //time += 32.;
    float t; // = tl(time, 1.0, 0.5);
    // uvToP = mix(0.0, 1.0, eio(t));
    
    // t = tl(time, 1.0, 1.0);
    
    // t = tl(time, 0.5, 1.0);
    // raySphereAlpha = mix(0., 1., t);
    // cornersAlpha = mix(0., 1., t);
    
    // t = tl(time, 0.5, 1.0);
    // cornersAnim = mix(0., 1., t);
    
    t = tl(time, 0.0, 0.5);
    // cameraPos = mix(cameraPos, vec3(5., 6., -3.) , t); // eio(t));
    ///cameraTarget = vec3(0.0, 0.0, 3.0); //mix(cameraTarget, , eio(t));
    
    // t = tl(time, 0.5, 0.5);
    // screenZAlpha = mix(0., 1., t);

	// t = tl(time, 0.5, .5);
    // colToUv = mix(colToUv, 0.0, eio(t));
    
    // t = tl(time, 0.5, 1.0);
    // screenZ = mix(screenZ, 5.0, eio(t));
    
    // t = tl(time, 1.0, 1.0);
    // screenZ = mix(screenZ, .75, eio(t));
    
    
    // t = tl(time, 1.0, 1.0);
    // screenZ = mix(screenZ, 2.5, eio(t));
    
    
    // t = tl(time, 0.5, 1.0);
    // cornersAlpha = mix(cornersAlpha, 0., t);
    // screenZAlpha = mix(screenZAlpha, 0., t);
    // colToUv = mix(colToUv, 1.0, eio(t));
    
    // t = tl(time, 0.0, 0.0);
    // cornersAnim = mix(cornersAnim, 0., t);
    
    // t = tl(time, 0.5, 1.0);
    // cameraPos = mix(cameraPos, vec3(5., 15., 6.0), eio(t));
    // cameraTarget = mix(cameraTarget, vec3(0.0, 0.0, 6.), eio(t));
    
    raySphereAlpha = 1.;
    
    for (int i=0; i<20; i++) {
        t = tl(time, 0., 0.5);
        radiusAnim = mix(radiusAnim, float(i) + 2., t);
        t = tl(time, 0., 0.5);
        routeAnim = mix(routeAnim, float(i) + 2., t);
        t = tl(time, 0., 0.5);
        sphereAnim = mix(sphereAnim, float(i) + 2., t);
    }
    
    
    // t = tl(time, 1.0, 1.0);
    // cameraPos = mix(cameraPos, vec3(2., 3., -3.), eio(t));
    // cameraTarget = mix(cameraTarget, vec3(0.0, 0.0, 3.0), eio(t));
    // radiusAlpha = mix(radiusAlpha, 0.0, t);
    // routeAnim = mix(routeAnim, 1.0, t);
    // sphereAnim = mix(sphereAnim, 1.0, t);
    // colToUv = mix(colToUv, 1.0, eio(t));
    // screenAlpha = mix(screenAlpha, 0.0, t);
    
    t = tl(time, 0.5, 0.0);
    radiusAnim = mix(radiusAnim, 0.0, t);
    routeAnim = mix(routeAnim, 0.0, t);
    sphereAnim = mix(sphereAnim, 0.0, t);

    
    t = tl(time, 0.5, 1.0);
    rayDirectionAnim2 = mix(rayDirectionAnim2, 1.0, t); // fio(t)); // eio(t));
    
    t = tl(time, 0.5, 3.0);
    rayDirectionAnim = mix(rayDirectionAnim, 1.0, fio(t));
    
    t = tl(time, 3.5, 0.0);
    rayDirectionAnim2 = mix(rayDirectionAnim2, 0.0, t);
    raySphereAlpha = mix(raySphereAlpha, 0.0, t);
    
    // rayDirectionAnim = mix(rayDirectionAnim, 0.0, t);
    // raySphereAlpha = mix(raySphereAlpha, 0.0, t);
    // rayDirectionAnim2 = mix(rayDirectionAnim2, 0.0, t);
    
    
    // t = tl(time, 0.5, 1.0);
    // cameraPos = mix(cameraPos, vec3(0.), eio(t));
    // cameraTarget = mix(cameraTarget, vec3(0.0, 0.0, 5.), eio(t));
}

vec2 opU(vec2 a, vec2 b)
{
    return a.x < b.x ? a : b;
}

float sdLine( vec3 p, vec3 a, vec3 b, float r )
{
    vec3 pa = p - a, ba = b - a;
    float h = clamp( dot(pa,ba)/dot(ba,ba), 0.0, 1.0 );
    return length( pa - ba*h ) - r;
}

float sdBox( vec3 p, vec3 b )
{
  vec3 d = abs(p) - b;
  return length(max(d,0.0))
         + min(max(d.x,max(d.y,d.z)),0.0); // remove this line for an only partially signed sdf 
}

float sdSphere(vec3 p, float s)
{
    return length(p) - s;
}

vec3 sunDir = normalize(vec3(.3, .25, .2));
vec3 skyColor(vec3 rd)
{
    float sundot = clamp(dot(rd,sunDir),0.0,1.0);
    // sky		
    vec3 col = mix(vec3(0.2,0.5,0.85)*1.1, vec3(0.0,0.15,0.7), rd.y);
    col = mix( col, 0.85*vec3(0.8,0.8,0.7), pow( 1.0-max(rd.y,0.0), 4.0 ) );
    // sun
    col += 0.3*vec3(1.0,0.7,0.4)*pow( sundot,5.0 );
    col += 0.5*vec3(1.0,0.8,0.6)*pow( sundot,64.0 );
    col += 6.0*vec3(1.0,0.8,0.6)*pow( sundot,1024.0 );
    return col * 3.0;
}

vec3 screenPos;

vec2 sceneSpheres(vec3 p)
{
    vec2 s1 = vec2(sdSphere(p - vec3(-2.0, 0.0, 6.0), 1.), MAT_SPHERE1);
    vec2 s2 = vec2(sdSphere(p - vec3(0.0, 0.0, 16.0), 1.), MAT_SPHERE2);
    vec2 s3 = vec2(sdSphere(p - vec3(2.0, 0.0, 9.0), 1.), MAT_SPHERE3);
    return opU(opU(s1, s2), s3);
}

vec2 sceneMap(vec3 p)
{
    vec2 s = sceneSpheres(p);
    vec2 p1 = vec2(p.y + 2.0, MAT_FLOOR);
    return opU(s, p1);
}

vec3 normal(vec3 p)
{
    vec2 e = vec2(0.001, 0.0);
    float d = sceneMap(p).x;
    return normalize(d - vec3(sceneMap(p - e.xyy).x, sceneMap(p - e.yxy).x, sceneMap(p - e.yyx).x));
}

void sceneTrace(inout vec3 pos, vec3 ray, out vec2 mat, inout float depth, float maxD)
{
    vec3 ro = pos;
    int i = 0;
    for(; i < 70; i++) {
        if (depth > maxD) {
        	depth = maxD;
            break;
        }
    	pos = ro + ray * depth;
        mat = sceneMap(pos);
        if (mat.x < 0.001) {
            break;
        }
        depth += mat.x;
    }
}

vec2 screenMap(vec3 p)
{
    vec2 screenSize = iResolution.xy / min(iResolution.x, iResolution.y);
    vec2 b1 = vec2(sdBox(p - vec3(0., 0., screenZ), vec3(screenSize, 0.01)), MAT_SCREEN);
    return b1;
}

void screenTrace(inout vec3 pos, vec3 ray, out vec2 mat, inout float depth, float maxD)
{
    vec3 ro = pos;
    for(int i = 0; i < 20; i++) {
        if (depth > maxD) {
        	depth = maxD;
            break;
        }
    	pos = ro + ray * depth;
        mat = screenMap(pos);
        if (mat.x < 0.001) {
            break;
        }
        depth += mat.x;
    }
}

float remap(float val, float im, float ix, float om, float ox)
{
    return clamp(om + (val - im) * (ox - om) / (ix - im), om, ox);
}

vec2 gizmoCorners(vec3 p)
{
    vec2 screenSize = iResolution.xy / min(iResolution.x, iResolution.y);
    
    float a1 = remap(cornersAnim, 0.0, 0.25, 0.0, 1.0);
    float a2 = remap(cornersAnim, 0.25, 0.5, 0.0, 1.0);
    float a3 = remap(cornersAnim, 0.5, 0.75, 0.0, 1.0);
    float a4 = remap(cornersAnim, 0.75, 1.0, 0.0, 1.0);
    
    vec2 c1 = vec2(sdLine(p, vec3(0.), mix(vec3(0.), vec3(screenSize, screenZ), a1), 0.02), MAT_CORNER);
    vec2 c2 = vec2(sdLine(p, vec3(0.), mix(vec3(0.), vec3(screenSize * vec2(1.,-1.), screenZ), a2), 0.02), MAT_CORNER);
    vec2 c3 = vec2(sdLine(p, vec3(0.), mix(vec3(0.), vec3(screenSize * vec2(-1.,-1.), screenZ), a3), 0.02), MAT_CORNER);
    vec2 c4 = vec2(sdLine(p, vec3(0.), mix(vec3(0.), vec3(screenSize * vec2(-1.,1.), screenZ), a4), 0.02), MAT_CORNER);
    return opU(c1, opU(c2, opU(c3, c4)));
}

vec2 gizmoScreenZ(vec3 p)
{
    vec2 c1 = vec2(sdLine(p, vec3(0.), vec3(0., 0., screenZ), 0.03), MAT_SCREENZ);
    return c1;
}

float sphereID = 0.0;

vec3 rayRoute = vec3(0., 1., 0.); // y must be init to non-zero
vec2 gizmoMarching(vec3 p)
{
    vec3 ray = vec3(0., 0., 1.);
    vec2 d = vec2(10000.);
    
    if(rayDirectionAnim2 > 0.0) {
        if(abs(rayRoute.y - 0.) <= 0.01) {
            ray = normalize(vec3(rayRoute.x, rayRoute.y, rayRoute.z));
            radiusAnim = 20.;
            routeAnim = 20.;
            sphereAnim = 20.;
        }
    }
    
    float t = 0.0;
    vec3 pos;
    int maxSteps = hitSphereID == 0. ? 20 : int(hitSphereID) + 1;
    for(int i=0; i<maxSteps; i++) {
    	pos = ray * t;
        vec2 s = vec2(sdSphere(p - pos, 0.15), MAT_MARCHSPHERE);
        if (s.x < d.x) {
            d = s;
            sphereID = float(i);   
        }
        
        float dist = sceneSpheres(pos).x;
        
        float anim = clamp(routeAnim - float(i) - 1., 0.0, 1.0);
        vec2 c1 = vec2(sdLine(p, pos, pos + mix(vec3(0.), ray * dist, anim), 0.03), MAT_MARCHROUTE);
        
        d = opU(d, c1);
        t += dist;
        
    }    
    return d;
}

vec2 gizmoMarchingRadius(vec3 p)
{
    vec2 d = vec2(p.y, MAT_MARCHRADIUS);
    return d;
}


vec2 gizmoRayDirection(vec3 p)
{
    float a1 = fract(rayDirectionAnim * 19.99999);
    float a2 = floor(rayDirectionAnim * 19.99999) / 20.;
    vec2 screenSize = iResolution.xy / min(iResolution.x, iResolution.y);
    screenSize.y *= -1.;
    screenSize = mix(-screenSize, screenSize, vec2(a1, a2));
    rayRoute = mix(vec3(0.), vec3(screenSize, screenZ) * 10.0, rayDirectionAnim2);
    vec2 c1 = vec2(sdLine(p, vec3(0.), mix(vec3(0.), vec3(screenSize, screenZ) * 10.0, rayDirectionAnim2), 0.02), MAT_RAYDIRECTION);
    
    vec3 ray = normalize(vec3(screenSize, screenZ));
    vec3 pos;
    float t = 0.01;
    for(int i = 0; i < 40; i++) {
    	pos = ray * t;
        vec2 d = sceneMap(pos);
        if (d.x < 0.001) {
            break;
        }
        t += d.x;
    }
    
    c1 = opU(c1, vec2(sdSphere(p - pos, 0.15), MAT_HITPOINT));
    
    return c1;
}

vec2 gizmoMap(vec3 p)
{
    vec2 d = opU(gizmoCorners(p), gizmoScreenZ(p));
    d = opU(d, gizmoMarching(p));
    d = opU(d, gizmoRayDirection(p));
    return d;
}

void gizmoTrace(inout vec3 pos, vec3 ray, out vec2 mat, inout float depth, float maxD)
{
    vec3 ro = pos;
    for(int i = 0; i < 60; i++) {
        if (depth > maxD) {
        	depth = maxD;
            break;
        }
    	pos = ro + ray * depth;
        mat = gizmoMap(pos);
        if (mat.x < 0.001) {
            break;
        }
        depth += mat.x;
    }
}

vec2 radiusMap(vec3 p)
{
    return gizmoMarchingRadius(p);
}

void radiusTrace(inout vec3 pos, vec3 ray, out vec2 mat, inout float depth, float maxD)
{
    vec3 ro = pos;
    for(int i = 0; i < 40; i++) {
        if (depth > maxD) {
        	depth = maxD;
            break;
        }
    	pos = ro + ray * depth;
        mat = radiusMap(pos);
        if (mat.x < 0.001 || depth > maxD) {
            break;
        }
        depth += mat.x;
    }
}

float shadow(in vec3 p, in vec3 l, float ma)
{
    float t = 0.1;
    float t_max = ma;
    
    float res = 1.0;
    for (int i = 0; i < 16; ++i)
    {
        if (t > t_max) break;
        vec3 pos = p + t*l;
        float d = opU(sceneMap(pos), screenMap(pos)).x; // sceneMap(pos).x; //opU(, screenMap(pos)).x;
        if (d < 0.001)
        {
            return 0.0;
        }
        t += d*1.0;
        res = min(res, 10.0 * d / t);
    }
    
    return res;
}

// checkerbord
// https://www.shadertoy.com/view/XlcSz2
float checkersTextureGradBox( in vec2 p, in vec2 ddx, in vec2 ddy )
{
    // filter kernel
    vec2 w = max(abs(ddx), abs(ddy)) + 0.01;  
    // analytical integral (box filter)
    vec2 i = 2.0*(abs(fract((p-0.5*w)/2.0)-0.5)-abs(fract((p+0.5*w)/2.0)-0.5))/w;
    // xor pattern
    return 0.5 - 0.5*i.x*i.y;                  
}

vec3 sceneShade(vec2 mat, vec3 pos, vec3 ray, float depth, float maxD)
{
    vec3 col;
    vec3 sky = skyColor(ray);
    if (depth > maxD - 0.01) {
    	return sky;
    }
    float sha = shadow(pos, sunDir, 10.);
    vec3 norm = normal(pos);
    vec3 albedo = vec3(0.);
    if(mat.y == MAT_SPHERE1) {
        albedo = vec3(1., 0., 0.);
    } else if (mat.y == MAT_SPHERE2) {
        albedo = vec3(0., 1., 0.);
    } else if (mat.y == MAT_SPHERE3) {
        albedo = vec3(0., 0., 1.);
    } else if(mat.y == MAT_FLOOR) {
        vec2 ddx_uvw = dFdx( pos.xz ); 
        vec2 ddy_uvw = dFdy( pos.xz ); 
        float checker = checkersTextureGradBox( pos.xz, ddy_uvw, ddy_uvw );
        albedo = vec3(max(0.2, checker)) * vec3(.8,0.8,0.7) * 2.0;
    }
    
    float diffuse = clamp(dot(norm, sunDir), 0.0, 1.0) * sha * 2.0;
    col = albedo * (diffuse + 0.05);
    
    float fo = 1.0 - exp2(-0.0001 * depth * depth);
    vec3 fco = 0.65*vec3(0.4,0.65,1.0);
    col = mix( col, sky, fo );
    return col;
}

vec3 screenShade(vec2 mat, vec3 pos)
{
    vec3 ro = vec3(0., 0., 0.);
    vec3 ta = vec3(0., 0., 3.);
    vec3 fo = normalize(ta - ro);
    vec3 ri = normalize(cross(vec3(0.,1.,0.), fo));
    vec3 up = normalize(cross(fo,ri));
    mat3 cam = mat3(ri,up,fo);
    vec3 ray = cam * normalize(pos);
    float depth = 0.01;
    vec3 p = ro;
    sceneTrace(p, ray, mat, depth, 100.);
    vec3 col = vec3(0.);
    col = sceneShade(mat, p, ray, depth, 100.);
    
    float a1 = fract(rayDirectionAnim * 19.99999);
    float a2 = floor(rayDirectionAnim * 19.99999 + 1.0) / 20.;
    float a3 = floor(rayDirectionAnim * 19.99999) / 20.;
    
    float aspect = iResolution.y / iResolution.x;
    float halfAspect = aspect * 0.5;
    
    a1 = step(pos.x * halfAspect + 0.5, a1);
    a2 = step(-pos.y * 0.49 + 0.51+ 0.0, a2);
    a3 = step(-pos.y * 0.49 + 0.51 + 0.0, a3);
    //col *= min(screenAlpha + min(a2 * a1 + a3, 1.0), 1.0);
    
    // vec3 uvCoord = vec3(pow(clamp(mix(pos.xy*0.5+0.5, pos.xy, uvToP), 0.0, 1.0), vec2(2.2)), 0.0);
    vec3 uvCoord = vec3(1.);
    col = mix(col, uvCoord, colToUv - min(screenAlpha + min(a2 * a1 + a3, 1.0), 1.0));
    return col;
    //return vec3(pow(clamp(pos.xy, 0.0, 1.0), vec2(2.2)), 0.0);
}

vec4 gizmoShade(vec2 mat, vec3 p)
{
    vec4 col = vec4(0.);

    if(mat.y == MAT_CORNER) {
        col = vec4(1.,0.,0., cornersAlpha);
    } else if (mat.y == MAT_SCREENZ) {
    	col = vec4(0.05, 0.05, 1., screenZAlpha);
    } else if (mat.y == MAT_MARCHSPHERE) {
        float alpha = clamp(sphereAnim - sphereID, 0.0, 1.0);
        vec3 sc = mix(vec3(.0, .1, 3.), vec3(1.02, 1., 0.02), float(sphereID == 0. || sphereID == hitSphereID));
        // vec3 sc = mix(vec3(.0, .1, 3.), vec3(.02, 1., 1.02), float(sphereID == 0. || sphereID == hitSphereID));
        float inRange = 1.0; //sphereID <= hitSphereID ? 1.0 : 0.0;
    	col = vec4(sc, alpha * raySphereAlpha * inRange);
    } else if (mat.y == MAT_MARCHROUTE) {
        col = vec4(1., 0., 0., .9);
    } else if (mat.y == MAT_RAYDIRECTION) {
        col = vec4(1., 0., 0., .9);
    } else if (mat.y == MAT_HITPOINT) {
    	col = vec4(1.02, 1., .02, raySphereAlpha);
    }
    return col;
}

vec4 radiusShade(vec2 mat, vec3 p)
{
    vec4 col = vec4(0.);

    vec3 ray = vec3(0., 0., 1.);
    vec2 d = vec2(10000.);
    
    if(rayDirectionAnim2 > 0.0) {
        if(abs(rayRoute.y - 0.) <= 0.01) {
            ray = normalize(vec3(rayRoute.x, rayRoute.y, rayRoute.z));
            radiusAnim = 20.;
            routeAnim = 20.;
            sphereAnim = 20.;
        }
    }
    
    float t = 0.0;
    for(int i=0; i<20; i++) {
        vec3 pos = ray * t;
        vec2 dd = sceneSpheres(pos);
        d = vec2(sdSphere(p - pos, dd.x), MAT_MARCHSPHERE);

        vec2 sceneDist = sceneMap(pos);
        if (sceneDist.x < 0.001) {
            if(hitSphereID == 0.) {
              hitSphereID = float(i);
            }
            break;
        }

        float alpha2 = step(radiusAnim, float(i) + 2.);
        float alpha = clamp(radiusAnim - float(i) - 1., 0.0, 1.0);
        
        vec4 cirCol = vec4(.0, 0.05, 0.1, 0.9); //mix(vec4(.0, 0.05, 0.1, 0.9), vec4(0.2, 1., 1.4, .6) , alpha2);
        // fill
        col = mix(col, cirCol, smoothstep(0.01, 0., d.x) * cirCol.a * alpha);
        // strike
        col = mix(col, vec4(0., 0., 0., 1.), smoothstep(0.02, 0., abs(d.x) - 0.01) * alpha);
        t += dd.x;
    }
    col *= radiusAlpha;
    return col;
}

float luminance2(vec3 col)
{
    return dot(vec3(0.298912, 0.586611, 0.114478), col);
}

vec3 reinhard(vec3 col, float exposure, float white) {
    col *= exposure;
    white *= exposure;
    float lum = luminance2(col);
    return (col * (lum / (white * white) + 1.0) / (lum + 1.0));
}

vec3 render(vec2 p) {
    screenPos = vec3(p, 2.5);
    vec3 ro = cameraPos;
    vec3 ta = cameraTarget;
    vec3 fo = normalize(ta - ro);
    vec3 ri = normalize(cross(vec3(0.,1.,0.), fo));
    vec3 up = normalize(cross(fo,ri));
    mat3 cam = mat3(ri,up,fo);
    vec3 ray = cam * normalize(screenPos);
    float depth = 0.01;
    vec2 mat;
    vec3 col = vec3(0.);
    
    vec3 pos = ro;
    sceneTrace(pos, ray, mat, depth, 100.);
    
    col = sceneShade(mat, pos, ray, depth, 100.);
    
    float sceneDepth = depth;
    depth = 0.01;
    pos = ro;
    screenTrace(pos, ray, mat, depth, sceneDepth);
    if (depth < sceneDepth) {
        col = screenShade(mat, pos);
    }
    
    float sceneAndScreenDepth = depth; // sceneDepth
    depth = 0.01;
    pos = ro;
    radiusTrace(pos, ray, mat, depth, sceneAndScreenDepth);
    if (depth < sceneAndScreenDepth) {
        vec4 radius = radiusShade(mat, pos);
        col = mix(col, radius.rgb, radius.a);
    }
    
    float sceneAndScreenAndGizmoDepth = sceneAndScreenDepth;
    depth = 0.01;
    pos = ro;
    gizmoTrace(pos, ray, mat, depth, sceneAndScreenAndGizmoDepth);
    if (depth < sceneAndScreenAndGizmoDepth) {
        vec4 gizmo = gizmoShade(mat, pos);
        col = mix(col, gizmo.rgb, gizmo.a);
    }
    return col;
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 uv = (iMouse.xy - iResolution.xy)/min(iResolution.x, iResolution.y);
    // timeLine(mod(length(uv.x), 41.));
    timeLine(mod(iTime/4., 4.5));

    vec2 p = (fragCoord * 2.0 - iResolution.xy)/min(iResolution.x, iResolution.y);
    vec2 dp = 1. / iResolution.xy;

    vec3 col = vec3(0.);
    
    
    // AA
    // https://www.shadertoy.com/view/Msl3Rr
    for(int y = 0; y < 2; y++) {
        for(int x = 0; x < 2; x++) {
            vec2 off = vec2(float(x),float(y))/2.;
            vec2 xy = (-iResolution.xy+2.0*(fragCoord + off * dp / 4.0)) / iResolution.y;
        	col += render(xy)*0.25;
        }
    }
    
    col = reinhard(col, 1.0, 1000.0);
    col = pow(col, vec3(1.0/2.2));
    
    //col = mix(col, vec3(mix(uv, p, uvToP), 0.), colToUv);
    
    fragColor = vec4(col,1.0);
}
                  adapted from kaneta's Raymarching in Raymarching
                  add the algorithm and the shader for sphere tracing

                  Theorem 8.  [gillespie2024ray, thm. 2] [ag-000U]

                    Given a function \(F: \mathbb {R} \rightarrow \mathbb {R}\) with Lipschitz bound \(\lambda \geq \) Lip \(F\), and an initial point \(t_0\), sphere tracing converges linearly to the smallest root greater than \(t_0\).

                    Remark 9.  [ag-000Q]

                      Theorem 8 means that sphere tracing requires the function to be Lipschitz, i.e. has a known Lipschitz bound.

                      Not all functions are Lipschitz, e.g. Harmonic Functions are not globally Lipschitz, many may have no local Lipschitz bound even within small neighborhoods, thus requires a new technique named Harnack tracing developed in [gillespie2024ray].

                      References

                        Reference. Ray tracing harmonic functions [gillespie2024ray]

                        • 2024
                        • Mark Gillespie, Denise Yang, Mario Botsch, Keenan Crane
                        @article{gillespie2024ray,
 title = {Ray Tracing Harmonic Functions},
 author = {Gillespie, Mark and Yang, Denise and Botsch, Mario and Crane, Keenan},
 year = {2024},
 journal = {ACM Trans. Graph},
 volume = {43},
 number = {4}
}

                        Reference. Ray-marching thurston geometries [coulon2022ray]

                        • 2022
                        • Rémi Coulon, Elisabetta A Matsumoto, Henry Segerman, Steve J Trettel
                        @article{coulon2022ray,
 title = {Ray-marching Thurston geometries},
 author = {Coulon, R{\'e}mi and Matsumoto, Elisabetta A and Segerman, Henry and Trettel, Steve J},
 year = {2022},
 journal = {Experimental Mathematics},
 volume = {31},
 number = {4},
 pages = {1197--1277},
 publisher = {Taylor \& Francis}
}

                        Reference. Sphere tracing: A geometric method for the antialiased ray tracing of implicit surfaces [hart1996sphere]

                        • 1996
                        • John C Hart
                        @article{hart1996sphere,
 title = {Sphere tracing: A geometric method for the antialiased ray tracing of implicit surfaces},
 author = {Hart, John C},
 year = {1996},
 journal = {The Visual Computer},
 volume = {12},
 number = {10},
 pages = {527--545},
 publisher = {Springer-Verlag Berlin/Heidelberg}
}

                        Context

                          Notes on ray-marching implicit surfaces [ag-000G]

                          Preface [ag-000H]

                            The note should contain:

                            • ray-marching ordinary SDFs
                            • ray-marching SDFs with singularities
                            • optimizations: bounded, adaptive etc.
                            • lighting
                            • shadow, including self-shadow
                            • edge detection
                            • auto-differentiation, curvature and geodesics
                            • non-Euclidean space
                            • dynamics, chaos, fractals
                            • volumetric ray marching
                            • meshes, polygonization, reconstruction
                            • higher dimensions
                            • PGA, CGA
                            • geometric construction based on topology and constraints
                            • neural implicit surfaces [Takikawa et al. 2021, Section 2]
                            • limitation

                            This is also an experiment to see how well the following 4 mix together in Forester:

                            • formal math notes of the concepts and formulas
                            • illustrative diagrams in TikZ
                            • psuedocode of the algorithms, then syntax highlighted acutal shade code
                            • the rendered GL canvas

                            See some links for rendering implicit surfaces, some links for learning shaders and 0910~0911 in Learning diary for some links as references. See test Three.js and shaders, test compute shaders for experiments on shaders.

                            Implicit surfaces [ag-000M]

                              Definition. Level set [wiki-level-set] [ag-000Y]

                                 A level set of a real-valued function \(f\) (a.k.a. an iso-contour of a scalar field in 3D) \[L_c(f)=\left \{\boldsymbol {x} \mid \forall \boldsymbol {x} \in X, f(\boldsymbol {x})=c \right \} \] where \(c\) is a constant.

                                Definition. Implicit surface [winchenbach2024lipschitz, sec. 1] [ag-000X]

                                  Let \(\Omega \) be a subset of a topological space \(X\), and \(\partial \Omega \) be its boundary.

                                  \(\partial \Omega \) as an implicit surface can be defined by a level set \(L_0(f)\) where \[ f(\boldsymbol {p})\begin {cases} \le 0 & \text { if } \boldsymbol {p} \in \Omega \\ > 0 & \text { if } \boldsymbol {p} \notin \Omega \end {cases} \]

                                  SDFs [ag-000Z]

                                    Definition. Signed distance function [wiki-sdf] [ag-000I]

                                      Let \(\Omega \) be a subset of a metric space \(X\) with metric \(d\), and \(\partial \Omega \) be its boundary. The distance between a point \(\boldsymbol {p}\) of \(X\) and the subset \(\partial \Omega \) of \(X\) is defined as usual as \[ d(\boldsymbol {p}, \partial \Omega )=\inf _{\boldsymbol {q} \in \partial \Omega } d(\boldsymbol {p}, \boldsymbol {q}) \] where \(\inf \) denotes the infimum, i.e. the greatest lower bound.

                                      The signed distance function or signed distance field (SDF) from a point \(\boldsymbol {p}\) of \(X\) to \(\Omega \) is defined by \[ f(\boldsymbol {p})= \begin {cases} -d(\boldsymbol {p}, \partial \Omega ) & \text { if } \boldsymbol {p} \in \Omega \\ d(\boldsymbol {p}, \partial \Omega ) & \text { if } \boldsymbol {p} \notin \Omega \end {cases} \]

                                      Convention. Sign [ag-000J]

                                        Simply put, SDFs are the minimum possible distance from a point to an implicit surface defined by \(f(\boldsymbol {p})=0\).

                                        The convention adopted in this note is that the \(+\) and \(-\) signs indicate whether the point is outside or inside the surface, respectively, so that when a ray marches towards the surface from the outside, the distance is positive, becomes smaller when approaching the surface.

                                        [wiki-sdf] uses a convention with opposite signs.

                                        Remark. SDFs in a Euclidean space [ag-0012]

                                          If \(X\) is also a normed space, i.e. for a point \(\boldsymbol {p} \in X\), its norm can be defined, e.g. in a Euclidean space \(\mathbb R^n\) where \(\boldsymbol {p}\) can be expressed by basis vectors as \((p_1, \ldots , p_n)\), then its norm \(\lVert \boldsymbol {p}\rVert \) can be defined as the Euclidean norm \[ \lVert \boldsymbol {p}\rVert _2:=\sqrt {p_1^2+\cdots +p_n^2} \] which is essentially the distance from the origin \(\boldsymbol {o}=(0, \ldots , 0)\) to the point \(\boldsymbol {p}\) \[ \lVert \boldsymbol {p}\rVert =d(\boldsymbol {p}, \boldsymbol {o})=\lVert \boldsymbol {p} - \boldsymbol {o}\rVert \] where \(d\) is the Euclidean distance function.

                                          In the following, we will be working with Euclidean space \(\mathbb R^n\), which is both a metric space and a normed space. Most of the time, we will focus on \(\mathbb R^3\) for visualization purposes.

                                          Example. SDF of a sphere [ag-000L]

                                            The SDF of a sphere with radius \(r\) centered at the origin \(\boldsymbol {o}\) is \[ f(\boldsymbol {p})=\lVert \boldsymbol {p}\rVert -r \] as depicted below:

                                            Figure. SDF for a sphere [ag-0017]

                                              The SDF above for a sphere translates to the following GLSL:

                                              float sdSphere(vec3 p, float r) { return length(p) - r; }

                                              Example. SDFs of geometric primitives [ag-0011]

                                                SDFs of some geometric primitives in \(\mathbb R^3\), can be found in Appendix A: Distance to Natural Quadrics and Torus of [hart1996sphere].

                                                Many more can be found in Inigo Quilez's article Distance functions (written in GLSL), which also includes an example that renders them by ray-marching.

                                                Figure. Raymarching primitives [ag-0014]

                                                   
                                                  #define HW_PERFORMANCE 0
// The MIT License
// Copyright © 2013 Inigo Quilez
// Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
//
// The license is here only not because I want to (can one
// license pieces of math?), but because people get upset
// if I don't add one...

// A list of useful distance function to simple primitives. All
// these functions (except for ellipsoid) return an exact
// euclidean distance, meaning they produce a better SDF than
// what you'd get if you were constructing them from boolean
// operations (such as cutting an infinite cylinder with two planes).

// List of other 3D SDFs:
//    https://www.shadertoy.com/playlist/43cXRl
// and
//    https://iquilezles.org/articles/distfunctions

#if HW_PERFORMANCE==0
#define AA 1
#else
#define AA 1   // make this 2 or 3 for antialiasing
#endif

//------------------------------------------------------------------
float dot2( in vec2 v ) { return dot(v,v); }
float dot2( in vec3 v ) { return dot(v,v); }
float ndot( in vec2 a, in vec2 b ) { return a.x*b.x - a.y*b.y; }

float sdPlane( vec3 p )
{
	return p.y;
}

float sdSphere( vec3 p, float s )
{
    return length(p)-s;
}

float sdBox( vec3 p, vec3 b )
{
    vec3 d = abs(p) - b;
    return min(max(d.x,max(d.y,d.z)),0.0) + length(max(d,0.0));
}

float sdBoxFrame( vec3 p, vec3 b, float e )
{
       p = abs(p  )-b;
  vec3 q = abs(p+e)-e;

  return min(min(
      length(max(vec3(p.x,q.y,q.z),0.0))+min(max(p.x,max(q.y,q.z)),0.0),
      length(max(vec3(q.x,p.y,q.z),0.0))+min(max(q.x,max(p.y,q.z)),0.0)),
      length(max(vec3(q.x,q.y,p.z),0.0))+min(max(q.x,max(q.y,p.z)),0.0));
}
float sdEllipsoid( in vec3 p, in vec3 r ) // approximated
{
    float k0 = length(p/r);
    float k1 = length(p/(r*r));
    return k0*(k0-1.0)/k1;
}

float sdTorus( vec3 p, vec2 t )
{
    return length( vec2(length(p.xz)-t.x,p.y) )-t.y;
}

float sdCappedTorus(in vec3 p, in vec2 sc, in float ra, in float rb)
{
    p.x = abs(p.x);
    float k = (sc.y*p.x>sc.x*p.y) ? dot(p.xy,sc) : length(p.xy);
    return sqrt( dot(p,p) + ra*ra - 2.0*ra*k ) - rb;
}

float sdHexPrism( vec3 p, vec2 h )
{
    vec3 q = abs(p);

    const vec3 k = vec3(-0.8660254, 0.5, 0.57735);
    p = abs(p);
    p.xy -= 2.0*min(dot(k.xy, p.xy), 0.0)*k.xy;
    vec2 d = vec2(
       length(p.xy - vec2(clamp(p.x, -k.z*h.x, k.z*h.x), h.x))*sign(p.y - h.x),
       p.z-h.y );
    return min(max(d.x,d.y),0.0) + length(max(d,0.0));
}

float sdOctogonPrism( in vec3 p, in float r, float h )
{
  const vec3 k = vec3(-0.9238795325,   // sqrt(2+sqrt(2))/2 
                       0.3826834323,   // sqrt(2-sqrt(2))/2
                       0.4142135623 ); // sqrt(2)-1 
  // reflections
  p = abs(p);
  p.xy -= 2.0*min(dot(vec2( k.x,k.y),p.xy),0.0)*vec2( k.x,k.y);
  p.xy -= 2.0*min(dot(vec2(-k.x,k.y),p.xy),0.0)*vec2(-k.x,k.y);
  // polygon side
  p.xy -= vec2(clamp(p.x, -k.z*r, k.z*r), r);
  vec2 d = vec2( length(p.xy)*sign(p.y), p.z-h );
  return min(max(d.x,d.y),0.0) + length(max(d,0.0));
}

float sdCapsule( vec3 p, vec3 a, vec3 b, float r )
{
	vec3 pa = p-a, ba = b-a;
	float h = clamp( dot(pa,ba)/dot(ba,ba), 0.0, 1.0 );
	return length( pa - ba*h ) - r;
}

float sdRoundCone( in vec3 p, in float r1, float r2, float h )
{
    vec2 q = vec2( length(p.xz), p.y );
    
    float b = (r1-r2)/h;
    float a = sqrt(1.0-b*b);
    float k = dot(q,vec2(-b,a));
    
    if( k < 0.0 ) return length(q) - r1;
    if( k > a*h ) return length(q-vec2(0.0,h)) - r2;
        
    return dot(q, vec2(a,b) ) - r1;
}

float sdRoundCone(vec3 p, vec3 a, vec3 b, float r1, float r2)
{
    // sampling independent computations (only depend on shape)
    vec3  ba = b - a;
    float l2 = dot(ba,ba);
    float rr = r1 - r2;
    float a2 = l2 - rr*rr;
    float il2 = 1.0/l2;
    
    // sampling dependant computations
    vec3 pa = p - a;
    float y = dot(pa,ba);
    float z = y - l2;
    float x2 = dot2( pa*l2 - ba*y );
    float y2 = y*y*l2;
    float z2 = z*z*l2;

    // single square root!
    float k = sign(rr)*rr*rr*x2;
    if( sign(z)*a2*z2 > k ) return  sqrt(x2 + z2)        *il2 - r2;
    if( sign(y)*a2*y2 < k ) return  sqrt(x2 + y2)        *il2 - r1;
                            return (sqrt(x2*a2*il2)+y*rr)*il2 - r1;
}

float sdTriPrism( vec3 p, vec2 h )
{
    const float k = sqrt(3.0);
    h.x *= 0.5*k;
    p.xy /= h.x;
    p.x = abs(p.x) - 1.0;
    p.y = p.y + 1.0/k;
    if( p.x+k*p.y>0.0 ) p.xy=vec2(p.x-k*p.y,-k*p.x-p.y)/2.0;
    p.x -= clamp( p.x, -2.0, 0.0 );
    float d1 = length(p.xy)*sign(-p.y)*h.x;
    float d2 = abs(p.z)-h.y;
    return length(max(vec2(d1,d2),0.0)) + min(max(d1,d2), 0.);
}

// vertical
float sdCylinder( vec3 p, vec2 h )
{
    vec2 d = abs(vec2(length(p.xz),p.y)) - h;
    return min(max(d.x,d.y),0.0) + length(max(d,0.0));
}

// arbitrary orientation
float sdCylinder(vec3 p, vec3 a, vec3 b, float r)
{
    vec3 pa = p - a;
    vec3 ba = b - a;
    float baba = dot(ba,ba);
    float paba = dot(pa,ba);

    float x = length(pa*baba-ba*paba) - r*baba;
    float y = abs(paba-baba*0.5)-baba*0.5;
    float x2 = x*x;
    float y2 = y*y*baba;
    float d = (max(x,y)<0.0)?-min(x2,y2):(((x>0.0)?x2:0.0)+((y>0.0)?y2:0.0));
    return sign(d)*sqrt(abs(d))/baba;
}

// vertical
float sdCone( in vec3 p, in vec2 c, float h )
{
    vec2 q = h*vec2(c.x,-c.y)/c.y;
    vec2 w = vec2( length(p.xz), p.y );
    
	vec2 a = w - q*clamp( dot(w,q)/dot(q,q), 0.0, 1.0 );
    vec2 b = w - q*vec2( clamp( w.x/q.x, 0.0, 1.0 ), 1.0 );
    float k = sign( q.y );
    float d = min(dot( a, a ),dot(b, b));
    float s = max( k*(w.x*q.y-w.y*q.x),k*(w.y-q.y)  );
	return sqrt(d)*sign(s);
}

float sdCappedCone( in vec3 p, in float h, in float r1, in float r2 )
{
    vec2 q = vec2( length(p.xz), p.y );
    
    vec2 k1 = vec2(r2,h);
    vec2 k2 = vec2(r2-r1,2.0*h);
    vec2 ca = vec2(q.x-min(q.x,(q.y < 0.0)?r1:r2), abs(q.y)-h);
    vec2 cb = q - k1 + k2*clamp( dot(k1-q,k2)/dot2(k2), 0.0, 1.0 );
    float s = (cb.x < 0.0 && ca.y < 0.0) ? -1.0 : 1.0;
    return s*sqrt( min(dot2(ca),dot2(cb)) );
}

float sdCappedCone(vec3 p, vec3 a, vec3 b, float ra, float rb)
{
    float rba  = rb-ra;
    float baba = dot(b-a,b-a);
    float papa = dot(p-a,p-a);
    float paba = dot(p-a,b-a)/baba;

    float x = sqrt( papa - paba*paba*baba );

    float cax = max(0.0,x-((paba<0.5)?ra:rb));
    float cay = abs(paba-0.5)-0.5;

    float k = rba*rba + baba;
    float f = clamp( (rba*(x-ra)+paba*baba)/k, 0.0, 1.0 );

    float cbx = x-ra - f*rba;
    float cby = paba - f;
    
    float s = (cbx < 0.0 && cay < 0.0) ? -1.0 : 1.0;
    
    return s*sqrt( min(cax*cax + cay*cay*baba,
                       cbx*cbx + cby*cby*baba) );
}

// c is the sin/cos of the desired cone angle
float sdSolidAngle(vec3 pos, vec2 c, float ra)
{
    vec2 p = vec2( length(pos.xz), pos.y );
    float l = length(p) - ra;
	float m = length(p - c*clamp(dot(p,c),0.0,ra) );
    return max(l,m*sign(c.y*p.x-c.x*p.y));
}

float sdOctahedron(vec3 p, float s)
{
    p = abs(p);
    float m = p.x + p.y + p.z - s;

    // exact distance
    #if 0
    vec3 o = min(3.0*p - m, 0.0);
    o = max(6.0*p - m*2.0 - o*3.0 + (o.x+o.y+o.z), 0.0);
    return length(p - s*o/(o.x+o.y+o.z));
    #endif
    
    // exact distance
    #if 1
 	vec3 q;
         if( 3.0*p.x < m ) q = p.xyz;
    else if( 3.0*p.y < m ) q = p.yzx;
    else if( 3.0*p.z < m ) q = p.zxy;
    else return m*0.57735027;
    float k = clamp(0.5*(q.z-q.y+s),0.0,s); 
    return length(vec3(q.x,q.y-s+k,q.z-k)); 
    #endif
    
    // bound, not exact
    #if 0
	return m*0.57735027;
    #endif
}

float sdPyramid( in vec3 p, in float h )
{
    float m2 = h*h + 0.25;
    
    // symmetry
    p.xz = abs(p.xz);
    p.xz = (p.z>p.x) ? p.zx : p.xz;
    p.xz -= 0.5;
	
    // project into face plane (2D)
    vec3 q = vec3( p.z, h*p.y - 0.5*p.x, h*p.x + 0.5*p.y);
   
    float s = max(-q.x,0.0);
    float t = clamp( (q.y-0.5*p.z)/(m2+0.25), 0.0, 1.0 );
    
    float a = m2*(q.x+s)*(q.x+s) + q.y*q.y;
	float b = m2*(q.x+0.5*t)*(q.x+0.5*t) + (q.y-m2*t)*(q.y-m2*t);
    
    float d2 = min(q.y,-q.x*m2-q.y*0.5) > 0.0 ? 0.0 : min(a,b);
    
    // recover 3D and scale, and add sign
    return sqrt( (d2+q.z*q.z)/m2 ) * sign(max(q.z,-p.y));;
}

// la,lb=semi axis, h=height, ra=corner
float sdRhombus(vec3 p, float la, float lb, float h, float ra)
{
    p = abs(p);
    vec2 b = vec2(la,lb);
    float f = clamp( (ndot(b,b-2.0*p.xz))/dot(b,b), -1.0, 1.0 );
	vec2 q = vec2(length(p.xz-0.5*b*vec2(1.0-f,1.0+f))*sign(p.x*b.y+p.z*b.x-b.x*b.y)-ra, p.y-h);
    return min(max(q.x,q.y),0.0) + length(max(q,0.0));
}

float sdHorseshoe( in vec3 p, in vec2 c, in float r, in float le, vec2 w )
{
    p.x = abs(p.x);
    float l = length(p.xy);
    p.xy = mat2(-c.x, c.y, 
              c.y, c.x)*p.xy;
    p.xy = vec2((p.y>0.0 || p.x>0.0)?p.x:l*sign(-c.x),
                (p.x>0.0)?p.y:l );
    p.xy = vec2(p.x,abs(p.y-r))-vec2(le,0.0);
    
    vec2 q = vec2(length(max(p.xy,0.0)) + min(0.0,max(p.x,p.y)),p.z);
    vec2 d = abs(q) - w;
    return min(max(d.x,d.y),0.0) + length(max(d,0.0));
}

float sdU( in vec3 p, in float r, in float le, vec2 w )
{
    p.x = (p.y>0.0) ? abs(p.x) : length(p.xy);
    p.x = abs(p.x-r);
    p.y = p.y - le;
    float k = max(p.x,p.y);
    vec2 q = vec2( (k<0.0) ? -k : length(max(p.xy,0.0)), abs(p.z) ) - w;
    return length(max(q,0.0)) + min(max(q.x,q.y),0.0);
}

//------------------------------------------------------------------

vec2 opU( vec2 d1, vec2 d2 )
{
	return (d1.x<d2.x) ? d1 : d2;
}

//------------------------------------------------------------------

#define ZERO (min(iFrame,0))

//------------------------------------------------------------------

vec2 map( in vec3 pos )
{
    vec2 res = vec2( pos.y, 0.0 );

    // bounding box
    if( sdBox( pos-vec3(-2.0,0.3,0.25),vec3(0.3,0.3,1.0) )<res.x )
    {
      res = opU( res, vec2( sdSphere(    pos-vec3(-2.0,0.25, 0.0), 0.25 ), 26.9 ) );
	  res = opU( res, vec2( sdRhombus(  (pos-vec3(-2.0,0.25, 1.0)).xzy, 0.15, 0.25, 0.04, 0.08 ),17.0 ) );
    }

    // bounding box
    if( sdBox( pos-vec3(0.0,0.3,-1.0),vec3(0.35,0.3,2.5) )<res.x )
    {
	res = opU( res, vec2( sdCappedTorus((pos-vec3( 0.0,0.30, 1.0))*vec3(1,-1,1), vec2(0.866025,-0.5), 0.25, 0.05), 25.0) );
    res = opU( res, vec2( sdBoxFrame(    pos-vec3( 0.0,0.25, 0.0), vec3(0.3,0.25,0.2), 0.025 ), 16.9 ) );
	res = opU( res, vec2( sdCone(        pos-vec3( 0.0,0.45,-1.0), vec2(0.6,0.8),0.45 ), 55.0 ) );
    res = opU( res, vec2( sdCappedCone(  pos-vec3( 0.0,0.25,-2.0), 0.25, 0.25, 0.1 ), 13.67 ) );
    res = opU( res, vec2( sdSolidAngle(  pos-vec3( 0.0,0.00,-3.0), vec2(3,4)/5.0, 0.4 ), 49.13 ) );
    }

    // bounding box
    if( sdBox( pos-vec3(1.0,0.3,-1.0),vec3(0.35,0.3,2.5) )<res.x )
    {
	res = opU( res, vec2( sdTorus(      (pos-vec3( 1.0,0.30, 1.0)).xzy, vec2(0.25,0.05) ), 7.1 ) );
    res = opU( res, vec2( sdBox(         pos-vec3( 1.0,0.25, 0.0), vec3(0.3,0.25,0.1) ), 3.0 ) );
    res = opU( res, vec2( sdCapsule(     pos-vec3( 1.0,0.00,-1.0),vec3(-0.1,0.1,-0.1), vec3(0.2,0.4,0.2), 0.1  ), 31.9 ) );
	res = opU( res, vec2( sdCylinder(    pos-vec3( 1.0,0.25,-2.0), vec2(0.15,0.25) ), 8.0 ) );
    res = opU( res, vec2( sdHexPrism(    pos-vec3( 1.0,0.2,-3.0), vec2(0.2,0.05) ), 18.4 ) );
    }

    // bounding box
    if( sdBox( pos-vec3(-1.0,0.35,-1.0),vec3(0.35,0.35,2.5))<res.x )
    {
	res = opU( res, vec2( sdPyramid(    pos-vec3(-1.0,-0.6,-3.0), 1.0 ), 13.56 ) );
	res = opU( res, vec2( sdOctahedron( pos-vec3(-1.0,0.15,-2.0), 0.35 ), 23.56 ) );
    res = opU( res, vec2( sdTriPrism(   pos-vec3(-1.0,0.15,-1.0), vec2(0.3,0.05) ),43.5 ) );
    res = opU( res, vec2( sdEllipsoid(  pos-vec3(-1.0,0.25, 0.0), vec3(0.2, 0.25, 0.05) ), 43.17 ) );
    res = opU( res, vec2( sdHorseshoe(  pos-vec3(-1.0,0.25, 1.0), vec2(cos(1.3),sin(1.3)), 0.2, 0.3, vec2(0.03,0.08) ), 11.5 ) );
    }

    // bounding box
    if( sdBox( pos-vec3(2.0,0.3,-1.0),vec3(0.35,0.3,2.5) )<res.x )
    {
    res = opU( res, vec2( sdOctogonPrism(pos-vec3( 2.0,0.2,-3.0), 0.2, 0.05), 51.8 ) );
    res = opU( res, vec2( sdCylinder(    pos-vec3( 2.0,0.14,-2.0), vec3(0.1,-0.1,0.0), vec3(-0.2,0.35,0.1), 0.08), 31.2 ) );
	res = opU( res, vec2( sdCappedCone(  pos-vec3( 2.0,0.09,-1.0), vec3(0.1,0.0,0.0), vec3(-0.2,0.40,0.1), 0.15, 0.05), 46.1 ) );
    res = opU( res, vec2( sdRoundCone(   pos-vec3( 2.0,0.15, 0.0), vec3(0.1,0.0,0.0), vec3(-0.1,0.35,0.1), 0.15, 0.05), 51.7 ) );
    res = opU( res, vec2( sdRoundCone(   pos-vec3( 2.0,0.20, 1.0), 0.2, 0.1, 0.3 ), 37.0 ) );
    }
    
    return res;
}

// https://iquilezles.org/articles/boxfunctions
vec2 iBox( in vec3 ro, in vec3 rd, in vec3 rad ) 
{
    vec3 m = 1.0/rd;
    vec3 n = m*ro;
    vec3 k = abs(m)*rad;
    vec3 t1 = -n - k;
    vec3 t2 = -n + k;
	return vec2( max( max( t1.x, t1.y ), t1.z ),
	             min( min( t2.x, t2.y ), t2.z ) );
}

vec2 raycast( in vec3 ro, in vec3 rd )
{
    vec2 res = vec2(-1.0,-1.0);

    float tmin = 1.0;
    float tmax = 20.0;

    // raytrace floor plane
    float tp1 = (0.0-ro.y)/rd.y;
    if( tp1>0.0 )
    {
        tmax = min( tmax, tp1 );
        res = vec2( tp1, 1.0 );
    }
    //else return res;
    
    // raymarch primitives   
    vec2 tb = iBox( ro-vec3(0.0,0.4,-0.5), rd, vec3(2.5,0.41,3.0) );
    if( tb.x<tb.y && tb.y>0.0 && tb.x<tmax)
    {
        //return vec2(tb.x,2.0);
        tmin = max(tb.x,tmin);
        tmax = min(tb.y,tmax);

        float t = tmin;
        for( int i=0; i<70 && t<tmax; i++ )
        {
            vec2 h = map( ro+rd*t );
            if( abs(h.x)<(0.0001*t) )
            { 
                res = vec2(t,h.y); 
                break;
            }
            t += h.x;
        }
    }
    
    return res;
}

// https://iquilezles.org/articles/rmshadows
float calcSoftshadow( in vec3 ro, in vec3 rd, in float mint, in float tmax )
{
    // bounding volume
    float tp = (0.8-ro.y)/rd.y; if( tp>0.0 ) tmax = min( tmax, tp );

    float res = 1.0;
    float t = mint;
    for( int i=ZERO; i<24; i++ )
    {
		float h = map( ro + rd*t ).x;
        float s = clamp(8.0*h/t,0.0,1.0);
        res = min( res, s );
        t += clamp( h, 0.01, 0.2 );
        if( res<0.004 || t>tmax ) break;
    }
    res = clamp( res, 0.0, 1.0 );
    return res*res*(3.0-2.0*res);
}

// https://iquilezles.org/articles/normalsSDF
vec3 calcNormal( in vec3 pos )
{
#if 0
    vec2 e = vec2(1.0,-1.0)*0.5773*0.0005;
    return normalize( e.xyy*map( pos + e.xyy ).x + 
					  e.yyx*map( pos + e.yyx ).x + 
					  e.yxy*map( pos + e.yxy ).x + 
					  e.xxx*map( pos + e.xxx ).x );
#else
    // inspired by tdhooper and klems - a way to prevent the compiler from inlining map() 4 times
    vec3 n = vec3(0.0);
    for( int i=ZERO; i<4; i++ )
    {
        vec3 e = 0.5773*(2.0*vec3((((i+3)>>1)&1),((i>>1)&1),(i&1))-1.0);
        n += e*map(pos+0.0005*e).x;
      //if( n.x+n.y+n.z>100.0 ) break;
    }
    return normalize(n);
#endif    
}

// https://iquilezles.org/articles/nvscene2008/rwwtt.pdf
float calcAO( in vec3 pos, in vec3 nor )
{
	float occ = 0.0;
    float sca = 1.0;
    for( int i=ZERO; i<5; i++ )
    {
        float h = 0.01 + 0.12*float(i)/4.0;
        float d = map( pos + h*nor ).x;
        occ += (h-d)*sca;
        sca *= 0.95;
        if( occ>0.35 ) break;
    }
    return clamp( 1.0 - 3.0*occ, 0.0, 1.0 ) * (0.5+0.5*nor.y);
}

// https://iquilezles.org/articles/checkerfiltering
float checkersGradBox( in vec2 p, in vec2 dpdx, in vec2 dpdy )
{
    // filter kernel
    vec2 w = abs(dpdx)+abs(dpdy) + 0.001;
    // analytical integral (box filter)
    vec2 i = 2.0*(abs(fract((p-0.5*w)*0.5)-0.5)-abs(fract((p+0.5*w)*0.5)-0.5))/w;
    // xor pattern
    return 0.5 - 0.5*i.x*i.y;                  
}

vec3 render( in vec3 ro, in vec3 rd, in vec3 rdx, in vec3 rdy )
{ 
    // background
    vec3 col = vec3(0.7, 0.7, 0.9) - max(rd.y,0.0)*0.3;
    
    // raycast scene
    vec2 res = raycast(ro,rd);
    float t = res.x;
	float m = res.y;
    if( m>-0.5 )
    {
        vec3 pos = ro + t*rd;
        vec3 nor = (m<1.5) ? vec3(0.0,1.0,0.0) : calcNormal( pos );
        vec3 ref = reflect( rd, nor );
        
        // material        
        col = 0.2 + 0.2*sin( m*2.0 + vec3(0.0,1.0,2.0) );
        float ks = 1.0;
        
        if( m<1.5 )
        {
            // project pixel footprint into the plane
            vec3 dpdx = ro.y*(rd/rd.y-rdx/rdx.y);
            vec3 dpdy = ro.y*(rd/rd.y-rdy/rdy.y);

            float f = checkersGradBox( 3.0*pos.xz, 3.0*dpdx.xz, 3.0*dpdy.xz );
            col = 0.15 + f*vec3(0.05);
            ks = 0.4;
        }

        // lighting
        float occ = calcAO( pos, nor );
        
		vec3 lin = vec3(0.0);

        // sun
        {
            vec3  lig = normalize( vec3(-0.5, 0.4, -0.6) );
            vec3  hal = normalize( lig-rd );
            float dif = clamp( dot( nor, lig ), 0.0, 1.0 );
          //if( dif>0.0001 )
        	      dif *= calcSoftshadow( pos, lig, 0.02, 2.5 );
			float spe = pow( clamp( dot( nor, hal ), 0.0, 1.0 ),16.0);
                  spe *= dif;
                  spe *= 0.04+0.96*pow(clamp(1.0-dot(hal,lig),0.0,1.0),5.0);
                //spe *= 0.04+0.96*pow(clamp(1.0-sqrt(0.5*(1.0-dot(rd,lig))),0.0,1.0),5.0);
            lin += col*2.20*dif*vec3(1.30,1.00,0.70);
            lin +=     5.00*spe*vec3(1.30,1.00,0.70)*ks;
        }
        // sky
        {
            float dif = sqrt(clamp( 0.5+0.5*nor.y, 0.0, 1.0 ));
                  dif *= occ;
            float spe = smoothstep( -0.2, 0.2, ref.y );
                  spe *= dif;
                  spe *= 0.04+0.96*pow(clamp(1.0+dot(nor,rd),0.0,1.0), 5.0 );
          //if( spe>0.001 )
                  spe *= calcSoftshadow( pos, ref, 0.02, 2.5 );
            lin += col*0.60*dif*vec3(0.40,0.60,1.15);
            lin +=     2.00*spe*vec3(0.40,0.60,1.30)*ks;
        }
        // back
        {
        	float dif = clamp( dot( nor, normalize(vec3(0.5,0.0,0.6))), 0.0, 1.0 )*clamp( 1.0-pos.y,0.0,1.0);
                  dif *= occ;
        	lin += col*0.55*dif*vec3(0.25,0.25,0.25);
        }
        // sss
        {
            float dif = pow(clamp(1.0+dot(nor,rd),0.0,1.0),2.0);
                  dif *= occ;
        	lin += col*0.25*dif*vec3(1.00,1.00,1.00);
        }
        
		col = lin;

        col = mix( col, vec3(0.7,0.7,0.9), 1.0-exp( -0.0001*t*t*t ) );
    }

	return vec3( clamp(col,0.0,1.0) );
}

mat3 setCamera( in vec3 ro, in vec3 ta, float cr )
{
	vec3 cw = normalize(ta-ro);
	vec3 cp = vec3(sin(cr), cos(cr),0.0);
	vec3 cu = normalize( cross(cw,cp) );
	vec3 cv =          ( cross(cu,cw) );
    return mat3( cu, cv, cw );
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 mo = iMouse.xy/iResolution.xy;
	float time = 32.0 + iTime*1.5;

    // camera	
    vec3 ta = vec3( 0.25, -0.75, -0.25 );
    vec3 ro = ta + vec3( 6.0*cos(0.1*time + 7.0*mo.x), 4.5, 6.0*sin(0.1*time + 7.0*mo.x) );
    // camera-to-world transformation
    mat3 ca = setCamera( ro, ta, 0.0 );

    vec3 tot = vec3(0.0);
#if AA>1
    for( int m=ZERO; m<AA; m++ )
    for( int n=ZERO; n<AA; n++ )
    {
        // pixel coordinates
        vec2 o = vec2(float(m),float(n)) / float(AA) - 0.5;
        vec2 p = (2.0*(fragCoord+o)-iResolution.xy)/iResolution.y;
#else    
        vec2 p = (2.0*fragCoord-iResolution.xy)/iResolution.y;
#endif

        // focal length
        const float fl = 2.5;
        
        // ray direction
        vec3 rd = ca * normalize( vec3(p,fl) );

         // ray differentials
        vec2 px = (2.0*(fragCoord+vec2(1.0,0.0))-iResolution.xy)/iResolution.y;
        vec2 py = (2.0*(fragCoord+vec2(0.0,1.0))-iResolution.xy)/iResolution.y;
        vec3 rdx = ca * normalize( vec3(px,fl) );
        vec3 rdy = ca * normalize( vec3(py,fl) );
        
        // render	
        vec3 col = render( ro, rd, rdx, rdy );

        // gain
        // col = col*3.0/(2.5+col);
        
		// gamma
        col = pow( col, vec3(0.4545) );

        tot += col;
#if AA>1
    }
    tot /= float(AA*AA);
#endif
    
    fragColor = vec4( tot, 1.0 );
}
                                                  adapted from Inigo Quilez's Raymarching - Primitives

                                                  Example. SDFs of exotic implicit surfaces [ag-0013]

                                                    SDFs of many exotic implicit surfaces can be found in:

                                                    The following example renders many exotic implicit surfaces by ray-marching:

                                                    Figure. Implicit surfaces [uts-000J]

                                                       
                                                      #define AS_LIB 1

int get_shape() {
    return int(iTime) % 52;
}

#include "/forest/shader/implicit.glsl"

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 uv = 2.*(fragCoord-iResolution.xy/2.)/iResolution.y; // contains [-1,1]^2
    vec3 col = vec3(0.);
    
    // Camera rays
    vec3 camPos = vec3(4.,0.,0.);
    vec3 camDir = - normalize(camPos);
    vec3 rayPos, rayDir;
    float zoom = 1.3; // 1.8*cos(iTime);
    // if (checkKey(KEY_E)) zoom = 0.5;
    float fov = 0.4*zoom;
    float fov_ortho = 1.5*zoom;
    #if perspective
        // perspective cam
        rayPos = camPos;
        rayDir = normalize(camDir + fov*vec3(0., uv.x, uv.y));
    #else
        // orthographic cam
        rayPos = camPos + fov_ortho*vec3(0., uv.x, uv.y);
        rayDir = camDir;
    #endif
    // for perspective background in orthographic mode
    vec3 cubemapDir = normalize(camDir + fov*vec3(0., uv.x, uv.y));
    
    // Mouse-controlled rotation
    vec2 mouse = initMouse + vec2(0.015625*sin(iTime*PI), 0.0); // initMouse; // iMouse.xy == vec2(0.,0.) ? initMouse : (iMouse.xy/iResolution.xy - 0.5);
    float yaw = clamp(- mouse.x * 2.*PI * 1., -PI,PI);
    float pitch = clamp( mouse.y * PI * 1.2, -PI*0.5, PI*0.5);
    // pitch and yaw rotations (column-wise matrices)
    mat3 rot = mat3(cos(yaw), sin(yaw), 0., -sin(yaw), cos(yaw), 0., 0., 0., 1.);
    rot = rot * mat3(cos(pitch), 0., -sin(pitch), 0., 1., 0., sin(pitch), 0., cos(pitch));
    // apply
    camPos = rot*camPos;
    camDir = rot*camDir;
    rayPos = rot*rayPos;
    rayDir = rot*rayDir;
    cubemapDir = rot*cubemapDir;
    
    //cubemapDir = vec3(cubemapDir.x, cubemapDir.z, cubemapDir.y);
    
    vec3 hitPoint = raycast(rayPos, rayDir);
    
    if (hitPoint == BINGO) { fragColor = vec4(BINGO,1.0); return; }
    //if (hitPoint == NOHIT) { fragColor = vec4(NOHIT,1.0); return; }
    //if (hitPoint == NOBOUNDHIT) { fragColor = vec4(NOBOUNDHIT,1.0); return; }
    //if (hitPoint == ESCAPEDBOUNDS) { fragColor = vec4(ESCAPEDBOUNDS,1.0); return; }
    //if (hitPoint == MAXDISTREACHED) { fragColor = vec4(MAXDISTREACHED,1.0); return; }
    //if (hitPoint == MAXITERREACHED) { fragColor = vec4(MAXITERREACHED,1.0); return; }
    
    if (hitPoint == NOBOUNDHIT  || hitPoint == NOHIT || hitPoint == ESCAPEDBOUNDS || hitPoint == MAXITERREACHED) {
        //fragColor = vec4(vec3(0.2),1.0); return;

        // make background transparent
        fragColor = vec4(0.0,0.0,0.0,0.0); return;
        
        col = with_background(cubemapDir);
        
        #if showBoundingCube
        // darken bounding cube
        if (hitPoint != NOBOUNDHIT) { col *= vec3(0.7); }
        #endif

        fragColor = vec4(col,1.0); return;
    }

    vec3 grad = gradf(hitPoint+1.1*EPS*(-rayDir));
    float s = -sign(dot(grad,rayDir));
    
    col = with_color_mode(grad, s, hitPoint, camPos);

    col = clamp(col, 0., 1.);

    col = with_surface_pattern(col, hitPoint);

    col = with_shading(col, grad, s, rayDir);
    
    col = clamp(col, 0., 1.);

    fragColor = vec4(col,1.0);
}
                                                      adapted from Rogolop's Implicit surfaces

                                                      Example. SDFs of implicit surfaces in exotic geometries [ag-0016]

                                                        Implicit surface could be defined in geometries other than Euclidean space, for some examples of such geometries, see [szirmay2022adapting] [coulon2022ray] [kopczynski2022real] and more.

                                                        The following example renders exotic implicit surfaces in exotic geometries by ray-marching:

                                                        Figure. Geodesic raytracing [ag-0015]

                                                           
                                                          //********************************************************************
//
// Begin of automatic differentiation header
// Full code with additional functions (gradients, jacobians, ...) can be found
// at: https://github.com/sibaku/glsl-autodiff
//
//********************************************************************
#ifndef HESSNUM_3_H_
#define HESSNUM_3_H_

// This file contains methods to compute the gradient and hessian 
// of a scalar valued 3 dimensional function using automatic forward differentiation

//--------------------------------
// Types
//--------------------------------

// Data type to hold information about a scalar valued 3 dimensional function
// These should be created by the constH3 (for constants) and varH3 (for variables) helpers
struct HNum3
{
    // The current value
    float val;
    // The current gradient
    vec3 g;
    // The current hessian
    mat3 h;
};

//--------------------------------
// Prototypes
//--------------------------------

/**
* Creates a constant HNum3
* @param val The current value of the constant
*/
HNum3 constH3(in float val);
/**
* Creates a HNum3 corresponding to the variable with the given index
* @param val The current value of the variable
* @param index The variable's index
*/
HNum3 varH3(in float val, in int index);
/**
* Creates a HNum3 corresponding to the variable x (index = 0)
* @param val The current value of the variable
*/
HNum3 varH3x(in float val);
/**
* Creates a HNum3 corresponding to the variable y (index = 1)
* @param val The current value of the variable
*/
HNum3 varH3y(in float val);
/**
* Creates a HNum3 corresponding to the variable z (index = 2)
* @param val The current value of the variable
*/
HNum3 varH3z(in float val);
HNum3 add(in HNum3 a, in HNum3 b);
HNum3 add(in HNum3 a, in float b);
HNum3 add(in float a, in HNum3 b);
HNum3 sub(in HNum3 a, in HNum3 b);
HNum3 sub(in HNum3 a, in float b);
HNum3 sub(in float a, in HNum3 b);
HNum3 mult(in HNum3 a, in HNum3 b);
HNum3 mult(in HNum3 a, in float b);
HNum3 mult(in float a, in HNum3 b);
HNum3 neg(in HNum3 a);
HNum3 div(in HNum3 a, in HNum3 b);
HNum3 div(in HNum3 a, in float b);
HNum3 div(in float a, in HNum3 b);
HNum3 inv(in HNum3 a);
HNum3 a_pow(in HNum3 a, in HNum3 b);
HNum3 a_pow(in HNum3 a, in float b);
HNum3 a_pow(in float a, in HNum3 b);
HNum3 a_min(in HNum3 a, in HNum3 b);
HNum3 a_max(in HNum3 a, in HNum3 b);
HNum3 a_exp2(in HNum3 a);
HNum3 a_inversesqrt(in HNum3 a);
HNum3 a_atan(in HNum3 a);
HNum3 a_sqrt(in HNum3 a);
HNum3 a_sinh(in HNum3 a);
HNum3 a_ceil(in HNum3 a);
HNum3 a_tan(in HNum3 a);
HNum3 a_asinh(in HNum3 a);
HNum3 a_asin(in HNum3 a);
HNum3 a_acosh(in HNum3 a);
HNum3 a_abs(in HNum3 a);
HNum3 a_exp(in HNum3 a);
HNum3 a_cosh(in HNum3 a);
HNum3 a_floor(in HNum3 a);
HNum3 a_log(in HNum3 a);
HNum3 a_atanh(in HNum3 a);
HNum3 a_log2(in HNum3 a);
HNum3 a_acos(in HNum3 a);
HNum3 a_tanh(in HNum3 a);
HNum3 a_cos(in HNum3 a);
HNum3 a_sin(in HNum3 a);
HNum3 a_atan2(in HNum3 y, in HNum3 x);
HNum3 a_atan2(in HNum3 y, in float x);
HNum3 a_atan2(in float y, in HNum3 x);
HNum3 a_mix(in HNum3 a, in HNum3 b, in HNum3 t);
HNum3 a_mix(in HNum3 a, in HNum3 b, in float t);
HNum3 a_mix(in HNum3 a, in float b, in HNum3 t);
HNum3 a_mix(in HNum3 a, in float b, in float t);
HNum3 a_mix(in float a, in HNum3 b, in HNum3 t);
HNum3 a_mix(in float a, in HNum3 b, in float t);
HNum3 a_mix(in float a, in float b, in HNum3 t);

//--------------------------------
// Macros
//--------------------------------

#define HESSIAN3(f,x, y, z,result)  {     result = f(varH3x(x), varH3y(y), varH3z(z)); }

//--------------------------------
// Utilities prototypes
//--------------------------------

mat3 a_outerProduct(in vec3 a, in vec3 b);

//--------------------------------
// Implementation
//--------------------------------

HNum3 constH3(in float val)
{
    return HNum3(val, vec3(0.0), mat3(0.0));
}
//--------------------------------
HNum3 varH3(in float val, in int index)
{   
    vec3 g = vec3(0.0);
    g[index] = 1.0;
    return HNum3(val, g, mat3(0.0));
}
//--------------------------------
HNum3 varH3x(in float val)
{   
    vec3 g = vec3(0.0);
    g[0] = 1.0;
    return HNum3(val, g, mat3(0.0));
}
//--------------------------------
HNum3 varH3y(in float val)
{   
    vec3 g = vec3(0.0);
    g[1] = 1.0;
    return HNum3(val, g, mat3(0.0));
}
//--------------------------------
HNum3 varH3z(in float val)
{   
    vec3 g = vec3(0.0);
    g[2] = 1.0;
    return HNum3(val, g, mat3(0.0));
}
//--------------------------------
HNum3 add(in HNum3 a, in HNum3 b)
{
    return HNum3(a.val + b.val , a.g + b.g, a.h + b.h);
}
//--------------------------------
HNum3 add(in HNum3 a, in float b)
{
    return HNum3(a.val + b , a.g, a.h);
}
//--------------------------------
HNum3 add(in float a, in HNum3 b)
{
    return HNum3(a + b.val , b.g, b.h);
}
//--------------------------------
HNum3 sub(in HNum3 a, in HNum3 b)
{
    return HNum3(a.val - b.val , a.g - b.g, a.h - b.h);
}
//--------------------------------
HNum3 sub(in HNum3 a, in float b)
{
    return HNum3(a.val - b , a.g, a.h);
}
//--------------------------------
HNum3 sub(in float a, in HNum3 b)
{
    return HNum3(a - b.val , - b.g, - b.h);
}
//--------------------------------
HNum3 mult(in HNum3 a, in HNum3 b)
{
    return HNum3(a.val * b.val, 
        a.val*b.g + b.val*a.g, 
        a.val*b.h + b.val*a.h + a_outerProduct(b.g,a.g) + a_outerProduct(a.g,b.g)
    );
}
//--------------------------------
HNum3 mult(in HNum3 a, in float b)
{
    return HNum3(a.val * b, b*a.g, b*a.h);
}
//--------------------------------
HNum3 mult(in float a, in HNum3 b)
{
    return HNum3(a * b.val, a*b.g, a*b.h);
}
//--------------------------------
HNum3 neg(in HNum3 a)
{
    return mult(-1.0,a);
}
//--------------------------------
HNum3 div(in HNum3 a, in HNum3 b)
{
    float b1 = b.val;
    float b2 = b1*b1;
    float b3 = b2*b1;

    return HNum3(a.val / b.val , 
        (b.val*a.g - a.val*b.g)/b2, 
        2.0*a.val/b3*a_outerProduct(b.g,b.g) 
        - a.val/b2*b.h
        + a.h/b1 
        - a_outerProduct(b.g/b2, a.g)
        - a_outerProduct(a.g/b2, b.g)
    );
}
//--------------------------------
HNum3 div(in HNum3 a, in float b)
{
    return HNum3(a.val / b, a.g/b, a.h/b);
}
//--------------------------------
HNum3 div(in float a, in HNum3 b)
{
    float b1 = b.val;
    float b2 = b1*b1;
    float b3 = b2*b1;

    return HNum3(a / b.val, 
        -a*b.g/b2, 
        2.0*a/b3*a_outerProduct(b.g,b.g) - a/b2*b.h
    );
}
//--------------------------------
HNum3 inv(in HNum3 a)
{
    return div(1.0, a);
}
//--------------------------------
HNum3 a_pow(in HNum3 a, in HNum3 b)
{
    return a_exp(mult(b,a_log(a)));
}
//--------------------------------
HNum3 a_pow(in HNum3 a, in float b)
{
    // constant exponent -> make special case
    float v = pow(a.val, b); // value f(a(x))
    float da = b*pow(a.val,b-1.0); // first derivative f'(a(x))
    float dda = b*(b-1.0)*pow(a.val,b-2.0); // second derivative f''(a(x))
    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_pow(in float a, in HNum3 b)
{
    return a_exp(mult(b,log(a)));
}
//--------------------------------
HNum3 a_min(in HNum3 a, in HNum3 b)
{
    if(a.val < b.val)
    {
        return a;
    }
    return b;
}
//--------------------------------
HNum3 a_max(in HNum3 a, in HNum3 b)
{
    if(a.val > b.val)
    {
        return a;
    }
    return b;
}
//--------------------------------
HNum3 a_exp2(in HNum3 a)
{
    float v = exp2(a.val); // value f(a(x))
    float da = log(2.0)*exp2(a.val); // first derivative f'(a(x))
    float dda = log(2.0)*log(2.0)*exp2(a.val); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_inversesqrt(in HNum3 a)
{
    float v = inversesqrt(a.val); // value f(a(x))
    float da = -0.5/pow(sqrt(a.val),3.0); // first derivative f'(a(x))
    float dda = 0.75/pow(sqrt(a.val),5.0); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_atan(in HNum3 a)
{
    float v = atan(a.val); // value f(a(x))
    float da = 1.0/(1.0 + a.val * a.val); // first derivative f'(a(x))
    float dda = -2.0*a.val/pow(1.0 + a.val * a.val, 2.0); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_sqrt(in HNum3 a)
{
    float v = sqrt(a.val); // value f(a(x))
    float da = 0.5/sqrt(a.val); // first derivative f'(a(x))
    float dda = -0.25/pow(sqrt(a.val),3.0); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_sinh(in HNum3 a)
{
    float v = sinh(a.val); // value f(a(x))
    float da = cosh(a.val); // first derivative f'(a(x))
    float dda = sinh(a.val); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_ceil(in HNum3 a)
{
    float v = ceil(a.val); // value f(a(x))
    float da = 0.0; // first derivative f'(a(x))
    float dda = 0.0; // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_tan(in HNum3 a)
{
    float v = tan(a.val); // value f(a(x))
    float da = 1.0 + pow(tan(a.val),2.0); // first derivative f'(a(x))
    float dda = 2.0*tan(a.val)*(1.0 + pow(tan(a.val),2.0)); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_asinh(in HNum3 a)
{
    float v = asinh(a.val); // value f(a(x))
    float da = 1.0/sqrt(1.0 + a.val * a.val); // first derivative f'(a(x))
    float dda = -a.val/pow(sqrt(1.0 + a.val * a.val),3.0); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_asin(in HNum3 a)
{
    float v = asin(a.val); // value f(a(x))
    float da = 1.0/sqrt(1.0 - a.val * a.val); // first derivative f'(a(x))
    float dda = a.val/pow(sqrt(1.0 - a.val * a.val),3.0); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_acosh(in HNum3 a)
{
    float v = acosh(a.val); // value f(a(x))
    float da = 1.0/sqrt(-1.0 + a.val * a.val); // first derivative f'(a(x))
    float dda = -a.val/pow(sqrt(-1.0 + a.val * a.val),3.0); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_abs(in HNum3 a)
{
    float v = abs(a.val); // value f(a(x))
    float da = a.val < 0.0 ? -1.0 : 1.0; // first derivative f'(a(x))
    float dda = 0.0; // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_exp(in HNum3 a)
{
    float v = exp(a.val); // value f(a(x))
    float da = exp(a.val); // first derivative f'(a(x))
    float dda = exp(a.val); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_cosh(in HNum3 a)
{
    float v = cosh(a.val); // value f(a(x))
    float da = sinh(a.val); // first derivative f'(a(x))
    float dda = cosh(a.val); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_floor(in HNum3 a)
{
    float v = floor(a.val); // value f(a(x))
    float da = 0.0; // first derivative f'(a(x))
    float dda = 0.0; // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_log(in HNum3 a)
{
    float v = log(a.val); // value f(a(x))
    float da = 1.0/a.val; // first derivative f'(a(x))
    float dda = -1.0/(a.val * a.val); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_atanh(in HNum3 a)
{
    float v = atanh(a.val); // value f(a(x))
    float da = 1.0/(1.0 - a.val * a.val); // first derivative f'(a(x))
    float dda = 2.0*a.val/pow(1.0 - a.val * a.val,2.0); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_log2(in HNum3 a)
{
    float v = log2(a.val); // value f(a(x))
    float da = 1.0/(a.val * log(2.0)); // first derivative f'(a(x))
    float dda = -1.0/(a.val * a.val * log(2.0)); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_acos(in HNum3 a)
{
    float v = acos(a.val); // value f(a(x))
    float da = -1.0/sqrt(1.0 - a.val * a.val); // first derivative f'(a(x))
    float dda = -a.val/pow(sqrt(1.0 - a.val * a.val),3.0); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_tanh(in HNum3 a)
{
    float v = tanh(a.val); // value f(a(x))
    float da = 1.0 - pow(tanh(a.val),2.0); // first derivative f'(a(x))
    float dda = -2.0*tanh(a.val)*(1.0 - pow(tanh(a.val),2.0)); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_cos(in HNum3 a)
{
    float v = cos(a.val); // value f(a(x))
    float da = -sin(a.val); // first derivative f'(a(x))
    float dda = -cos(a.val); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_sin(in HNum3 a)
{
    float v = sin(a.val); // value f(a(x))
    float da = cos(a.val); // first derivative f'(a(x))
    float dda = -sin(a.val); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_atan2(in HNum3 y, in HNum3 x)
{
    const float pi = 3.14159265; 
    // from https://en.wikipedia.org/wiki/Atan2
    if(x.val > 0.0)
    {
        HNum3 n = a_sqrt(add(mult(x,x),mult(y,y)));
        HNum3 inner = div(y, add(n,x));
        
        return mult(2.0,a_atan(inner));
        
    }else if(x.val <= 0.0 && abs(y.val) > 1E-6)
    {
        HNum3 n = a_sqrt(add(mult(x,x),mult(y,y)));
        HNum3 inner = div(sub(n,x),y);
         return mult(2.0,a_atan(inner));
    }else if(x.val < 0.0 && abs(y.val) <= 1E-6)
    {
        return constH3(pi);
    }
    // return 0 for undefined
    return constH3(0.0); 
}
//--------------------------------
HNum3 a_atan2(in HNum3 y, in float x)
{
    return a_atan2(y,constH3(x));
}
//--------------------------------
HNum3 a_atan2(in float y, in HNum3 x)
{
    return a_atan2(constH3(y),x);
}
//--------------------------------
HNum3 a_mix(in HNum3 a, in HNum3 b, in HNum3 t)
{
    return add(mult(a, sub(1.0, t)), mult(b, t));
}

//--------------------------------
HNum3 a_mix(in HNum3 a, in HNum3 b, in float t)
{
    return add(mult(a, 1.0 - t), mult(b, t));
}

//--------------------------------
HNum3 a_mix(in HNum3 a, in float b, in HNum3 t)
{
    return add(mult(a, sub(1.0, t)), mult(b, t));
}

//--------------------------------
HNum3 a_mix(in HNum3 a, in float b, in float t)
{
    return add(mult(a, 1.0 - t), b*t);
}

//--------------------------------
HNum3 a_mix(in float a, in HNum3 b, in HNum3 t)
{
    return add(mult(a, sub(1.0, t)), mult(b, t));
}

//--------------------------------
HNum3 a_mix(in float a, in HNum3 b, in float t)
{
    return add(a * (1.0 - t), mult(b, t));
}

//--------------------------------
HNum3 a_mix(in float a, in float b, in HNum3 t)
{
    return add(mult(a, sub(1.0, t)), mult(b, t));
}


//--------------------------------
// Implementation prototypes
//--------------------------------

mat3 a_outerProduct(in vec3 a, in vec3 b)
{
    return mat3(a * b[0], a * b[1], a * b[2]);
}


#endif // HESSNUM_3_H_
//********************************************************************
//
// End automatic differentiation header
//
//********************************************************************

//********************************************************************
//
// Metric functions
//
//********************************************************************
#define PI 3.14159265359

// float pSphere(vec3 p, float s) {return length(p)-s;}
HNum3 fSphere(in HNum3 x, in HNum3 y, in HNum3 z,
                 float globalTime) {             
                 
   return sub(a_sqrt(add(add(mult(x,x), mult(y,y)), mult(z,z))), 5.0)  ;
}

// Predefined functions!
// Just comment in the one you want to see! (comment out the other ones)
// If none is defined, the default identity operation is used (f(x,y,z) = 0).
// You can comment out all predefined ones and implement your own in fCustom

// #define EXPONENTIAL
// #define SADDLE// 
// #define WAVES
// #define PARABOLA
// #define IDENTITY

HNum3 fExponential(in HNum3 x, in HNum3 y, in HNum3 z,
                      float globalTime);
HNum3 fSaddle(in HNum3 x, in HNum3 y, in HNum3 z, float globalTime);
HNum3 fWaves(in HNum3 x, in HNum3 y, in HNum3 z, float globalTime);
HNum3 fParabola(in HNum3 x, in HNum3 y, in HNum3 z,
                   float globalTime);
HNum3 fIdentity(in HNum3 x, in HNum3 y, in HNum3 z,
                   float globalTime);
                   
// Implement your custom height function here and comment out all the defines
// above! You can use all the mathematical operations defined in the automatic
// differentiation header above
HNum3 fCustom(in HNum3 x, in HNum3 y, in HNum3 z,
                 float globalTime) {
  return div(fExponential(x, y, z, globalTime), 1.0);
  // return a_max(fExponential(x, y, z, globalTime), fSphere(x, y, z, globalTime));
}

// Height function used for the metric. More information in compute_update
HNum3 f(in HNum3 x, in HNum3 y, in HNum3 z, float globalTime) {

#if defined EXPONENTIAL
  return fExponential(x, y, z, globalTime);
#elif defined SADDLE
  return fSaddle(x, y, z, globalTime);
#elif defined WAVES
  return fWaves(x, y, z, globalTime);
#elif defined PARABOLA
  return fParabola(x, y, z, globalTime);
#elif defined IDENTITY
  return fIdentity(x, y, z, globalTime);
#else
  return fCustom(x, y, z, globalTime);
#endif
}

HNum3 hessF(vec3 p, float globalTime) {
  // compute hessian with addition time parameter

  vec3 uGrad = vec3(1., 0., 0.);
  HNum3 uHessian = HNum3(p.x, uGrad, mat3(0.));
  vec3 vGrad = vec3(0., 1., 0.);
  HNum3 vHessian = HNum3(p.y, vGrad, mat3(0.));
  vec3 wGrad = vec3(0., 0., 1.);
  HNum3 wHessian = HNum3(p.z, wGrad, mat3(0.));
  return f(uHessian, vHessian, wHessian, globalTime);
}
vec3 compute_update(vec3 p, vec3 v, float globalTime) {

  // explicit form of the christoffel symbols for the metric given by graph of
  // the function f(x_1,x_1,x_3) df/dx_i = f_i d^2f/(dx_j dx_j) = f_{i,j}
  // gamma_{i,j}^k = f_k*f_{i,j}/(1 + |grad(f)|^2)

  // This is the update given by the differential geodesic equation used for
  // numerical integration in local coordinates (x_i) for a curve c_i(t) = x_i
  // d^2 c_k/dt^2 = gamma_{i,j}^k * d c_i/dt * d c_j/dt

  // The local position's velocity is given by the ray's current direction

  // compute first and second order derivatives
  HNum3 r = hessF(p, globalTime);

  vec3 vp = vec3(0.0);

  for (int k = 0; k < 3; k++) {
    for (int i = 0; i < 3; i++) {
      for (int j = 0; j < 3; j++) {
        vp[k] += r.g[k] * r.h[j][i] * v[i] * v[j];
      }
    }
  }

  return -vp / (1.0 + dot(r.g, r.g));
}

//********************************************************************
//
// Data
//
//********************************************************************

const float deg2Rad = PI / 180.0;
const vec2 D_lastMouse = vec2(0., 0.);
vec4 D_thetaPhi;
const vec3 D_pos = vec3(0.0);

//********************************************************************
//
// Camera functions
//
//********************************************************************
// This assumes the pixel position px to be in [0,1],
// which can be done by (x+0.5)/w or (y+0.5)/h (or h-y +0.5 for screens
// with top left origin) to sample pixel centers
vec3 createRay(vec2 px, mat4 PInv, mat4 VInv);

mat4 translate(vec3 t);
mat4 translateInv(vec3 t);

mat4 scale(vec3 s);
mat4 scaleInv(vec3 s);

mat4 rightToLeft();
mat4 rightToLeftInv();

mat4 ortho(float l, float r, float b, float t, float n, float f);
mat4 orthoInv(float l, float r, float b, float t, float n, float f);

mat4 projection(float n, float f);
mat4 projectionInv(float n, float f);

mat4 perspective(float fov, float aspect, float n, float f);
mat4 perspectiveInv(float fov, float aspect, float n, float f);

mat4 lookAt(vec3 eye, vec3 center, vec3 up);
mat4 lookAtInv(vec3 eye, vec3 center, vec3 up);

//********************************************************************
//
// Implementation
//
//********************************************************************
HNum3 fExponential(in HNum3 x, in HNum3 y, in HNum3 z,
                      float globalTime) {
  HNum3 rx = sub(x, 0.0);
  HNum3 ry = sub(y, 1.0);
  HNum3 rz = sub(z, 2.0);

  HNum3 x2 = mult(rx, rx);
  HNum3 y2 = mult(ry, ry);
  HNum3 z2 = mult(rz, rz);

  HNum3 sum = neg(add(x2, add(y2, z2)));

  HNum3 f1 = a_exp(mult(sum, 1.0 / (abs(sin(globalTime / 5.0)) + 0.1)));

  rx = sub(x, -6.0);
  ry = sub(y, -0.25);
  rz = sub(z, 3.5);

  x2 = mult(rx, rx);
  y2 = mult(ry, ry);
  z2 = mult(rz, rz);

  sum = neg(add(x2, add(y2, z2)));
  sum = mult(sum, 0.1);
  HNum3 f2 = mult(a_exp(sum), abs(sin(globalTime / 5.0)) * 5.0);
  return add(f1, f2);
}
HNum3 fSaddle(in HNum3 x, in HNum3 y, in HNum3 z,
                 float globalTime) {
  return mult(0.25 * pow(sin(globalTime / 2.0), 2.0),
              sub(add(mult(x, x), mult(z, z)), mult(y, y)));
}
HNum3 fWaves(in HNum3 x, in HNum3 y, in HNum3 z, float globalTime) {
  float frequency = 2.0 * PI / (10.0 + 3.0 * sin(0.1 * globalTime));
  return a_cos(add(mult(x, frequency), mult(z, frequency)));
}
HNum3 fParabola(in HNum3 x, in HNum3 y, in HNum3 z,
                   float globalTime) {
  HNum3 rel = sub(x, 0.0);
  HNum3 sum = mult(rel, rel);

  rel = sub(y, 0.0);
  sum = add(sum, mult(rel, rel));

  rel = sub(z, -10.0);
  sum = add(sum, mult(rel, rel));
  return mult(sum, 0.1);
}
HNum3 fIdentity(in HNum3 x, in HNum3 y, in HNum3 z,
                   float globalTime) {
  return constH3(0.0);
}


//********************************************************************
vec3 createRay(vec2 px, mat4 PInv, mat4 VInv) {

  // convert pixel to NDS
  // [0,1] -> [-1,1]
  vec2 pxNDS = px * 2. - 1.;

  // choose an arbitrary point in the viewing volume
  // z = -1 equals a point on the near plane, i.e. the screen
  vec3 pointNDS = vec3(pxNDS, -1.);

  // as this is in homogenous space, add the last homogenous coordinate
  vec4 pointNDSH = vec4(pointNDS, 1.0);
  // transform by inverse projection to get the point in view space
  vec4 dirEye = PInv * pointNDSH;

  // since the camera is at the origin in view space by definition,
  // the current point is already the correct direction (dir(0,P) = P - 0 = P
  // as a direction, an infinite point, the homogenous component becomes 0
  // the scaling done by the w-division is not of interest, as the direction
  // in xyz will stay the same and we can just normalize it later
  dirEye.w = 0.;

  // compute world ray direction by multiplying the inverse view matrix
  vec3 dirWorld = (VInv * dirEye).xyz;

  // now normalize direction
  return normalize(dirWorld);
}

// matrix operations
mat4 translate(vec3 t) {
  return mat4(vec4(1., 0., 0., 0.), vec4(0., 1., 0., 0.), vec4(0., 0., 1., 0.),
              vec4(t, 1.));
}
mat4 translateInv(vec3 t) { return translate(-t); }

mat4 scale(vec3 s) {
  return mat4(vec4(s.x, 0., 0., 0.), vec4(0., s.y, 0., 0.),
              vec4(0., 0., s.z, 0.), vec4(0., 0., 0., 1.));
}
mat4 scaleInv(vec3 s) { return scale(1. / s); }

mat4 rightToLeft() {
  // 1 0 0  0
  // 0 1 0  0
  // 0 0 -1 0
  // 0 0 0  1
  return scale(vec3(1., 1., -1.));
}

mat4 rightToLeftInv() {
  // same matrix
  return rightToLeft();
}

mat4 ortho(float l, float r, float b, float t, float n, float f) {

  // translation and scale
  return scale(vec3(2. / (r - l), 2. / (t - b), 2. / (f - n))) *
         translate(vec3(-(l + r) / 2., -(t + b) / 2., -(f + n) / 2.));
}

mat4 orthoInv(float l, float r, float b, float t, float n, float f) {
  return translateInv(vec3(-(l + r) / 2., -(t + b) / 2., -(f + n) / 2.)) *
         scaleInv(vec3(2. / (r - l), 2. / (t - b), 2. / (f - n)));
}

mat4 projection(float n, float f) {
  // n 0 0 0	0
  // 0 n 0 0	0
  // 0 0 n+f	-fn
  // 0 0 1	0
  return mat4(vec4(n, 0., 0., 0.), vec4(0., n, 0., 0.), vec4(0., 0., n + f, 1.),
              vec4(0., 0., -f * n, 0.));
}

mat4 projectionInv(float n, float f) {
  // 1/n 	0 	0 		0
  // 0 	1/n	0 		0
  // 0	0	0 		1
  // 0	0	-1/fn	(f+n)/fn

  return mat4(vec4(1. / n, 0., 0., 0.), vec4(0., 1. / n, 0., 0.),
              vec4(0., 0., 0., -1. / (f * n)),
              vec4(0., 0., 1., (f + n) / (f * n)));
}

mat4 perspective(float fov, float aspect, float n, float f) {
  float l = tan(fov / 2.) * n;
  float b = l / aspect;

  return ortho(-l, l, -b, b, n, f) * projection(n, f) * rightToLeft();
}

mat4 perspectiveInv(float fov, float aspect, float n, float f) {
  float l = tan(fov / 2.) * n;
  float b = l / aspect;

  return rightToLeftInv() * projectionInv(n, f) * orthoInv(-l, l, -b, b, n, f);
}

mat4 lookAt(vec3 eye, vec3 center, vec3 up) {

  vec3 z = normalize(eye - center);
  vec3 x = normalize(cross(up, z));
  vec3 y = cross(z, x);

  mat4 v = mat4(vec4(x.x, y.x, z.x, 0.), vec4(x.y, y.y, z.y, 0.),
                vec4(x.z, y.z, z.z, 0.), vec4(0., 0., 0., 1.));

  return v * translate(-eye);
}

mat4 lookAtInv(vec3 eye, vec3 center, vec3 up) {
  vec3 z = normalize(eye - center);
  vec3 x = normalize(cross(up, z));
  vec3 y = cross(z, x);

  return translateInv(-eye) *
         mat4(vec4(x, 0.), vec4(y, 0.), vec4(z, 0.), vec4(0., 0., 0., 1.));
}

//********************************************************************
//
// INSTRUCTIONS:
//
// This shader implements basic distancefield raytracing of a regular grid on
// a 3D manifold given as the graph of a function f(x,y,z).
//
// Click and drag the mouse to look around. WASD for movement. You are moving
// through distorted space, so you might not move in a straight line, but
// instead move along the shortest paths in curved space.
// Press "R" to reset the view.
//
// WARNING: The space-bending and movement might be a bit nauseating. I tried to
// not make stuff move too fast, but wanted to give a short warning anyways.

// DISPLAY OPTIONS:
// Below you can find a few defines to control the look of the grid, background
// and general tracing
//
// FUNCTIONAL OPTIONS
// The actual function is defined in the common shader part under "Metric
// functions".  A few functions are predefined, but you can implement your own
// in "HessNum3 fCustom(in HessNum3 x, in HessNum3 y, in HessNum3 z, float
// globalTime)". Just comment in the one you want to see or comment all of them
// out to use the custom function. Additional information can be found there as
// well.
//
// Note: As calculating the metric from the function requires 1st and 2nd order
// derivatives and it would be cumbersome to calculate and implement them for
// each function, I use automatic differentiation, which is why you need to use
// the special types and operations defined in the automatic differentiation
// header below to implement your own f.
//
//
//********************************************************************

// Maximum number of trace steps
// Higher step number will let you look farther, especially with higher
// deformations
#define MAX_STEPS 100
// The maximum length traversed per step. Lower values produce higher quality
// paths through the curved space
#define STEP_LENGTH 0.2

// Grid parameters
#define GRID_SIZE 1.0
#define GRID_THICKNESS 0.0001

// Linear fog parameters
#define FOG_MAX_DIST 50.0
#define FOG_MIN_DIST 10.0
#define FOG_COLOR vec3(0.1)

// Surface epsilon
#define EPS 0.005

//********************************************************************
//
// Distance helpers and tracing
//
//********************************************************************
#define INF 3.402823466e+38

float sdf(in vec3 p);

vec3 g_eye;

// Tracing function
bool trace(inout vec3 p, in vec3 dir, out float t, out int steps) {

  steps = 0;
  t = 0.0;
  
  for (int n = 0; n < MAX_STEPS && t < FOG_MAX_DIST; n++) {

    float d = sdf(p);

    if (d < EPS) {
    
     if(length(p - g_eye) > 10.0 * GRID_SIZE) return false;
    
      return true;
    }


    // make step length dependent on distance field to not miss geometry, but
    // also use EPS as a min step
    float s = max(EPS, min(d, STEP_LENGTH));

    // "acceleration" of a particle traveling along the geodesic
    vec3 vp = compute_update(p, dir, iTime);

    // euler integration
    dir += vp * s;
    // position euler step
    p += dir * s;

    // maybe not fully correct, but try to preserve arc length by normalizing
    // dir
    dir = normalize(dir);

    // length traveled, used for fog
    t += s;

    steps++;
  }

  return false;
}

float sdfBox2d(vec2 p, vec2 b);
float sdfGrid(vec3 p, vec3 cellSize, vec3 gridHalfLength);
vec3 normal(in vec3 p);

//********************************************************************
//
// Image formation
//
//********************************************************************

vec3 shade(in vec3 P, in vec3 N, in vec3 color, in vec3 LPos);

float simpleColor(float coord, float gridSize, float repetition);

// vec4 loadValue(in ivec2 re) { return texelFetch(iChannel0, re, 0); }


void mainImage(out vec4 fragColor, in vec2 fragCoord) {

  vec2 uv = fragCoord.xy / iResolution.xy;

  float aspect = iResolution.x / iResolution.y;

  D_thetaPhi = vec4(90.0 * deg2Rad - clamp(2.0 * PI * iTime, 0., 2.0*PI) / 20., 0.0, 0.0, 0.0);

  vec3 rayDir;
  vec3 p;

  // {
    vec2 thetaPhi = D_thetaPhi.xy; // loadValue(D_thetaPhi).xy;
    float theta = thetaPhi.x;
    float phi = thetaPhi.y;

    vec3 eye = D_pos.xyz; // loadValue(D_pos).xyz;
    // vec3 eye = vec3(0.0);
    vec3 center =
        eye + vec3(sin(theta) * sin(phi), cos(theta), sin(theta) * cos(phi));
        
    g_eye = eye;
    // inverse projection and view matrices
    mat4 PInv = perspectiveInv(radians(90.), aspect, 0.1, 100.);
    mat4 VInv = lookAtInv(eye, center, vec3(0., 1., 0.));

    rayDir = createRay(uv, PInv, VInv);
    p = eye;
  // }

  float t = 0.0;
  int steps = 100;

  vec3 col = vec3(0.0);
  
  if (trace(p, rayDir, t, steps)) {
  
    vec3 N = normal(p);

    col =
        vec3(simpleColor(p.x, GRID_SIZE, 3.0), simpleColor(p.y, GRID_SIZE, 7.0),
             simpleColor(p.z, GRID_SIZE, 11.0));

    vec3 LPos = vec3(1.0, 2.0, 0.0);
    col = shade(p, N, col, LPos);

  } else {
    t = INF;
  }

  // simple linear fog
  float fogFactor = (FOG_MAX_DIST - t) / (FOG_MAX_DIST - FOG_MIN_DIST);
  fogFactor = clamp(fogFactor, 0.0, 1.0);
  col = mix(FOG_COLOR, col, fogFactor);

  fragColor = vec4(col, t == INF ? 0.0 : 1.0);
}

//***************

float sdfBox2d(vec2 p, vec2 b) {
  vec2 d = abs(p) - b;
  return min(max(d.x, d.y), 0.0) + length(max(d, 0.0));
}

float sdfGrid(vec3 p, vec3 cellSize, vec3 gridHalfLength) {
  vec3 pg = floor(p / cellSize) * cellSize;

  vec3 center = pg + cellSize * 0.5;

  float d = INF;
  
  d = min(d, sdfBox2d(p.xz - center.xz, gridHalfLength.xz));
  d = min(d, sdfBox2d(p.xy - center.xy, gridHalfLength.xy));
  d = min(d, sdfBox2d(p.yz - center.yz, gridHalfLength.yz));
  
  return d;
}

float pSphere(vec3 p, float s) {return length(p)-s;}

float sdf(in vec3 p) {
  float d = sdfGrid(p, vec3(GRID_SIZE), vec3(GRID_THICKNESS));
  
  // d = max(d, pSphere(p - g_center, 5.0 * GRID_SIZE));
  
  return d;
}

vec3 normal(in vec3 p) {
  return normalize(
      vec3(sdf(p + vec3(EPS, 0., 0.)) - sdf(p - vec3(EPS, 0., 0.)),
           sdf(p + vec3(0., EPS, 0.)) - sdf(p - vec3(0., EPS, 0.)),
           sdf(p + vec3(0., 0., EPS)) - sdf(p - vec3(0., 0., EPS))));
}

//********************************************************************
//
// Shading and main
//
//********************************************************************
vec3 shade(in vec3 P, in vec3 N, in vec3 color, in vec3 LPos) {
  vec3 L = normalize(LPos - P);
  // simple two sided diffuse shading
  return vec3(0.1, 0.8, 0.1) * clamp(abs(dot(L, N)), 0.0, 1.0);
}

float simpleColor(float coord, float gridSize, float repetition) {

  float cr = (gridSize * 0.5 - coord) / gridSize;

  float v = pow(sin(PI * cr / 2.0), 2.0);

  return v;
}
                                                          adapted from sibaku's Geodesic raytracer

                                                          Ray-marching [ag-000N]

                                                            Definition. Ray [ag-001A]

                                                              A ray is a half-line that starts at a point and extends indefinitely in one direction.

                                                              In the context of ray marching etc., a ray could be a view ray, i.e. a ray that starts at the camera, extending through the screen, possibly intersecting with objects in the scene, or a secondary ray, i.e. a ray generated from the intersection of other rays, e.g. shadow rays, reflection rays, or refraction rays.

                                                              Definition. Ray equation [gillespie2024ray, sec. 1.1 eq. 3] [ag-000V]

                                                                A ray anchored at origin \(\boldsymbol {r}_o\) in the direction of the unit vector \(\boldsymbol {r}_d\) can be parametrically defined as a ray equation \[ \boldsymbol {r}(t)=\boldsymbol {r}_o+t \boldsymbol {r}_d \]

                                                                Definition. Ray intersection [gillespie2024ray, sec. 1.1 eq. 4] [ag-000W]

                                                                  Plugging the ray equation \(r: \mathbb {R} \rightarrow \mathbb {R}^n\) into the function \(f: \mathbb {R}^n \rightarrow \mathbb {R}\) that defines the implicit surface produces the composite real function \(F: \mathbb {R} \rightarrow \mathbb {R}\) where \(F=f \circ \boldsymbol {r}\) such that the solutions to \[ F(t)=0 \] correspond to ray intersections with the implicit surface. Theses solutions are also called the roots of the function \(F\).

                                                                  Remark. Ray-casting/ray-tracing/ray-marching [ag-001E]

                                                                    The ray-casting/ray-tracing/ray-marching algorithms simply apply one of the multitude of root finding methods to solve the ray intersection equation.

                                                                    In general, ray-tracing typically uses a bunch of geometry (e.g. closed-form formulas) to calculate the intersection point, while ray-marching marches along the ray, and uses numerical methods to find the intersection point.

                                                                    The term ray-casting was introduced concurrently with ray-tracing in the 1980s, nowadays it could be used to mean the general ray intersection process, or specifically the ray-marching process with fixed step size.

                                                                    Definition. Ray marching [ag-001C]

                                                                      To render a scene, the Ray marching algorithm marches a ray from an origin towards a direction. At each step, the algorithm marches a short (possibly changing) distance and checks if a ray intersection occurs. This process continues until a stopping condition is met. The information obtained from the ray intersection is then used to determine the color of a pixel on the screen.

                                                                      Example. Ray-marching for different types of rays [ag-001D]

                                                                        A view ray has the origin at the camera, the direction is determined by the pixel on the screen that it passes through. It may intersect with a surface in the scene, and the information obtained from the intersection is used to determine the color of the corresponding pixel on the screen.

                                                                        A secondary ray typically originates from a point on a surface and its direction is determined by the type of ray. For instance, a shadow ray is cast towards a light source to determine if the point is in shadow, while a reflection ray and a refraction ray follow directions based on optics. The information obtained from secondary rays helps determine the color of the corresponding pixel that the view ray passes through.

                                                                        Figure. Ray marching [ag-001F]

                                                                          Inspired by Ray trace diagram from Wikipedia

                                                                          Definition. Ray marching (naïve) [ag-001B]

                                                                            The ray marching (naïve) algorithm is to march a ray by a fixed step size, check if a ray intersection occurs at each step, until the ray reaches a maximum distance or step count.

                                                                            Formally, the root found by ray marching (naïve) is the limit point of the sequence defined by the recurrence equation \[t_{i+1} = t_i + \Delta t\] where \(t_0 = 0\) and \(\Delta t\) is the fixed step size.

                                                                            At each step, the ray marches to \(\boldsymbol {r}(t_i)\), which is usually denoted \(\boldsymbol {r}_i\).

                                                                            Figure. Ray marching (naïve) [ag-000O]

                                                                              Inspired by Ray trace diagram from Wikipedia

                                                                              Algorithm. Ray marching (naïve) [ag-0018]

                                                                                Example. Ray marching (naïve) [ag-0019]

                                                                                  The following renders a scene with a unit sphere at the origin, the camera at \((0,0,8)\) and looking at the origin, through a screen of height 1.0, centered at 5.0 from the camera.

                                                                                  The color of pixels are adjusted so that they are brighter the closer they are to the camera, as if there is a light source at the camera. As it takes a fixed step size of 0.1 (not small enough), the sphere is not rendered smoothly, and clearly shows the stepping artifacts. 

                                                                                  /* The SDF for a sphere with radius r */
float sdSphere(vec3 p, float r)
{
  return length(p) - r;
}

/* The SDF for a scene with a unit sphere at the origin */
float sdScene(in vec3 pos) {
    return sdSphere(pos, 1.0);
}

#define EPSILON 0.01
#define T_MAX 10.0
#define DELTA_T 0.1
#define SDF_FUNCTION sdScene

/* Ray marching (naïve) */
float rayMarch(in vec3 ro, in vec3 rd) {

    float t = 0.0;

    while(t < T_MAX) {
        vec3 rt = ro + t * rd;
        float ft = SDF_FUNCTION(rt);
        
        if(ft < EPSILON) return t;

        t += DELTA_T;
    }
    
    return T_MAX;
}

/* The main function to render the scene per pixel */
void mainImage(out vec4 pixelColor, in vec2 pixelCoordinate) {
    
    // set up uv coordinates on the screen
    float aspect = iResolution.x/iResolution.y;
    vec2 uv = pixelCoordinate/iResolution.xy;
    uv -= vec2(0.5, 0.5);
    uv *= 2.0 * vec2(aspect, 1.0);

    // set up camera and view ray
    float observerDistance = 8.0;
    vec3 ro = vec3(0.0, 0.0, observerDistance);
    vec3 ta = ro - vec3(uv, 5.0);
    vec3 rd = normalize(ta - ro);

    // ray march
    float t = rayMarch(ro, rd);
    float hit = t != T_MAX ? 1.0 : 0.0;

    // adjust the color to look just bright enough
    float baseRatio = 0.65;
    // the closer, the brighter
    float distRatio = baseRatio + smoothstep(0.0, 1.0 - baseRatio, 1.0 - t / observerDistance);
    vec3 col = distRatio * vec3(1.0);

    // transparent if not hit
    pixelColor = vec4(col, hit);
}
                                                                                  simplified from zmcbeth's Raymarch 3.0
                                                                                  Code /* The SDF for a sphere with radius r */ float sdSphere(vec3 p, float r) { return length(p) - r; } /* The SDF for a scene with a unit sphere at the origin */ float sdScene(in vec3 pos) { return sdSphere(pos, 1.0); } #define EPSILON 0.01 #define T_MAX 10.0 #define DELTA_T 0.1 #define SDF_FUNCTION sdScene /* Ray marching (naïve) */ float rayMarch(in vec3 ro, in vec3 rd) { float t = 0.0; while(t < T_MAX) { vec3 rt = ro + t * rd; float ft = SDF_FUNCTION(rt); if(ft < EPSILON) return t; t += DELTA_T; } return T_MAX; } /* The main function to render the scene per pixel */ void mainImage(out vec4 pixelColor, in vec2 pixelCoordinate) { // set up uv coordinates on the screen float aspect = iResolution.x/iResolution.y; vec2 uv = pixelCoordinate/iResolution.xy; uv -= vec2(0.5, 0.5); uv *= 2.0 * vec2(aspect, 1.0); // set up camera and view ray float observerDistance = 8.0; vec3 ro = vec3(0.0, 0.0, observerDistance); vec3 ta = ro - vec3(uv, 5.0); vec3 rd = normalize(ta - ro); // ray march float t = rayMarch(ro, rd); float hit = t != T_MAX ? 1.0 : 0.0; // adjust the color to look just bright enough float baseRatio = 0.65; // the closer, the brighter float distRatio = baseRatio + smoothstep(0.0, 1.0 - baseRatio, 1.0 - t / observerDistance); vec3 col = distRatio * vec3(1.0); // transparent if not hit pixelColor = vec4(col, hit); }

                                                                                  Remark.  [ag-001G]

                                                                                    By convention, the GLSL community uses two-to-three letter variable names. Here are some common abbreviations:

                                                                                    • ro = ray origin
                                                                                    • rd = ray direction
                                                                                    • ca = camera
                                                                                    • ta = target (where the camera is pointing)
                                                                                    • uv = texture coordinates

                                                                                    Source: ⧉

                                                                                    Sphere tracing [ag-0010]
                                                                                    ✍️source

                                                                                      Definition. Point-to-set distance [hart1996sphere, def. 1] [ag-000R]

                                                                                        The point-to-set distance defines the distance from a point \(\boldsymbol {x} \in \mathbb {R}^3\) to a set \(A \subset \mathbb {R}^3\) as the distance from \(x\) to the closest point in \(A\), \[ d(\boldsymbol {x}, A)=\min _{y \in A}\lVert \boldsymbol {x}-\boldsymbol {y}\rVert \]

                                                                                        Definition. Signed distance bound [hart1996sphere, def. 2] [ag-000S]

                                                                                          A function \(f: \mathbb {R}^3 \rightarrow \mathbb {R}\) is a signed distance bound of its implicit surface \(f^{-1}(0)\) if and only if \[ |f(\boldsymbol {x})| \leq d\left (\boldsymbol {x}, f^{-1}(0)\right ) \] If equality holds, then \(f\) is a signed distance function.

                                                                                          Definition. Lipschitz bound [hart1996sphere, def. 3] [ag-000P]

                                                                                            A function \(f: \mathbb {R}^3 \rightarrow \mathbb {R}\) is Lipschitz over a domain \(D\) if and only if for all \(\boldsymbol {x}, \boldsymbol {y} \in D\), there exists a positive finite constant \(\lambda \), called Lipschitz bound, such that \[ |f(\boldsymbol {x})-f(\boldsymbol {y})| \leq \lambda \lVert \boldsymbol {x}-\boldsymbol {y}\rVert . \] The Lipschitz constant, denoted \(\operatorname {Lip} f\), is the minimum \(\lambda \) satisfying the condition above.

                                                                                            Theorem.  [gillespie2024ray, thm. 1] [ag-000T]

                                                                                              Let \(f\) be Lipschitz with Lipschitz bound \(\lambda \geq \) Lip \(f\). Then the function \(f / \lambda \) is a signed distance bound of its implicit surface.

                                                                                              Definition. Sphere tracing [hart1996sphere, sec. 2.3, eq. 12] [ag-001I]

                                                                                                The sphere tracing algorithm is to march a ray by an adaptive safe step size, which is the absolute value of the signed distance bound calculated by \(F/\lambda \) per Theorem [ag-000T], the rest is the same as ray marching (naïve).

                                                                                                Formally, the root found by sphere tracing is the limit point of the sequence defined by the recurrence equation \[t_{i+1} = t_i + \frac {|F(t_i)|}{\lambda } = t_i + \frac {f(\boldsymbol {r}(t_i))}{\lambda } \] where \(F\) is the ray intersection, and \(f\) is the SDF.

                                                                                                Usually \(\lambda = 1\), and \(F(t_i)\) stays positive until the stopping condition \(F(t_i) < \epsilon \) is met, so the above can be simplified to \[t_{i+1} = t_i + F(t_i) = t_i + f(\boldsymbol {r}(t_i)) \]

                                                                                                Figure. Sphere tracing [ag-001H]

                                                                                                  Inspired by [coulon2022ray, fig. 2.2]

                                                                                                  Figure. Raymarching in raymarching (sphere tracing) [ag-001J]

                                                                                                    #define MAT_SCREEN 0.
#define MAT_SPHERE1 1.
#define MAT_FLOOR  2.
#define MAT_CORNER 3.
#define MAT_SCREENZ 4.
#define MAT_MARCHSPHERE 5.
#define MAT_MARCHROUTE 6.
#define MAT_MARCHRADIUS 7.
#define MAT_SPHERE2 8.
#define MAT_SPHERE3 9.
#define MAT_RAYDIRECTION 10.
#define MAT_HITPOINT 11.

float uvToP = 0.0;
float colToUv = 1.0;
float screenZ = 2.5;

float sphereAnim = 1.0;
float radiusAnim = 1.0;
float routeAnim = 1.0;

float raySphereAlpha = 0.0;

float cornersAlpha = 0.0;
float cornersAnim = 0.0;

float screenZAlpha = 0.0;

float radiusAlpha = 1.0;

float rayDirectionAnim = 0.0;
float rayDirectionAnim2 = 0.0;

float screenAlpha = 0.0;

float hitSphereID = 0.0;

// =====================Camera========================
vec3 cameraPos, cameraTarget;

//================================

float sum = 0.0;

float tl(float val, float offset, float range)
{
    float im = sum + offset;
    float ix = im + range;
    sum += offset + range;
    return clamp((val - im) / (ix - im), 0.0, 1.0);
}

float cio(float t) {
	return t < 0.5
	? 0.5 * (1.0 - sqrt(1.0 - 4.0 * t * t))
	: 0.5 * (sqrt((3.0 - 2.0 * t) * (2.0 * t - 1.0)) + 1.0);
}

float eio(float t) {
	return t == 0.0 || t == 1.0
		? t
		: t < 0.5
		? +0.5 * pow(2.0, (20.0 * t) - 10.0)
		: -0.5 * pow(2.0, 10.0 - (t * 20.0)) + 1.0;
}

float fio(float t) {
	return t == 0.0 || t == 1.0
		? t
		: t <= 0.5
		? -0.5 * pow(2.0, 5.0 - ((t + 0.5) * 10.0)) + 1.0 - 0.5
		: +0.5 * pow(2.0, (10.0 * (t - 0.5)) - 5.0) + 0.5;
}

void timeLine(float time)
{
    cameraPos = vec3(5., 5., -3.); // mix(vec3(0.0), , eio(t));
    cameraTarget = vec3(0.0, 0.0, 5.);
    
    //time += 32.;
    float t; // = tl(time, 1.0, 0.5);
    // uvToP = mix(0.0, 1.0, eio(t));
    
    // t = tl(time, 1.0, 1.0);
    
    // t = tl(time, 0.5, 1.0);
    // raySphereAlpha = mix(0., 1., t);
    // cornersAlpha = mix(0., 1., t);
    
    // t = tl(time, 0.5, 1.0);
    // cornersAnim = mix(0., 1., t);
    
    t = tl(time, 0.0, 0.5);
    // cameraPos = mix(cameraPos, vec3(5., 6., -3.) , t); // eio(t));
    ///cameraTarget = vec3(0.0, 0.0, 3.0); //mix(cameraTarget, , eio(t));
    
    // t = tl(time, 0.5, 0.5);
    // screenZAlpha = mix(0., 1., t);

	// t = tl(time, 0.5, .5);
    // colToUv = mix(colToUv, 0.0, eio(t));
    
    // t = tl(time, 0.5, 1.0);
    // screenZ = mix(screenZ, 5.0, eio(t));
    
    // t = tl(time, 1.0, 1.0);
    // screenZ = mix(screenZ, .75, eio(t));
    
    
    // t = tl(time, 1.0, 1.0);
    // screenZ = mix(screenZ, 2.5, eio(t));
    
    
    // t = tl(time, 0.5, 1.0);
    // cornersAlpha = mix(cornersAlpha, 0., t);
    // screenZAlpha = mix(screenZAlpha, 0., t);
    // colToUv = mix(colToUv, 1.0, eio(t));
    
    // t = tl(time, 0.0, 0.0);
    // cornersAnim = mix(cornersAnim, 0., t);
    
    // t = tl(time, 0.5, 1.0);
    // cameraPos = mix(cameraPos, vec3(5., 15., 6.0), eio(t));
    // cameraTarget = mix(cameraTarget, vec3(0.0, 0.0, 6.), eio(t));
    
    raySphereAlpha = 1.;
    
    for (int i=0; i<20; i++) {
        t = tl(time, 0., 0.5);
        radiusAnim = mix(radiusAnim, float(i) + 2., t);
        t = tl(time, 0., 0.5);
        routeAnim = mix(routeAnim, float(i) + 2., t);
        t = tl(time, 0., 0.5);
        sphereAnim = mix(sphereAnim, float(i) + 2., t);
    }
    
    
    // t = tl(time, 1.0, 1.0);
    // cameraPos = mix(cameraPos, vec3(2., 3., -3.), eio(t));
    // cameraTarget = mix(cameraTarget, vec3(0.0, 0.0, 3.0), eio(t));
    // radiusAlpha = mix(radiusAlpha, 0.0, t);
    // routeAnim = mix(routeAnim, 1.0, t);
    // sphereAnim = mix(sphereAnim, 1.0, t);
    // colToUv = mix(colToUv, 1.0, eio(t));
    // screenAlpha = mix(screenAlpha, 0.0, t);
    
    t = tl(time, 0.5, 0.0);
    radiusAnim = mix(radiusAnim, 0.0, t);
    routeAnim = mix(routeAnim, 0.0, t);
    sphereAnim = mix(sphereAnim, 0.0, t);

    
    t = tl(time, 0.5, 1.0);
    rayDirectionAnim2 = mix(rayDirectionAnim2, 1.0, t); // fio(t)); // eio(t));
    
    t = tl(time, 0.5, 3.0);
    rayDirectionAnim = mix(rayDirectionAnim, 1.0, fio(t));
    
    t = tl(time, 3.5, 0.0);
    rayDirectionAnim2 = mix(rayDirectionAnim2, 0.0, t);
    raySphereAlpha = mix(raySphereAlpha, 0.0, t);
    
    // rayDirectionAnim = mix(rayDirectionAnim, 0.0, t);
    // raySphereAlpha = mix(raySphereAlpha, 0.0, t);
    // rayDirectionAnim2 = mix(rayDirectionAnim2, 0.0, t);
    
    
    // t = tl(time, 0.5, 1.0);
    // cameraPos = mix(cameraPos, vec3(0.), eio(t));
    // cameraTarget = mix(cameraTarget, vec3(0.0, 0.0, 5.), eio(t));
}

vec2 opU(vec2 a, vec2 b)
{
    return a.x < b.x ? a : b;
}

float sdLine( vec3 p, vec3 a, vec3 b, float r )
{
    vec3 pa = p - a, ba = b - a;
    float h = clamp( dot(pa,ba)/dot(ba,ba), 0.0, 1.0 );
    return length( pa - ba*h ) - r;
}

float sdBox( vec3 p, vec3 b )
{
  vec3 d = abs(p) - b;
  return length(max(d,0.0))
         + min(max(d.x,max(d.y,d.z)),0.0); // remove this line for an only partially signed sdf 
}

float sdSphere(vec3 p, float s)
{
    return length(p) - s;
}

vec3 sunDir = normalize(vec3(.3, .25, .2));
vec3 skyColor(vec3 rd)
{
    float sundot = clamp(dot(rd,sunDir),0.0,1.0);
    // sky		
    vec3 col = mix(vec3(0.2,0.5,0.85)*1.1, vec3(0.0,0.15,0.7), rd.y);
    col = mix( col, 0.85*vec3(0.8,0.8,0.7), pow( 1.0-max(rd.y,0.0), 4.0 ) );
    // sun
    col += 0.3*vec3(1.0,0.7,0.4)*pow( sundot,5.0 );
    col += 0.5*vec3(1.0,0.8,0.6)*pow( sundot,64.0 );
    col += 6.0*vec3(1.0,0.8,0.6)*pow( sundot,1024.0 );
    return col * 3.0;
}

vec3 screenPos;

vec2 sceneSpheres(vec3 p)
{
    vec2 s1 = vec2(sdSphere(p - vec3(-2.0, 0.0, 6.0), 1.), MAT_SPHERE1);
    vec2 s2 = vec2(sdSphere(p - vec3(0.0, 0.0, 16.0), 1.), MAT_SPHERE2);
    vec2 s3 = vec2(sdSphere(p - vec3(2.0, 0.0, 9.0), 1.), MAT_SPHERE3);
    return opU(opU(s1, s2), s3);
}

vec2 sceneMap(vec3 p)
{
    vec2 s = sceneSpheres(p);
    vec2 p1 = vec2(p.y + 2.0, MAT_FLOOR);
    return opU(s, p1);
}

vec3 normal(vec3 p)
{
    vec2 e = vec2(0.001, 0.0);
    float d = sceneMap(p).x;
    return normalize(d - vec3(sceneMap(p - e.xyy).x, sceneMap(p - e.yxy).x, sceneMap(p - e.yyx).x));
}

void sceneTrace(inout vec3 pos, vec3 ray, out vec2 mat, inout float depth, float maxD)
{
    vec3 ro = pos;
    int i = 0;
    for(; i < 70; i++) {
        if (depth > maxD) {
        	depth = maxD;
            break;
        }
    	pos = ro + ray * depth;
        mat = sceneMap(pos);
        if (mat.x < 0.001) {
            break;
        }
        depth += mat.x;
    }
}

vec2 screenMap(vec3 p)
{
    vec2 screenSize = iResolution.xy / min(iResolution.x, iResolution.y);
    vec2 b1 = vec2(sdBox(p - vec3(0., 0., screenZ), vec3(screenSize, 0.01)), MAT_SCREEN);
    return b1;
}

void screenTrace(inout vec3 pos, vec3 ray, out vec2 mat, inout float depth, float maxD)
{
    vec3 ro = pos;
    for(int i = 0; i < 20; i++) {
        if (depth > maxD) {
        	depth = maxD;
            break;
        }
    	pos = ro + ray * depth;
        mat = screenMap(pos);
        if (mat.x < 0.001) {
            break;
        }
        depth += mat.x;
    }
}

float remap(float val, float im, float ix, float om, float ox)
{
    return clamp(om + (val - im) * (ox - om) / (ix - im), om, ox);
}

vec2 gizmoCorners(vec3 p)
{
    vec2 screenSize = iResolution.xy / min(iResolution.x, iResolution.y);
    
    float a1 = remap(cornersAnim, 0.0, 0.25, 0.0, 1.0);
    float a2 = remap(cornersAnim, 0.25, 0.5, 0.0, 1.0);
    float a3 = remap(cornersAnim, 0.5, 0.75, 0.0, 1.0);
    float a4 = remap(cornersAnim, 0.75, 1.0, 0.0, 1.0);
    
    vec2 c1 = vec2(sdLine(p, vec3(0.), mix(vec3(0.), vec3(screenSize, screenZ), a1), 0.02), MAT_CORNER);
    vec2 c2 = vec2(sdLine(p, vec3(0.), mix(vec3(0.), vec3(screenSize * vec2(1.,-1.), screenZ), a2), 0.02), MAT_CORNER);
    vec2 c3 = vec2(sdLine(p, vec3(0.), mix(vec3(0.), vec3(screenSize * vec2(-1.,-1.), screenZ), a3), 0.02), MAT_CORNER);
    vec2 c4 = vec2(sdLine(p, vec3(0.), mix(vec3(0.), vec3(screenSize * vec2(-1.,1.), screenZ), a4), 0.02), MAT_CORNER);
    return opU(c1, opU(c2, opU(c3, c4)));
}

vec2 gizmoScreenZ(vec3 p)
{
    vec2 c1 = vec2(sdLine(p, vec3(0.), vec3(0., 0., screenZ), 0.03), MAT_SCREENZ);
    return c1;
}

float sphereID = 0.0;

vec3 rayRoute = vec3(0., 1., 0.); // y must be init to non-zero
vec2 gizmoMarching(vec3 p)
{
    vec3 ray = vec3(0., 0., 1.);
    vec2 d = vec2(10000.);
    
    if(rayDirectionAnim2 > 0.0) {
        if(abs(rayRoute.y - 0.) <= 0.01) {
            ray = normalize(vec3(rayRoute.x, rayRoute.y, rayRoute.z));
            radiusAnim = 20.;
            routeAnim = 20.;
            sphereAnim = 20.;
        }
    }
    
    float t = 0.0;
    vec3 pos;
    int maxSteps = hitSphereID == 0. ? 20 : int(hitSphereID) + 1;
    for(int i=0; i<maxSteps; i++) {
    	pos = ray * t;
        vec2 s = vec2(sdSphere(p - pos, 0.15), MAT_MARCHSPHERE);
        if (s.x < d.x) {
            d = s;
            sphereID = float(i);   
        }
        
        float dist = sceneSpheres(pos).x;
        
        float anim = clamp(routeAnim - float(i) - 1., 0.0, 1.0);
        vec2 c1 = vec2(sdLine(p, pos, pos + mix(vec3(0.), ray * dist, anim), 0.03), MAT_MARCHROUTE);
        
        d = opU(d, c1);
        t += dist;
        
    }    
    return d;
}

vec2 gizmoMarchingRadius(vec3 p)
{
    vec2 d = vec2(p.y, MAT_MARCHRADIUS);
    return d;
}


vec2 gizmoRayDirection(vec3 p)
{
    float a1 = fract(rayDirectionAnim * 19.99999);
    float a2 = floor(rayDirectionAnim * 19.99999) / 20.;
    vec2 screenSize = iResolution.xy / min(iResolution.x, iResolution.y);
    screenSize.y *= -1.;
    screenSize = mix(-screenSize, screenSize, vec2(a1, a2));
    rayRoute = mix(vec3(0.), vec3(screenSize, screenZ) * 10.0, rayDirectionAnim2);
    vec2 c1 = vec2(sdLine(p, vec3(0.), mix(vec3(0.), vec3(screenSize, screenZ) * 10.0, rayDirectionAnim2), 0.02), MAT_RAYDIRECTION);
    
    vec3 ray = normalize(vec3(screenSize, screenZ));
    vec3 pos;
    float t = 0.01;
    for(int i = 0; i < 40; i++) {
    	pos = ray * t;
        vec2 d = sceneMap(pos);
        if (d.x < 0.001) {
            break;
        }
        t += d.x;
    }
    
    c1 = opU(c1, vec2(sdSphere(p - pos, 0.15), MAT_HITPOINT));
    
    return c1;
}

vec2 gizmoMap(vec3 p)
{
    vec2 d = opU(gizmoCorners(p), gizmoScreenZ(p));
    d = opU(d, gizmoMarching(p));
    d = opU(d, gizmoRayDirection(p));
    return d;
}

void gizmoTrace(inout vec3 pos, vec3 ray, out vec2 mat, inout float depth, float maxD)
{
    vec3 ro = pos;
    for(int i = 0; i < 60; i++) {
        if (depth > maxD) {
        	depth = maxD;
            break;
        }
    	pos = ro + ray * depth;
        mat = gizmoMap(pos);
        if (mat.x < 0.001) {
            break;
        }
        depth += mat.x;
    }
}

vec2 radiusMap(vec3 p)
{
    return gizmoMarchingRadius(p);
}

void radiusTrace(inout vec3 pos, vec3 ray, out vec2 mat, inout float depth, float maxD)
{
    vec3 ro = pos;
    for(int i = 0; i < 40; i++) {
        if (depth > maxD) {
        	depth = maxD;
            break;
        }
    	pos = ro + ray * depth;
        mat = radiusMap(pos);
        if (mat.x < 0.001 || depth > maxD) {
            break;
        }
        depth += mat.x;
    }
}

float shadow(in vec3 p, in vec3 l, float ma)
{
    float t = 0.1;
    float t_max = ma;
    
    float res = 1.0;
    for (int i = 0; i < 16; ++i)
    {
        if (t > t_max) break;
        vec3 pos = p + t*l;
        float d = opU(sceneMap(pos), screenMap(pos)).x; // sceneMap(pos).x; //opU(, screenMap(pos)).x;
        if (d < 0.001)
        {
            return 0.0;
        }
        t += d*1.0;
        res = min(res, 10.0 * d / t);
    }
    
    return res;
}

// checkerbord
// https://www.shadertoy.com/view/XlcSz2
float checkersTextureGradBox( in vec2 p, in vec2 ddx, in vec2 ddy )
{
    // filter kernel
    vec2 w = max(abs(ddx), abs(ddy)) + 0.01;  
    // analytical integral (box filter)
    vec2 i = 2.0*(abs(fract((p-0.5*w)/2.0)-0.5)-abs(fract((p+0.5*w)/2.0)-0.5))/w;
    // xor pattern
    return 0.5 - 0.5*i.x*i.y;                  
}

vec3 sceneShade(vec2 mat, vec3 pos, vec3 ray, float depth, float maxD)
{
    vec3 col;
    vec3 sky = skyColor(ray);
    if (depth > maxD - 0.01) {
    	return sky;
    }
    float sha = shadow(pos, sunDir, 10.);
    vec3 norm = normal(pos);
    vec3 albedo = vec3(0.);
    if(mat.y == MAT_SPHERE1) {
        albedo = vec3(1., 0., 0.);
    } else if (mat.y == MAT_SPHERE2) {
        albedo = vec3(0., 1., 0.);
    } else if (mat.y == MAT_SPHERE3) {
        albedo = vec3(0., 0., 1.);
    } else if(mat.y == MAT_FLOOR) {
        vec2 ddx_uvw = dFdx( pos.xz ); 
        vec2 ddy_uvw = dFdy( pos.xz ); 
        float checker = checkersTextureGradBox( pos.xz, ddy_uvw, ddy_uvw );
        albedo = vec3(max(0.2, checker)) * vec3(.8,0.8,0.7) * 2.0;
    }
    
    float diffuse = clamp(dot(norm, sunDir), 0.0, 1.0) * sha * 2.0;
    col = albedo * (diffuse + 0.05);
    
    float fo = 1.0 - exp2(-0.0001 * depth * depth);
    vec3 fco = 0.65*vec3(0.4,0.65,1.0);
    col = mix( col, sky, fo );
    return col;
}

vec3 screenShade(vec2 mat, vec3 pos)
{
    vec3 ro = vec3(0., 0., 0.);
    vec3 ta = vec3(0., 0., 3.);
    vec3 fo = normalize(ta - ro);
    vec3 ri = normalize(cross(vec3(0.,1.,0.), fo));
    vec3 up = normalize(cross(fo,ri));
    mat3 cam = mat3(ri,up,fo);
    vec3 ray = cam * normalize(pos);
    float depth = 0.01;
    vec3 p = ro;
    sceneTrace(p, ray, mat, depth, 100.);
    vec3 col = vec3(0.);
    col = sceneShade(mat, p, ray, depth, 100.);
    
    float a1 = fract(rayDirectionAnim * 19.99999);
    float a2 = floor(rayDirectionAnim * 19.99999 + 1.0) / 20.;
    float a3 = floor(rayDirectionAnim * 19.99999) / 20.;
    
    float aspect = iResolution.y / iResolution.x;
    float halfAspect = aspect * 0.5;
    
    a1 = step(pos.x * halfAspect + 0.5, a1);
    a2 = step(-pos.y * 0.49 + 0.51+ 0.0, a2);
    a3 = step(-pos.y * 0.49 + 0.51 + 0.0, a3);
    //col *= min(screenAlpha + min(a2 * a1 + a3, 1.0), 1.0);
    
    // vec3 uvCoord = vec3(pow(clamp(mix(pos.xy*0.5+0.5, pos.xy, uvToP), 0.0, 1.0), vec2(2.2)), 0.0);
    vec3 uvCoord = vec3(1.);
    col = mix(col, uvCoord, colToUv - min(screenAlpha + min(a2 * a1 + a3, 1.0), 1.0));
    return col;
    //return vec3(pow(clamp(pos.xy, 0.0, 1.0), vec2(2.2)), 0.0);
}

vec4 gizmoShade(vec2 mat, vec3 p)
{
    vec4 col = vec4(0.);

    if(mat.y == MAT_CORNER) {
        col = vec4(1.,0.,0., cornersAlpha);
    } else if (mat.y == MAT_SCREENZ) {
    	col = vec4(0.05, 0.05, 1., screenZAlpha);
    } else if (mat.y == MAT_MARCHSPHERE) {
        float alpha = clamp(sphereAnim - sphereID, 0.0, 1.0);
        vec3 sc = mix(vec3(.0, .1, 3.), vec3(1.02, 1., 0.02), float(sphereID == 0. || sphereID == hitSphereID));
        // vec3 sc = mix(vec3(.0, .1, 3.), vec3(.02, 1., 1.02), float(sphereID == 0. || sphereID == hitSphereID));
        float inRange = 1.0; //sphereID <= hitSphereID ? 1.0 : 0.0;
    	col = vec4(sc, alpha * raySphereAlpha * inRange);
    } else if (mat.y == MAT_MARCHROUTE) {
        col = vec4(1., 0., 0., .9);
    } else if (mat.y == MAT_RAYDIRECTION) {
        col = vec4(1., 0., 0., .9);
    } else if (mat.y == MAT_HITPOINT) {
    	col = vec4(1.02, 1., .02, raySphereAlpha);
    }
    return col;
}

vec4 radiusShade(vec2 mat, vec3 p)
{
    vec4 col = vec4(0.);

    vec3 ray = vec3(0., 0., 1.);
    vec2 d = vec2(10000.);
    
    if(rayDirectionAnim2 > 0.0) {
        if(abs(rayRoute.y - 0.) <= 0.01) {
            ray = normalize(vec3(rayRoute.x, rayRoute.y, rayRoute.z));
            radiusAnim = 20.;
            routeAnim = 20.;
            sphereAnim = 20.;
        }
    }
    
    float t = 0.0;
    for(int i=0; i<20; i++) {
        vec3 pos = ray * t;
        vec2 dd = sceneSpheres(pos);
        d = vec2(sdSphere(p - pos, dd.x), MAT_MARCHSPHERE);

        vec2 sceneDist = sceneMap(pos);
        if (sceneDist.x < 0.001) {
            if(hitSphereID == 0.) {
              hitSphereID = float(i);
            }
            break;
        }

        float alpha2 = step(radiusAnim, float(i) + 2.);
        float alpha = clamp(radiusAnim - float(i) - 1., 0.0, 1.0);
        
        vec4 cirCol = vec4(.0, 0.05, 0.1, 0.9); //mix(vec4(.0, 0.05, 0.1, 0.9), vec4(0.2, 1., 1.4, .6) , alpha2);
        // fill
        col = mix(col, cirCol, smoothstep(0.01, 0., d.x) * cirCol.a * alpha);
        // strike
        col = mix(col, vec4(0., 0., 0., 1.), smoothstep(0.02, 0., abs(d.x) - 0.01) * alpha);
        t += dd.x;
    }
    col *= radiusAlpha;
    return col;
}

float luminance2(vec3 col)
{
    return dot(vec3(0.298912, 0.586611, 0.114478), col);
}

vec3 reinhard(vec3 col, float exposure, float white) {
    col *= exposure;
    white *= exposure;
    float lum = luminance2(col);
    return (col * (lum / (white * white) + 1.0) / (lum + 1.0));
}

vec3 render(vec2 p) {
    screenPos = vec3(p, 2.5);
    vec3 ro = cameraPos;
    vec3 ta = cameraTarget;
    vec3 fo = normalize(ta - ro);
    vec3 ri = normalize(cross(vec3(0.,1.,0.), fo));
    vec3 up = normalize(cross(fo,ri));
    mat3 cam = mat3(ri,up,fo);
    vec3 ray = cam * normalize(screenPos);
    float depth = 0.01;
    vec2 mat;
    vec3 col = vec3(0.);
    
    vec3 pos = ro;
    sceneTrace(pos, ray, mat, depth, 100.);
    
    col = sceneShade(mat, pos, ray, depth, 100.);
    
    float sceneDepth = depth;
    depth = 0.01;
    pos = ro;
    screenTrace(pos, ray, mat, depth, sceneDepth);
    if (depth < sceneDepth) {
        col = screenShade(mat, pos);
    }
    
    float sceneAndScreenDepth = depth; // sceneDepth
    depth = 0.01;
    pos = ro;
    radiusTrace(pos, ray, mat, depth, sceneAndScreenDepth);
    if (depth < sceneAndScreenDepth) {
        vec4 radius = radiusShade(mat, pos);
        col = mix(col, radius.rgb, radius.a);
    }
    
    float sceneAndScreenAndGizmoDepth = sceneAndScreenDepth;
    depth = 0.01;
    pos = ro;
    gizmoTrace(pos, ray, mat, depth, sceneAndScreenAndGizmoDepth);
    if (depth < sceneAndScreenAndGizmoDepth) {
        vec4 gizmo = gizmoShade(mat, pos);
        col = mix(col, gizmo.rgb, gizmo.a);
    }
    return col;
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 uv = (iMouse.xy - iResolution.xy)/min(iResolution.x, iResolution.y);
    // timeLine(mod(length(uv.x), 41.));
    timeLine(mod(iTime/4., 4.5));

    vec2 p = (fragCoord * 2.0 - iResolution.xy)/min(iResolution.x, iResolution.y);
    vec2 dp = 1. / iResolution.xy;

    vec3 col = vec3(0.);
    
    
    // AA
    // https://www.shadertoy.com/view/Msl3Rr
    for(int y = 0; y < 2; y++) {
        for(int x = 0; x < 2; x++) {
            vec2 off = vec2(float(x),float(y))/2.;
            vec2 xy = (-iResolution.xy+2.0*(fragCoord + off * dp / 4.0)) / iResolution.y;
        	col += render(xy)*0.25;
        }
    }
    
    col = reinhard(col, 1.0, 1000.0);
    col = pow(col, vec3(1.0/2.2));
    
    //col = mix(col, vec3(mix(uv, p, uvToP), 0.), colToUv);
    
    fragColor = vec4(col,1.0);
}
                                                                                                    adapted from kaneta's Raymarching in Raymarching
                                                                                                    add the algorithm and the shader for sphere tracing

                                                                                                    Theorem.  [gillespie2024ray, thm. 2] [ag-000U]

                                                                                                      Given a function \(F: \mathbb {R} \rightarrow \mathbb {R}\) with Lipschitz bound \(\lambda \geq \) Lip \(F\), and an initial point \(t_0\), sphere tracing converges linearly to the smallest root greater than \(t_0\).

                                                                                                      Remark.  [ag-000Q]

                                                                                                        Theorem [ag-000U] means that sphere tracing requires the function to be Lipschitz, i.e. has a known Lipschitz bound.

                                                                                                        Not all functions are Lipschitz, e.g. Harmonic Functions are not globally Lipschitz, many may have no local Lipschitz bound even within small neighborhoods, thus requires a new technique named Harnack tracing developed in [gillespie2024ray].