From ⧉:

float SDSphere(in vec3 worldCoordinate, in vec3 center, in float radius) {

    return distance(worldCoordinate, center) - radius;

}

float SDPlane(in vec3 worldCoordinate, in float planeHeight, in vec3 rayDirection) {

    return (worldCoordinate.y - planeHeight)/abs(rayDirection.y);

}

float SDOctahedron( vec3 p, float s)
{
  p = abs(p);
  return (p.x+p.y+p.z-s)*0.57735027;
}


float SDScene(in vec3 worldCoordinate, in vec3 rayDirection) {

    float sphereOne = SDSphere(worldCoordinate, vec3(0.0, 0.0, 0.0), 0.5);
    //float sphereTwo = SDSphere(worldCoordinate, vec3(cos(iTime/3.0), -0.3, sin(iTime/5.0)), 0.2);
    //float sphereThree = SDSphere(worldCoordinate, vec3(sin(iTime/5.0), 0.3, cos(iTime/5.0)), 0.2);
    float orbital = SDOctahedron(worldCoordinate - vec3(cos(iTime/2.0), -0.15, sin(iTime/2.0)), 0.2);
    
    float groundPlane = SDPlane(worldCoordinate, -0.5, rayDirection);

    return min(min(sphereOne, orbital), groundPlane);

}


# define EPSILON 0.001

vec3 estimateNormal(in vec3 worldCoordinate) {

    if (worldCoordinate.y < -0.49) return vec3(0.0, 1.0, 0.0);

    float partialX = SDScene(worldCoordinate + vec3(EPSILON, 0.0, 0.0), vec3(0.0, 0.0, 1.0)) - SDScene(worldCoordinate - vec3(EPSILON, 0.0, 0.0), vec3(0.0, 0.0, 1.0));
    float partialY = SDScene(worldCoordinate + vec3(0.0, EPSILON, 0.0), vec3(0.0, 0.0, 1.0)) - SDScene(worldCoordinate - vec3(0.0, EPSILON, 0.0), vec3(0.0, 0.0, 1.0));
    float partialZ = SDScene(worldCoordinate + vec3(0.0, 0.0, EPSILON), vec3(0.0, 0.0, 1.0)) - SDScene(worldCoordinate - vec3(0.0, 0.0, EPSILON), vec3(0.0, 0.0, 1.0));
    
    return normalize(vec3(partialX, partialY, partialZ));

}


# define MAX_STEPS 105

vec3 rayMarch(in vec3 rayOrigin, in vec3 rayDirection) {

    vec3 ray = rayOrigin;

    for (int index = 0; index < MAX_STEPS; index++) {
   
        float distanceToScene = SDScene(ray, rayDirection);
        
        if (distanceToScene < EPSILON) return ray;
        
        ray += distanceToScene * rayDirection;
    
    }
    
    return vec3(0.0, 0.0, 0.0);

}


void mainImage(out vec4 pixelColor, in vec2 pixelCoordinate) {
    
    float aspect = iResolution.x/iResolution.y;
    
    vec2 uv = pixelCoordinate/iResolution.xy;
    
    uv -= vec2(0.5, 0.5);
    uv *= 2.0 * vec2(aspect, 1.0);
    
    vec3 observerPosition = vec3(0.0, 0.0, -5.0);
    
    vec3 cameraBox = observerPosition + vec3(uv, 5.0);
    
    vec3 rayDirection = normalize(cameraBox - observerPosition);
    vec3 worldCoordinate = rayMarch(observerPosition, rayDirection);
    
    float renderable = length(worldCoordinate) > 0.0 ? 1.0 : 0.0;
    
    vec3 surfaceNormal = estimateNormal(worldCoordinate);
    float fresnel = 0.9 + dot(surfaceNormal, rayDirection);
    
    vec3 lightSource = vec3(7.0 * cos(iTime), 5.0, 7.0 * sin(iTime));
    
    vec3 photonDirection = normalize(worldCoordinate - lightSource);
    vec3 photonPosition = rayMarch(lightSource, photonDirection);
    
    float directLight = 1.2 - step(0.01, distance(photonPosition, worldCoordinate));
    
    float diffuseLight = 0.25 * smoothstep(0.0, 1.0, directLight * dot(surfaceNormal, -photonDirection));
    
    float ambientLight = 0.7;

    float specularLight = 0.1 * smoothstep(0.0, 1.0, dot(reflect(photonDirection, surfaceNormal), normalize(observerPosition - worldCoordinate)));

    float scatterFactor = 1.0 / sqrt(distance(worldCoordinate, lightSource)/5.0);

    pixelColor = vec4(renderable * (diffuseLight + ambientLight + specularLight) * scatterFactor * vec3(1.0, 1.0, 1.0), 1.0);
    
}

From ⧉:

#define AA 1

//------------------------------------------------------------------
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.75 );
    vec3 ro = ta + vec3( 4.5*cos(0.1*time + 7.0*mo.x), 2.2, 4.5*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 );
}

The following renders texts and a ray-marching fragment shader on a plane which always faces the camera, the cubes are outside the plane.

// Adapted from https://lygia.xyz/lighting/raymarch

#define LIGHT_COLOR vec3(1.000,0.794,0.066)
#define RESOLUTION          iResolution.xy
#define RAYMARCH_MULTISAMPLE 4
// #define RAYMARCH_BACKGROUND (RAYMARCH_AMBIENT + rayDirection.y * 0.8)
#define RAYMARCH_BACKGROUND     vec3(1.0, 1.0, 1.0)
#define RAYMARCH_AMBIENT    vec3(1.0, 1.0, 1.0)
// #define RAYMARCH_AMBIENT    vec3(0.7, 0.9, 1.0)
// #define RAYMARCH_AMBIENT    vec3(0.374,0.481,0.535)

#define EPSILON 0.0001

#include "lygia/space/ratio.glsl"
#include "lygia/sdf.glsl"
#include "lygia/lighting/raymarch.glsl"
#include "lygia/color/space/linear2gamma.glsl"

float checkBoard(vec2 uv, vec2 _scale) {
    uv = floor(fract(uv * _scale) * 2.0);
    return min(1.0, uv.x + uv.y) - (uv.x * uv.y);
}

float myCubeSDF(vec3 p) {
    float x = p.x;
    float y = p.y;
    float z = p.z;
    float x6 = pow(x, 6.);
    float y6 = pow(y, 6.);
    float z6 = pow(z, 6.);

    return x6+y6+z6-1.;
}

float mySDF2(vec3 p) {
    float x = p.x;
    float y = p.y;
    float z = p.z;

    float x2 = pow(x, 2.);
    float y2 = pow(y, 2.);
    float z2 = pow(z, 2.);

    float x3 = pow(x, 3.);
    float y3 = pow(y, 3.);
    float z3 = pow(z, 3.);

    float x4 = pow(x, 4.);
    float y4 = pow(y, 4.);
    float z4 = pow(z, 4.);

    // return (y-x2)*(z-x3);
    
    // return 5.*(z2+y3-y4-x2*y2);

    // return (y2-x2-z2);

    return (x2-z2*y2+y3);
}

float tapping_abs(float x) {
    return max(-0.01 * x, x);
}

Material raymarchMap( in vec3 pos ) {
    // float check = 0.5 + checkBoard(pos.xz, vec2(1.0, 1.0)) * 0.5;
    // Material res = materialNew(vec3(check), 0.0, 0.5, planeSDF(pos - vec3(0.0, -2.0, 0.0)));
    
    // res = opUnion( res, materialNew( vec3(1.0, 2.0, 1.0), myCubeSDF(pos - vec3(0.0, -1.0, 0.0))*0.004 ) );

    // res = opUnion( res, materialNew( vec3(1.0, 2.0, 1.0),opIntersection(0.004 * myCubeSDF(pos), cubeSDF(pos, 1.0))) );

    // res = opUnion(res, materialNew( vec3(1.0, 2.0, 1.0), 0.004 * abs(mySDF2(pos))));

    // res = opUnion(res, materialNew( vec3(1.0, 2.0, 1.0), opIntersection(0.05 * abs(mySDF2(pos)), cubeSDF(pos, 1.0))) );
    
    Material res = materialNew( vec3(1.0, 2.0, 1.0), opIntersection(0.05 * abs(mySDF2(pos)), sphereSDF(pos, 1.8)));

    // res = opUnion( res, materialNew( vec3(1.0, 1.0, 1.0), 1.0, 0.0, sphereSDF(pos-vec3( 0.0, 2.0, 0.0), 0.5 ) ) );
    // res = opUnion( res, materialNew( vec3(0.0, 1.0, 1.0), boxSDF(      pos-vec3( 2.0, 0.5, 0.0), vec3(0.4, 0.4, 0.4) ) ) );
    // res = opUnion( res, materialNew( vec3(0.3, 0.3, 1.0), torusSDF(    pos-vec3( 0.0, 0.5, 2.0), vec2(0.4,0.1) ) ) );
    // res = opUnion( res, materialNew( vec3(0.3, 0.1, 0.3), capsuleSDF(  pos,vec3(-2.3, 0.4,-0.2), vec3(-1.6,0.75,0.2), 0.2 ) ) );
    // res = opUnion( res, materialNew( vec3(0.5, 0.3, 0.4), triPrismSDF( pos-vec3(-2.0, 0.50,-2.0), vec2(0.5,0.1) ) ) );
    // res = opUnion( res, materialNew( vec3(0.2, 0.2, 0.8), cylinderSDF( pos-vec3( 2.0, 0.50,-2.0), vec2(0.2,0.4) ) ) );
    // res = opUnion( res, materialNew( vec3(0.7, 0.5, 0.2), coneSDF(     pos-vec3( 0.0, 0.75,-2.0), vec3(0.8,0.6,0.6) ) ) );
    // res = opUnion( res, materialNew( vec3(0.4, 0.2, 0.9), hexPrismSDF( pos-vec3(-2.0, 0.60, 2.0), vec2(0.5,0.1) ) ) );
    // res = opUnion( res, materialNew( vec3(0.1, 0.3, 0.6), pyramidSDF(  pos-vec3( 2.0, 0.10, 2.0), 1.0 ) ) );;
    return res;
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 u_resolution = iResolution.xy;
    vec2 u_mouse = iMouse.xy;
    float u_time = iTime;

    vec3 color = vec3(0.0);
    vec2 pixel = 1.0/u_resolution;
    vec2 st = fragCoord.xy * pixel;
    vec2 uv = ratio(st, u_resolution);

    vec2 mo = u_mouse * pixel;
    float time = 32.0 + u_time * 1.5;
    // vec3 cam = vec3( 12.5*cos(0.1*time), 2.2, 12.5*sin(0.1*time) );
    // vec3 cam = vec3( 8.5*cos(0.1*time - 7.0*mo.x), 2.2, 8.5*sin(0.1*time - 7.0*mo.x) );

    const vec2 initMouse = vec2(0.2, -0.1);
    vec2 mouse = initMouse + vec2(0.035625*sin(u_time*PI), 0.0); // initMouse; // iMouse.xy == vec2(0.,0.) ? initMouse : (iMouse.xy/iResolution.xy - 0.5);
    float pitch = clamp(- mouse.x * 2.*PI * 1., -PI,PI);
    float yaw = 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));
    vec3 camPos = 0.8*vec3(-4., 2.2, -2.);
    // apply
    vec3 cam = rot*camPos;

    color = raymarch(cam, vec3(0.0), uv).rgb;
    color = linear2gamma(color);

    // if color is almost white, make it transparent
    // if (color.r > 0.99 && color.g > 0.99 && color.b > 0.99) {
    //     fragColor = vec4(0.0,0.0,0.0,0.0);
    // } else {
        fragColor = vec4( color, 1.0 );
    // }
}

Adapted from ⧉:

#define EDGE_WIDTH 0.12
#define RAYMARCH_ITERATIONS 35
#define SHADOW_ITERATIONS 40
#define SHADOW_STEP 1.0
#define SHADOW_SMOOTHNESS 256.0
#define SHADOW_DARKNESS 0.75

// Distance functions from iquilezles.org
float fSubtraction(float a, float b) {return max(-a,b);}
float fIntersection(float d1, float d2) {return max(d1,d2);}
void fUnion(inout float d1, float d2) {d1 = min(d1,d2);}
float pSphere(vec3 p, float s) {return length(p)-s;}
float pRoundBox(vec3 p, vec3 b, float r) {return length(max(abs(p)-b,0.0))-r;}
float pTorus(vec3 p, vec2 t) {vec2 q = vec2(length(p.xz)-t.x,p.y); return length(q)-t.y;}
float pTorus2(vec3 p, vec2 t) {vec2 q = vec2(length(p.xy)-t.x,p.z); return length(q)-t.y;}
float pCapsule(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 mySDF(vec3 p) {
  float x = p.x;
  float y = p.y;
  float z = p.z;

  float x2 = pow(x, 2.);
  float y2 = pow(y, 2.);
  float z2 = pow(z, 2.);

  float x3 = pow(x, 3.);
  float y3 = pow(y, 3.);
  float z3 = pow(z, 3.);

  float x4 = pow(x, 4.);
  float y4 = pow(y, 4.);
  float z4 = pow(z, 4.);

  // return (y-x2)*(z-x3);

  // return 5.*(z2+y3-y4-x2*y2);

  // return (y2-x2-z2);

  return (x2-z2*y2+y3);
}

float distf(vec3 p)
{
  float d = 100000.0;

  fUnion(d, pRoundBox(vec3(0,0,10) + p, vec3(21,21,1), 1.0));
  fUnion(d, pSphere(vec3(10,10,0) + p, 8.0));
  fUnion(d, pSphere(vec3(16,0,4) + p, 4.0));
  fUnion(d, pCapsule(p, vec3(10,10,12), vec3(15,15,-6.5), 1.5));
  fUnion(d, pCapsule(p, vec3(10,10,12), vec3(5,15,-6.5), 1.5));
  fUnion(d, pCapsule(p, vec3(10,10,12), vec3(10,5,-6.5), 1.5));
  fUnion(d, pTorus(vec3(15,-15,0) + p, vec2(6,2)));
  fUnion(d, pTorus2(vec3(10,-15,0) + p, vec2(6,2)));
  fUnion(d, pRoundBox(vec3(-10,10,-2) + p, vec3(1,1,9), 1.0));
  fUnion(d, pRoundBox(vec3(-10,10,-4) + p, vec3(0.5,6,0.5), 1.0));
  fUnion(d, pRoundBox(vec3(-10,10,2) + p, vec3(6,0.5,0.5), 1.0));

  // d = mySDF(p);
  // d = fIntersection(d, pSphere(p,15.0));

  return d;
}


vec3 grad(vec3 p)
{
  const float eps = 0.01;
  float m;
    vec3 d_distf = vec3( (distf(vec3(p.x-eps,p.y,p.z)) - distf(vec3(p.x+eps,p.y,p.z))),
                   (distf(vec3(p.x,p.y-eps,p.z)) - distf(vec3(p.x,p.y+eps,p.z))),
                   (distf(vec3(p.x,p.y,p.z-eps)) - distf(vec3(p.x,p.y,p.z+eps)))
         );
    return d_distf / (2.*eps);
}

vec3 grad2(vec3 p)
{
  const float eps = 0.01;
  float m;
    vec3 d_grad = vec3( (grad(vec3(p.x-eps,p.y,p.z)) - grad(vec3(p.x+eps,p.y,p.z))).x,
                   (grad(vec3(p.x,p.y-eps,p.z)) - grad(vec3(p.x,p.y+eps,p.z))).y,
                   (grad(vec3(p.x,p.y,p.z-eps)) - grad(vec3(p.x,p.y,p.z+eps)).z)
         );
    return d_grad / (2.*eps);
}

vec3 normal(vec3 p)
{
    return normalize(grad(p));
}


vec4 raymarch(vec3 from, vec3 increment)
{
  const float maxDist = 200.0;
  const float minDist = 0.001;
  const int maxIter = RAYMARCH_ITERATIONS;

  float dist = 0.0;

  float lastDistEval = 1e10;
  float edge = 0.0;

  for(int i = 0; i < maxIter; i++) {
    vec3 pos = (from + increment * dist);
    float distEval = distf(pos);

    if (lastDistEval < EDGE_WIDTH && distEval > lastDistEval + 0.001) {
      edge = 1.0;
    }


    if (distEval < minDist) {
      break;
    }

    dist += distEval;
    if (distEval < lastDistEval) lastDistEval = distEval;
  }

  float mat = 1.0;
  if (dist >= maxDist) mat = 0.0;

  return vec4(dist, mat, edge, 0);
}

float shadow(vec3 from, vec3 increment)
{
  const float minDist = 1.0;

  float res = 1.0;
  float t = 1.0;
  for(int i = 0; i < SHADOW_ITERATIONS; i++) {
    float h = distf(from + increment * t);
    if(h < minDist)
      return 0.0;

      res = min(res, SHADOW_SMOOTHNESS * h / t);
      t += SHADOW_STEP;
    }
    return res;
}

float rand(float x)
{
  return fract(sin(x) * 43758.5453);
}

float triangle(float x)
{
  return abs(1.0 - mod(abs(x), 2.0)) * 2.0 - 1.0;
}

// Camera localized normal
vec3 campos, camup;
vec3 localNormal(vec3 p, vec3 rd) {
    vec3 n = normal(p), ln;
    vec3 side = cross(campos, camup);
    return vec3(dot(n,  side), dot(n,  camup), dot(n,  -rd));
}

float time;
vec4 getPixel(vec2 p, vec3 from, vec3 increment, vec3 light)
{
  vec4 c = raymarch(from, increment);
  vec3 hitPos = from + increment * c.x;
  vec3 normalDir = normal(hitPos);


  float diffuse = 1.0 + min(0.0, dot(normalDir, -light));
  float inshadow = (1.0 - shadow(hitPos, -light)) * SHADOW_DARKNESS;

  // diffuse = max(diffuse, inshadow);

  // if it's not edge, and no contact
  if (c.z != 1.0 && c.y == 0.0) return vec4(0.96, 0.94, 0.87,1);

  float low = 0.05;
  float high = 0.95;

  diffuse = diffuse > high ? 1.0 : (diffuse > low ? low : diffuse);

  vec4 mCol = mix(vec4(vec3(0.78, 0.76, 0.78) * 1.15,1), vec4(0.70, 0.68, 0.75,1), diffuse);

  float gridStep = 1.0;

  // optional chess style grid
  // mCol = mix(mCol,vec4(0.,0.,0.,1.0),0.1*mod(floor(hitPos.x/gridStep)+floor(hitPos.y/gridStep)+floor(hitPos.z/gridStep),2.));

  float dt = dot(vec3(0., 0., 1.), normalDir);

  float eps_high = 0.02;
  float eps_low = 0.00000; //1;

  // vec3 n = vec4(normalDir, 0.).xyz;
  vec3 n = localNormal(hitPos, increment); // grad(hitPos); // localNormal(hitPos, increment);
  vec3 ez = vec3(0.,0.,0.1);

  // https://www.shadertoy.com/view/fsGXzc

  const float MAX_DIST = 200.0;
  float depth = distance(from, hitPos); // /MAX_DIST;

    // I've mostly just copied and pasted Evan's code.
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	// Compute curvature
	vec3 dx = dFdx(n);
	vec3 dy = dFdy(n);
	vec3 xneg = n - dx;
	vec3 xpos = n + dx;
	vec3 yneg = n - dy;
	vec3 ypos = n + dy;
	float curvature = (cross(xneg, xpos).y - cross(yneg, ypos).x) * 4.0 / (depth);

  // curvature debug
  // return vec4(vec3(1,0.9,0.8) * 0.7 * (abs(curvature) * 500.),1);
  ;
  vec4 curve_color = vec4(mix(vec3(1,0.9,0.8) * 0.5, vec3(1.), 0.5), 1.); // * length(grad(hitPos).xy) //  * abs(curvature) * 200.

  float curvature_eps = 0.0001;
  if(c.z != 1. && abs(curvature) > curvature_eps && length(grad(hitPos).z) > eps_low && ( length(grad(hitPos).z) < eps_high * abs(curvature) * 100. ) ) {
      return curve_color;
  }

  if(c.z != 1. && abs(curvature) > curvature_eps && length(grad(hitPos).y) > eps_low && ( length(grad(hitPos).y) < eps_high * abs(curvature) * 100.) ) {
      return curve_color;
  }

  if(c.z != 1. && abs(curvature) > curvature_eps && length(grad(hitPos).x) > eps_low && ( length(grad(hitPos).x) < eps_high * abs(curvature) * 100.) ) {
      return curve_color;
  }

  if(c.z != 1. && length(grad(hitPos).z) > eps_low && ( length(grad(hitPos).z) < eps_high) ) {
      // return vec4(mix(vec3(1,0.9,0.8) * 0.5, vec3(1.), 0.5), 1.); // * (abs(curvature) * 1000.),1);
  }

  // .z is edge
  mCol = mix(mCol,vec4(vec3(1,0.9,0.8) * 0.5,1),c.z);

  return mCol;
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
  time = floor(iTime * 16.0) / 16.0;
  // pixel position
  vec2 q = fragCoord.xy / iResolution.xy;
  vec2 p = -1.0+2.0*q;
  p.x *= -iResolution.x/iResolution.y;

  // mouse
  vec2 mo = iMouse.xy/iResolution.xy;
  vec2 m = iMouse.xy / iResolution.xy;
  if (iMouse.x == 0.0 && iMouse.y == 0.0) {
    m = vec2(time * 0.06 + 1.67, 0.78);
  }
  m = -1.0 + 2.0 * m;
  m *= vec2(4.0,-0.75);
  m.y += 0.75;

  // camera position
  float dist = 65.0;
  vec3 ta = vec3(0,0,0);
  vec3 ro = vec3(cos(m.x) * cos(m.y) * dist, sin(m.x) * cos(m.y) * dist, sin(m.y) * dist);
  vec3 light = vec3(cos(m.x - 2.27) * 50.0, sin(m.x - 2.27) * 50.0, -20.0);

  // camera direction
  vec3 cw = normalize( ta-ro );
  vec3 cp = vec3( 0.0, 0.0, 1.0 );
  vec3 cu = normalize( cross(cw,cp) );
  vec3 cv = normalize( cross(cu,cw) );
  vec3 rd = normalize( p.x*cu + p.y*cv + 2.5*cw );
  campos = -cw, camup = cv;

  // calculate color
  vec4 col = getPixel(p, ro, rd, normalize(light));
  col = pow(col, vec4(1.0 / 2.2));
  col = col * 1.8 - 0.8;
  fragColor = col;
}

Color extracted with the help of ⧉.

See also:

- toon shader
- GLSL Fragment Shader: Sobel Edge Detection
- Cross hatching WebGL shader
- Penumbra Maps: Approximate Soft Shadows in Real-Time
- stackgl/glsl-lighting-walkthrough

Adpated from ⧉:

Rendered by WebGL:

Rendered by SVG:

This is an initial test for Use.GPU:

Adapted from ⧉:

@compute @workgroup_size(16, 16)
fn main_image(@builtin(global_invocation_id) id: vec3u) {
    // Viewport resolution (in pixels)
    let screen_size = textureDimensions(screen);

    // Prevent overdraw for workgroups on the edge of the viewport
    if (id.x >= screen_size.x || id.y >= screen_size.y) { return; }

    // Pixel coordinates (centre of pixel, origin at bottom left)
    let fragCoord = vec2f(f32(id.x) + .5, f32(screen_size.y - id.y) - .5);

    // Normalised pixel coordinates (from 0 to 1)
    let uv = fragCoord / vec2f(screen_size);

    // Time varying pixel colour
    var col = .5 + .5 * cos(time.elapsed + uv.xyx + vec3f(0.,2.,4.));

    // Convert from gamma-encoded to linear colour space
    col = pow(col, vec3f(2.2));

    // Output to screen (linear colour space)
    textureStore(screen, id.xy, vec4f(col, 1.));
}

From ⧉:

From ⧉:

- Styling considerations: Illustrating Geometry by Keenan Crane
- Workflow: Explanatory paper figures with Illustrator and Blender
- TikZ drawing tools
  - Quiver
  - TikZiT
  - TikzMaker
  - awesome-latex-drawing
- TikZ browser-side rendering
  - TikZJax
  - LaTeX.wasm

The following is adapted from exo7math/deepmath-exo7:

The following is adapted from bcharlier/HLMA410:

The following is adapted from giaan-nguyen/latex-notes:

From ⧉:

The following is adapted from 20squares/blog (rendered):

From tikz.org:

From ⧉:

From ⧉:

From ⧉:

From ⧉:

From ⧉:

From ⧉:

From ⧉:

From ⧉:

Reproducing \(\mathbf {V}\left (x^2-y^2 z^2+z^3\right )\) from [cox1997ideals, p. 7, p. 16] using addplot3

Adpated from Eric Wieser's snippet from the bivector discord (is it possible to combine this with kingdon and pylatex to render it with PGA etc.?):

From ⧉:

From ⧉:

The following are now used extensively in related links and backlinks of string diagram [nakahira2023diagrammatic, table 1.1].

The following is testing string diagram drawing using the example from Jon Sterling's String diagrams (extensively used in [sterling2023models] and [sterling2023monads]):

The following follows Algebraic description of natural models

The following is testing string diagram drawing using the example from Category theory using string diagrams:

From ⧉ (doc):

From ⧉:

- Tool: PGFPlotsEdt
- Three-dimensional graphics with PGF/TikZ by Keith Wolcott
- Draw a surface using Tikz on SE

This is native forester markup:

• external link: typst docs
• wiki link: test tikz drawing
• citations: [nakahira2023diagrammatic]
• cross-references: § [uts-000E]

This is the same markup rendered in Typst:

Adapted from ⧉, removed figures and bib for now, due to depending on external non-Typst files.