tinkering glsl

tinkering.glsl

// Licence CC0: Tunneling through appollian fractal

// -----------------------------------------------------------------------------
// COMMON
// -----------------------------------------------------------------------------

#define PI              3.141592654
#define TAU             (2.0*PI)
#define TIME            time
#define RESOLUTION      resolution
#define ROT(a)          mat2(cos(a), sin(a), -sin(a), cos(a))
#define PSIN(x)         (0.5+0.5*sin(x))
#define LESS(a,b,c)     mix(a,b,step(0.,c))
#define SABS(x,k)       LESS((.5/(k))*(x)*(x)+(k)*.5,abs(x),abs(x)-(k))
#define L2(x)           dot(x, x)

const vec3 std_gamma        = vec3(2.2, 2.2, 2.2);

float hash(float co) {
  co += 100.0;
  return fract(sin(co*12.9898) * 13758.5453);
}


float hash(vec2 co) {
  co += 100.0;
  return fract(sin(dot(co, vec2(12.9898,58.233))) * 13758.5453);
}



vec2 toPolar(vec2 p) {
  return vec2(length(p), atan(p.y, p.x));
}

vec2 toRect(vec2 p) {
  return vec2(p.x*cos(p.y), p.x*sin(p.y));
}

float modMirror1(inout float p, float size) {
  float halfsize = size*0.5;
  float c = floor((p + halfsize)/size);
  p = mod(p + halfsize,size) - halfsize;
  p *= mod(c, 2.0)*2.0 - 1.0;
  return c;
}


vec4 modMirror4(inout vec4 p, vec4 size) {
  vec4 halfsize = size*0.5;
  vec4 c = floor((p + halfsize)/size);
  p = mod(p + halfsize,size) - halfsize;
  p *= mod(c, 2.0)*2.0 - 1.0;
  return c;
}


float smoothKaleidoscope(inout vec2 p, float sm, float rep) {
  vec2 hp = p;

  vec2 hpp = toPolar(hp);
  float rn = modMirror1(hpp.y, TAU/rep);

  float sa = PI/rep - SABS(PI/rep - abs(hpp.y), sm);
  hpp.y = sign(hpp.y)*(sa);

  hp = toRect(hpp);

  p = hp;

  return rn;
}

vec4 alphaBlend(vec4 back, vec4 front) {
  float w = front.w + back.w*(1.0-front.w);
  vec3 xyz = (front.xyz*front.w + back.xyz*back.w*(1.0-front.w))/w;
  return w > 0.0 ? vec4(xyz, w) : vec4(0.0);
}

vec3 alphaBlend(vec3 back, vec4 front) {
  return mix(back, front.xyz, front.w);
}

float tanh_approx(float x) {
//  return tanh(x);
  float x2 = x*x;
  return clamp(x*(27.0 + x2)/(27.0+9.0*x2), -1.0, 1.0);
}

float emin(float a, float b, float k) {
  float res = exp2(-k*a) + exp2(-k*b);
  return -log2( res )/k;
}

float pmin(float a, float b, float k) {
  float h = clamp(0.5+0.5*(b-a)/k, 0.0, 1.0);
  
  return mix(b, a, h) - k*h*(1.0-h);
}

float pmax(float a, float b, float k) {
  return -pmin(-a, -b, k); 
}

float pwmin(float a, float b, float k) {
  a = pow(a, k); 
  b = pow(b, k);
  return pow((a*b)/(a+b), 1.0/k);
}

vec3 hsv2rgb(vec3 c) {
  const vec4 K = vec4(1.0, 2.0 / 3.0, 1.0 / 3.0, 3.0);
  vec3 p = abs(fract(c.xxx + K.xyz) * 6.0 - K.www);
  return c.z * mix(K.xxx, clamp(p - K.xxx, 0.0, 1.0), c.y);
}

float circle(vec2 p, float r) {
  return length(p) - r;
}

float sphere(vec4 p, float r) {
  return length(p) - r;
}


float hex(vec2 p, float r) {
  const vec3 k = vec3(-sqrt(3.0)/2.0,1.0/2.0,sqrt(3.0)/3.0);
  p = p.yx;
  p = abs(p);
  p -= 2.0*min(dot(k.xy,p),0.0)*k.xy;
  p -= vec2(clamp(p.x, -k.z*r, k.z*r), r);
  return length(p)*sign(p.y);
}

float apollian(vec4 p, float s, float koff) {
  float scale = 1.0;

  for(int i=0; i<7; ++i) {
    p = -1.0 + 2.0*fract(0.5*p+0.5);

//    float r2 = pow(dot(p,p), 0.125*i+1);
    float r2 = dot(p,p);
    
    float k  = s/r2;
// Weird stuff
//    const float ss = 50.0;
//    p       *= ss*tanh(k/ss)-0.225;
//    p       *= ss*tanh(k/ss)-0.225;
//    p       *= k+0.225*tanh_approx(k);
//    p       *= k*tanh(2.1/r2);
    p       *= k-koff;
    scale   *= k;
  }

  const float lw = 0.00125*.25;
  
  float d0 = abs(p.y)-lw*scale;
  float d1 = abs(circle(p.xz, 0.005*scale))-lw*scale;
//  float d0 = abs(p.x);
  float d = d0;
  d = min(d, d1);
  return (d/scale);
}

vec2 mod2_1(inout vec2 p) {
  vec2 c = floor(p + 0.5);
  p = fract(p + 0.5) - 0.5;
  return c;
}

// -----------------------------------------------------------------------------
// PATH
// -----------------------------------------------------------------------------

// The path function
vec3 offset(float z) {
  float a = z;
  vec2 p = -0.075*(vec2(cos(a), sin(a*sqrt(2.0))) + vec2(cos(a*sqrt(0.75)), sin(a*sqrt(0.5))));
  return vec3(p, z);
}

// The derivate of the path function
//  Used to generate where we are looking
vec3 doffset(float z) {
  float eps = 0.1;
  return 0.5*(offset(z + eps) - offset(z - eps))/eps;
}

// The second derivate of the path function
//  Used to generate tilt
vec3 ddoffset(float z) {
  float eps = 0.1;
  return 0.125*(doffset(z + eps) - doffset(z - eps))/eps;
}

// -----------------------------------------------------------------------------
// PLANE MARCHER
// -----------------------------------------------------------------------------

const float fixed_radius2 = 1.9;
const float min_radius2   = 0.5;
const float folding_limit = 1.0;
const float scale         = -2.8;

void sphere_fold(inout vec3 z, inout float dz) {
    float r2 = dot(z, z);
    if(r2 < min_radius2) {
        float temp = (fixed_radius2 / min_radius2);
        z *= temp;
        dz *= temp;
    } else if(r2 < fixed_radius2) {
        float temp = (fixed_radius2 / r2);
        z *= temp;
        dz *= temp;
    }
}

vec3 pmin(vec3 a, vec3 b, vec3 k) {
  vec3 h = clamp( 0.5+0.5*(b-a)/k, 0.0, 1.0);
  
  return mix(b, a, h) - k*h*(1.0-h);
}

void box_fold(inout vec3 z, inout float dz) {
  const float k = 0.05;
  // Soft clamp after suggestion from ollij
  vec3 zz = sign(z)*pmin(abs(z), vec3(folding_limit), vec3(k));
  // Hard clamp
  // z = clamp(z, -folding_limit, folding_limit);
  z = zz * 2.0 - z;
}

float sphere(vec3 p, float t) {
  return length(p)-t;
}

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

float mb(vec3 z) {
    vec3 offset = z;
    float dr = 1.0;
    float fd = 0.0;
    for(int n = 0; n < 5; ++n) {
        box_fold(z, dr);
        sphere_fold(z, dr);
        z = scale * z + offset;
        dr = dr * abs(scale) + 1.0;        
        float r1 = sphere(z, 5.0);
        float r2 = torus(z, vec2(8.0, 1));
        r2 = abs(r2) - 0.25;
        float r = n < 4 ? r2 : r1;        
        float dd = r / abs(dr);
        if (n < 3 || dd < fd) {
          fd = dd;
        }
    }
    return fd;
}

vec2 mod2(inout vec2 p, vec2 size) {
  vec2 c = floor((p + size*0.5)/size);
  p = mod(p + size*0.5,size) - size*0.5;
  return c;
}

float weird2(vec2 p, float h) {
  float tm = 0.25*TIME+100*h;
//  tm *= 0.0;
  const float s = 1.0-0.1;
  p /= s;
  float rep = 20.0;
  float ss = 0.05*6.0/rep;
  vec3 p3 = vec3(p.x, p.y, PSIN(tm*sqrt(0.5)));
  p3.yz *= ROT(tm);
  float n = smoothKaleidoscope(p3.xy, ss, rep);
  return mb(p3)*s;
}


float weird(vec2 p, float h) {
  float z = 4.0;
  float tm = 0.1*TIME+h*10;
  p *= ROT(tm*0.5);
  float r = 0.5;
  vec4 off = vec4(r*PSIN(tm*sqrt(3.0)), r*PSIN(tm*sqrt(1.5)), r*PSIN(tm*sqrt(2.0)), 0.0);
  vec4 pp = vec4(p.x, p.y, 0.0, 0.0)+off;
  pp.w = 0.125*(1.0-tanh_approx(length(pp.xyz)));
  pp.yz *= ROT(tm);
  pp.xz *= ROT(tm*sqrt(0.5));
  pp /= z;
  float d = apollian(pp, 0.8+1*h, 0.0);
  return d*z;
}

float circles(vec2 p) {
  vec2 pp = toPolar(p);
  const float ss = 2.0;
  pp.x = fract(pp.x/ss)*ss;
  p = toRect(pp);
  float d = circle(p, 1.0);
  return d;
}

vec2 df(vec2 p, float h) {
  vec2 wp = p;
  float rep = 2.0*int(mix(5.0, 10.0, h));
  rep = 10.0;
  float ss = 0.05*6.0/rep;
  float n = smoothKaleidoscope(wp, ss, rep);
  
  float d0 = weird(wp, h);
  float d1 = hex(p, 0.25)-0.1;
  float d2 = circles(p);
  const float lw = 0.0125;
  d2 = abs(d2)-lw;
  float d = d0;
  d = pmin(d, d2, 0.1);
  d  = pmin(d, abs(d1)-lw, 0.1);
  return vec2(d, d1+lw);
}

vec2 df(vec3 p, vec3 off, float s, mat2 rot, float h) {
  vec2 p2 = p.xy;
  p2 -= off.xy;
  p2 *= rot;
  return df(p2/s, h)*s;
}

// Plane generating function returns rgba
//  pp is point on plane
//  off is path at plane z
//  aa is estimated pixel size at the plane
//  n is plane number
vec4 plane1(vec3 ro, vec3 rd, vec3 pp, float pd, vec3 off, float aa, float n) {
  float h = hash(n);
  float s = 0.25*mix(0.5, 0.25, h);

  const vec3 nor  = vec3(0.0, 0.0, 1.0);
  const vec3 loff = vec3(0.25*0.5, 0.125*0.5, -0.125);
  vec3 lp1  = ro + loff;
  vec3 lp2  = ro + loff*vec3(-1.0, 1.0, 1.0);
  vec3 ld1  = normalize(pp - lp1);
  vec3 ld2  = normalize(pp - lp2);
  float lpw1= 0.2/L2(pp - lp1);
  float lpw2= 0.2/L2(pp - lp2);
  vec3 ref  = reflect(rd, nor);
  float ref1= pow(max(dot(nor, ld1), 0.0), 20.0);
  float ref2= pow(max(dot(nor, ld2), 0.0), 20.0);
  vec3  col1= vec3(0.75, 0.5, 1.0);
  vec3  col2= vec3(1.0, 0.5, 0.75);
  
  vec3 hn;
  mat2 rot = ROT(TAU*h);
//  p *= rot;
  vec2 d2 = df(pp, off, s, rot, h);   

  float ha = smoothstep(-aa, aa, d2.y);
  float d = d2.x;
  vec4 col = vec4(0.0);

  vec2 p = (pp-off*vec3(1.0, 1.0, 0.0)).xy;
  float l   = length(10*p);
  float pdf = 1.0/(1*(1+2.0*pd));
  float hue = fract(0.75*l-0.1*TIME)+0.3+0.15;
  float sat = 0.75*tanh_approx(2.0*l)*pdf;
  float vue = sqrt(pdf);
  vec3 hsv  = vec3(hue, sat, vue);
  vec3 bcol = hsv2rgb(hsv);
  col.xyz   = mix(col.xyz, bcol, smoothstep(-aa, aa, -d));  
  float glow = (exp(-(10.0+100.0*tanh_approx(l))*10*max(d, 0.0)*pdf));
  col.xyz   += 0.5*sqrt(bcol.zxy)*glow;
  col.w     = ha*mix(0.75, 1.0, ha*glow);
  col.xyz   += 0.125*col.w*(col1*ref1+col2*ref2);

  return col;
}

vec4 plane2(vec3 ro, vec3 rd, vec3 pp, float pd, vec3 off, float aa, float n) {
  float h = hash(n);
  float s = 0.25*mix(0.5, 0.25, h);
  const float lw = 0.0235;
  const float lh = 1.25;

  const vec3 nor  = vec3(0.0, 0.0, 1.0);
  const vec3 loff = 2*vec3(0.25*0.5, 0.125*0.5, -0.125);
  vec3 lp1  = ro + loff;
  vec3 lp2  = ro + loff*vec3(-2.0, 1.0, 1.0);

  vec2 p = pp.xy-off.xy;

  mat2 rot = ROT(TAU*h);
  const bool scale = true;

  vec2 d2 = df(pp, off, s, rot, h);

  vec3 ld1 = normalize(lp1 - pp);
  vec3 ld2 = normalize(lp2 - pp);

  const float boff = 0.0125*0.5;
  float dbt = boff/rd.z;
  
  vec3 bpp = ro + (pd + dbt)*rd;
  vec2 bp = bpp.xy - off.xy;

  vec3 srd1 = normalize(lp1-bpp);
  vec3 srd2 = normalize(lp2-bpp);
  float bl21= pow(L2(lp1-bpp), 1);
  float bl22= L2(lp2-bpp);

  float st1 = -boff/srd1.z;
  float st2 = -boff/srd2.z;

  vec3 spp1 = bpp + st1*srd1;
  vec3 spp2 = bpp + st2*srd2;
  
  vec2 sd1 = df(spp1, off, s, rot, h);
  vec2 sd2 = df(spp2, off, s, rot, h);

  vec3 col  = vec3(0.0);
  const float ss = 200.0;

  const vec3 lpw = vec3(0.125);
  col       += lpw*(1.0-exp(-ss*(max((sd1.x), 0.0))))/bl21;
  col       += lpw.zxy*0.5*(1.0-exp(-ss*(max((sd2.x), 0.0))))/bl22;
  col       = mix(col, vec3(1.0), smoothstep(-aa, aa, -(d2.x)));  

  float ha = smoothstep(-aa, aa, d2.y);

  return vec4(col, mix(0.0, 1.0, ha));
}

vec3 skyColor(vec3 ro, vec3 rd) {
  float ld = max(dot(rd, vec3(0.0, 0.0, 1.0)), 0.0);
  return vec3(1.25, 1.0, 1.1)*tanh_approx(3.0*pow(ld, 100.0));
}

vec3 color(vec3 ww, vec3 uu, vec3 vv, vec3 ro, vec2 p) {
  float lp = length(p);
  vec2 np = p + 1.0/RESOLUTION.xy;
  float rdd = (2.0+0.5*tanh_approx(lp));  // Playing around with rdd can give interesting distortions
  vec3 rd = normalize(p.x*uu + p.y*vv + rdd*ww);
  vec3 nrd = normalize(np.x*uu + np.y*vv + rdd*ww);

  const float planeDist = 1.0-0.75;
  const int furthest = 9;
  const int fadeFrom = max(furthest-4, 0);
  const float fadeDist = planeDist*(furthest - fadeFrom);
  float nz = floor(ro.z / planeDist);

  vec3 skyCol = skyColor(ro, rd);

  // Steps from nearest to furthest plane and accumulates the color

  vec4 acol = vec4(0.0);
  const float cutOff = 0.95;
  bool cutOut = false;
  
  for (int i = 1; i <= furthest; ++i) {
    float pz = planeDist*nz + planeDist*float(i);

    float pd = (pz - ro.z)/rd.z;

    if (pd > 0.0 && acol.w < cutOff) {
      vec3 pp = ro + rd*pd;
      vec3 npp = ro + nrd*pd;

      float aa = 3.0*length(pp - npp);

      vec3 off = offset(pp.z);

      vec4 pcol = plane1(ro, rd, pp, pd, off, aa, nz+float(i));

      float nz = pp.z-ro.z;
      float fadeIn = exp(-2.5*max((nz - planeDist*float(fadeFrom))/fadeDist, 0.0));
      float fadeOut = smoothstep(0.0, planeDist*0.1, nz);
      pcol.xyz = mix(skyCol, pcol.xyz, (fadeIn));
      pcol.w *= fadeOut;

      pcol = clamp(pcol, 0.0, 1.0);

      acol = alphaBlend(pcol, acol);
    } else {
      cutOut = true;
      break;
    }

  }

  vec3 col = alphaBlend(skyCol, acol);
// To debug cutouts due to transparency  
//  col += cutOut ? vec3(1.0, -1.0, 0.0) : vec3(0.0);
  return col;
}

// Classic post processing
vec3 postProcess(vec3 col, vec2 q) {
  col = clamp(col, 0.0, 1.0);
  col = pow(col, 1.0/std_gamma);
  col = col*0.6+0.4*col*col*(3.0-2.0*col);
  col = mix(col, vec3(dot(col, vec3(0.33))), -0.4);
  col *=0.5+0.5*pow(19.0*q.x*q.y*(1.0-q.x)*(1.0-q.y),0.7);
  return col;
}

vec3 effect(vec2 p, vec2 q) {
  float tm  = TIME*0.2+0.3;
  vec3 ro   = offset(tm);
  vec3 dro  = doffset(tm);
  vec3 ddro = ddoffset(tm);

  vec3 ww = normalize(dro);
  vec3 uu = normalize(cross(normalize(vec3(0.0,1.0,0.0)+ddro), ww));
  vec3 vv = normalize(cross(ww, uu));

  vec3 col = color(ww, uu, vv, ro, p);
  /*
  col = 1.0-col;
  col *= col;
  col *= col;
  */
//  col *= col;
//  col = 1.0-tanh(2.5*col);
//  col = 1.0-tanh(3.5*col-vec3(0.5, 0.25, .75)*(1.0-length(p)));
  col = postProcess(col, q);
  return col;
}

void main(void) {
  vec2 q = inData.v_texcoord;
  vec2 p = -1. + 2. * q;
  p.x *= RESOLUTION.x/RESOLUTION.y;

  vec3 col = effect(p, q);

  fragColor = vec4(col, 1.0);
}