RaymarchedHyperboloid.glsl

// raymarched hyperboloid
//boilerplate
//#ifdef GL_ES
precision highp float;
//#endif

uniform float time;
uniform vec2 mouse;
uniform vec4 color;
uniform vec2 resolution;

float scale = 3.3333;// .5*exp(1./max(.00125,color.r));
vec3 iResolution = vec3(resolution.xy * scale * .5, 1.);
float iTime = time;
vec4 iDate = vec4(color.r, iTime*.125, iTime*.5, iTime);
vec4 iMouse=vec4(mouse*iResolution.xy,1.,1.);

float sphererad = 10.0;

//Blakes graphing calculator - test ver
float sdfOfBound(vec3, vec3, vec3);
float sgn(float);
float LineHeight(float, vec2, float);
vec2 intersectionOfTwoLines(vec2, float, vec2, float);

// So the idea of this shader is to raymarch some function given some information about it.
// The relevant information consists of 3 functions, along with some bounds.
// The three functions below define the function we are going to raymarch. df is the first
// derivative of f, and ddf is the second derivative of f. Along with that, some bounds need to
// be defined as well, but we'll get to that.


//////// HACKS!!!!!
float yscale=.7; // underestimate f(x) by this factor and scale point to counteract.
float ddffactor=40.; // overestimate ddf/dx by this factor (why?!)

float f(float x)
{
    // return the function when fed an x value
    // given the roots of the cubic, calculate the function of x
    return yscale*sqrt(1.+x*x);
}

float ff(float x)
{
    // return the function when fed an x value
    // given the roots of the cubic, calculate the function of x
    return sqrt(1.+x*x);
}
float df(float x)
{
    // return the derivative of the function when fed an x value
    // given the roots of the cubic, calculate the derivative of x
    return (2.)*x/ff(x);
}
float ddf(float x)
{
    return (ff(x)*2.-2.*x*(ddffactor)*x/ff(x))/(1.+x*x);
}
vec3 GetBoundary(int i, vec3 pt)
{

    float ex0 = -1000.0;
    float ex1 = -1000.;//-120000.;//-2.5/2./fscale;//3.1415926;
    float ex2 = 1000.;//2.5/2./fscale;//3.1415926;
    float ex3 = 1000.0;

    float i0 = -1000.0;
    float i1 = pt.z>0.?1./30.:-1./30.;
    float i2 = 1000.0;

    if (i == 0)
    {
        return vec3(ex0, i0, ex1);
    }
    if (i == 1)
    {
        return vec3(ex1, i1, ex2);
    }
    if (i == 2)
    {
        return vec3(ex2, i2, ex3);
    }
    return vec3(0.0,0.0,0.0);
}

float sdfFunction(vec3 Point, vec3 Center)
{
    // The number of extrema (not counting infinities
    const int EXTREMA_COUNT = 2;
        vec3 pnt = Point - Center;
    //pnt.xz = mod(pnt.xz + vec2(2.0), vec2(4.0))-vec2(2.0);
    float h = length(pnt.xz);
    h = h * sgn(pnt.x)*sgn(pnt.z);
        Point = vec3(0.,yscale*pnt.y, h) + Center;
vec3 invPoint = vec3(0.,-yscale*pnt.y, h) + Center;
    float d = 100000.0;
    for (int i = 0; i < EXTREMA_COUNT + 1; i++)
    {
        d = min(d, min(sdfOfBound(Point, Center, GetBoundary(i,Point)), sdfOfBound(invPoint, Center, GetBoundary(i,invPoint))));
    }
    return d;
}

///////////////////////

// End of user input //

///////////////////////

float sdfOfBound(vec3 Point, vec3 Center, vec3 Keys)
{
    vec3 pnt = Point - Center;
    //pnt.xz = mod(pnt.xz + vec2(1.0), vec2(2.0))-vec2(1.0);
    float h = pnt.z;//length(pnt.x);//pnt.x + pnt.z; //max(abs(pnt.x),abs(pnt.z));
    //h = h * sgn(pnt.x)*sgn(pnt.z);
    float k = pnt.y;

    float a = Keys.x;
    float b = Keys.z;
    float I = Keys.y;

    vec2 P_temp = vec2(h, k);
    h = clamp(h, a, b);
    vec2 P = vec2(h,k);

    bool tangent_method = sgn(f(h) - k) == sgn(ddf(h));

    float tan_at_I = df(I);
    float nor_at_I = -1.0 / tan_at_I;
    vec2 pnt_at_I = vec2(I, f(I));
    float line_height = LineHeight(nor_at_I, pnt_at_I, h);
    bool chord_method = (k*ddf(h) < line_height*ddf(h)) && !tangent_method;
    bool clamped_tangent = !chord_method && !tangent_method;

    float ct_slope = tan_at_I;
    vec2 ct_pnt = pnt_at_I;

    float t_slope = df(h);
    vec2 t_pnt = vec2(h, f(h));

    float ch_samp_d = abs(f(h) - k) * sgn(df(h)) * sgn(ddf(h));
    if (P_temp.x < a || P_temp.x > b)
    {
         ch_samp_d = abs(f(h) - k) * -sgn(P_temp.x-h);
    }
    float ch_raw_Q = h + ch_samp_d;
    float ch_clamped_Q = clamp(ch_raw_Q, min(I, ch_raw_Q), max(I, ch_raw_Q));
    vec2 ch_P = vec2(h, f(h));
    vec2 ch_Q = vec2(ch_clamped_Q, f(ch_clamped_Q));
    float ch_slope = (ch_P.y - ch_Q.y) / (ch_P.x - ch_Q.x);
    vec2 ch_pnt = ch_P;
    vec2 ct_closest = intersectionOfTwoLines(ct_pnt, ct_slope, P, -1.0 / ct_slope);
    vec2 t_closest = intersectionOfTwoLines(t_pnt, t_slope, P, -1.0 / t_slope);
    vec2 ch_closest = intersectionOfTwoLines(ch_pnt, ch_slope, P, -1.0 / ch_slope);

    float ct_dist = distance(ct_closest, P);
    float t_dist = distance(t_closest, P);
    float ch_dist = distance(ch_closest, P);

    float d = -1.0;
    if (clamped_tangent){ d = ct_dist;}
    if (tangent_method) { d = t_dist; }
    if (chord_method)   { d = ch_dist;}

   float dx = d;
    d = max(d , -(P_temp.x - a));
    d = max(d , (P_temp.x - b));
    //if (f(h) > k)         { d *= -1.0;  }
    return dx > 0. ? d: dx;//max(dx,-1111111.);
}


// Convert degrees to radians
float Radians(float deg)
{
    return deg / 360.0 * 2.0 * 3.14159;
}

// Write a float4 function for some of the HLSL Code conversion
vec4 float4(float x, float y, float z, float w)
{
    return vec4(x,y,z,w);
}

// Write a float3 function for the same purpose
vec3 float3(float x, float y, float z)
{
    return vec3(x,y,z);
}

// and a float2 function as well
vec2 float2(float x, float y)
{
    return vec2(x, y);
}

// A method for intersecting two lines in R-two.
vec2 intersectionOfTwoLines(vec2 A, float sa, vec2 B, float sb)
{
    float buffer = (sa*A.x - A.y - sb * B.x + B.y) / (sa - sb);
    vec2 intersection = float2(buffer, sa*(buffer - A.x) + A.y);
    return intersection;
}

// Exact SDF for a sphere
float dSphere(vec3 pos, vec3 center, float radius)
{
    // find the distance to the center
    vec3 v = pos - center;

    // return that, minus the radius
    return length(v) - radius;
}

// Exact intersection of a sphere. Resolves a quatratic equation. Returns the
// min distance, max distance, and discriminant to determine if the intersections
// actually exist.
vec3 intersections_of_sphere(vec3 pos_vector, vec3 dir_vector, float sphere_radius)
{
    // Derivation for formula:
    //      Let the ray be represented as a point P plus a scalar multiple t of the direction vector v,
    //      The ray can then be expressed as P + vt
    //
    //      The point of intersection I = (x, y, z) must be expressed as this, but must also be some distance r
    //      from the center of the sphere, thus x*x + y*y + z*z = r*r, or in vector notation, I*I = r*r
    //
    //      It therefore follows that (P + vt)*(P + vt) = r*r, or when expanded and rearranged,
    //      (v*v)t^2 + (2P*v)t + (P*P - r*r) = 0. For this we will use the quadratic equation for the points of
    //      intersection

    // a, b, and c correspond to the second, first, and zeroth order terms of t, the parameter we are trying to solve for.
    float a = dot(dir_vector, dir_vector);
    float b = 2.0 * dot(pos_vector, dir_vector);
    float c = dot(pos_vector, pos_vector) - sphere_radius * sphere_radius;

    // to avoid imaginary number, we will find the absolute value of the discriminant.
    float discriminant = b * b - 4.0 * a*c;
    float abs_discriminant = abs(discriminant);
    float min_dist = (-b - sqrt(abs_discriminant)) / (2.0 * a);
    float max_dist = (-b + sqrt(abs_discriminant)) / (2.0 * a);

    // return the two intersections, along with the discriminant to determine if
    // the intersections actually exist.
    return float3(min_dist, max_dist, discriminant);

}

float sgn(float x)
{
    if (x < 0.0)
    {
        return -1.0;
    }
    return 1.0;
}

float LineHeight(float slope, vec2 point, float x)
{
    return slope*(x - point.x) + point.y;
}

// Distance estimation for the scene
float DE(vec3 p, vec3 c)
{
    // Raymarch a bounding sphere
    float sphere = dSphere(p, c, sphererad);

    // Raymarch generic sdf for a function
    float function = sdfFunction(p,c);

    // cut out the bounding sphere
    float d = max(sphere, function);

    // return the distance
    return d;
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    // Define the iterations for the marcher.
    const int ITERATIONS = 60;

    // Define the roation speed. Set to 0 to disable
    const float ROTATION_SPEED = 0.2;

    // Define the start angle for the rotation (in degrees)
    const float START_ANGLE = 0.0;

    // Define the orbit radius
    const float ORBIT_RADIUS = 16.0;

    // Define the epsilon value for closeness to be considered a hit
    const float EPSILON = 0.001;

    const bool HIDE_BACKGROUND = false;



    // Define the center of the julia set
    vec3 sdf_center = vec3(0.0, 0.0, 0.0);

    // Calculate the starting angles for the orbit
    float theta = iTime * ROTATION_SPEED;
    float phi = Radians(START_ANGLE);

    // Take some mouse input
    vec4 mouse = iMouse / iResolution.xyxx;

    // If the mouse is being held down
    //if (mouse.z > 0.0)
    {
        // convert the mouse input to angles
        theta = mouse.x * 2.0 * 3.14159;
        phi = (mouse.y - 0.5) * 1.0 * 3.14159;
    }

    // Define an orbital path based on time
    vec3 orbit = vec3(cos(theta)*cos(phi), sin(phi), sin(theta)*cos(phi));

    // Cacluate the normal of the path. Since its a circle, it will just
    // be back down into the center
    vec3 normal = -normalize(orbit);

    // Calculate the tangent of the path
    // A circle consists of <cost, sint>, which when differentiated yields
    // <-sint, cost>. since z is already sint, and x is already cost, the equation
    // is as follows.
    vec3 tangent = normalize(vec3(-normal.z, 0.0, normal.x));

    // Calculate the UV coordinates
    vec2 uv = fragCoord/iResolution.xy;

    // Convert the UV coordinates to a range between -1 and 1
    vec2 range = uv*2.0 - vec2(1.0,1.0);

    //// Define the Camera position
    //vec3 cam_pos = vec3(0,0,-2);

    //// Define the forward, up, and right vectors (needs rework)
    //vec3 forward = normalize(vec3(0,0,1));
    //vec3 up = normalize(vec3(0,1,0));
    //vec3 right = normalize(vec3(1,0,0));

    // Define the Camera position
    vec3 cam_pos = orbit*ORBIT_RADIUS;

    // Define the forward, up, and right vectors (needs rework)
    vec3 forward = normal;
    vec3 up = normalize(cross(normal, tangent));
    vec3 right = tangent;

    // Calculate the aspect ratio of the screen
    float aspect = float(iResolution.y) / float(iResolution.x);

    // Calculate the ray as a normalized combination of the forward, right, and up vectors.
    // Note that the purely forward + horizonal combination yield vectors 45 degrees outward
    // for a 90 degree field of view. This may be updated with a fov option
    vec3 ray = normalize(forward + range.x*right + range.y*up*aspect);

    // Initialize the ray marched point p
    vec3 p = cam_pos;


    // Initialize the distance
    float dist = 1.0;

    // Calculate the exact distance from a sphere of radius 2 using a raytracing function
    vec3 init_distance = intersections_of_sphere(p - sdf_center, ray, sphererad);

    // If we are outside a bubble around the raymarched fractal
    if (init_distance.z > 0.0)
    {
        // Step onto the sphere so we start off a bit closer.
        p += ray * clamp(init_distance.x, 0.0, init_distance.x);
    }

    // declare a dummy variable to store the number of iterations into.
    // I'm doing it this way because on my phone it didnt let me use an
    // already declared variable as the loop iterator.
    int j;

    float min_dist = 1.0;
    // Begin the raymarch
    for (int i = 0; i < ITERATIONS; i++)
    {
        // Estimate the distance to the julia set
        dist = DE(p, sdf_center);

        min_dist = min(dist, min_dist);

        // Move forward that distance
        p += ray*dist;

        // Record the number of iterations we are on
        j = i;

        // If we hit the julia set, or get too far away form it
        if (dist < EPSILON || dot(p - sdf_center, p-sdf_center) > sphererad * sphererad + 0.1)
        {
            // Break the loop.
            break;
        }

    }

        // determine if we hit the fractal or not
    float hit = step(dist, EPSILON);


    // calculate the brightness based on iterations used
    float di = (1.0 - (float(j) + (dist / EPSILON)*hit) / float(ITERATIONS));
    float glow = 1.0 - (min_dist / 1.0);

    //di = (di * hit) + (glow*(1.0-hit));
    if (HIDE_BACKGROUND)
    {
        di *= hit;
    }




    // define some phase angle
    float psi = Radians(70.0);

    // Time varying pixel color (included in default shadertoy project)
    //vec3 col = 0.8 + 0.2*cos(iTime*0.5+uv.xyx+vec3(0,2,4) + psi*hit);

    // Boring old white instead of the above commented code. Will tweak rendering later
    vec3 col = vec3(1.0,1.0,1.0);


    // Output to screen. Modifiy the color with the brightness calculated as di.
    fragColor = vec4(col*di,1.0);
}

void main( void ) {
    gl_FragColor = vec4(0.,10.,0.,0.);
    vec2 zoomedCoord = (color.gb - .5) * .125 * iResolution.y - (mouse.xy-vec2(.5,.5)) * iResolution.xy + gl_FragCoord.xy - resolution*.5;
zoomedCoord = gl_FragCoord.xy + resolution*.33;
    mainImage( gl_FragColor, zoomedCoord.xy );
}
// P.S. vrchat is blessed. ~Ly