diff --git a/GLMakie/assets/shader/shadertoy.frag b/GLMakie/assets/shader/shadertoy.frag
new file mode 100644
index 00000000000..a7afe50a676
--- /dev/null
+++ b/GLMakie/assets/shader/shadertoy.frag
@@ -0,0 +1,15 @@
+{{GLSL_VERSION}}
+
+{{SHADERTOY_INPUTS}}
+{{TOY_SHADER}}
+
+in vec2 f_uv;
+layout (location = 0) out vec4 fragment_color;
+layout (location = 1) out uvec2 fragment_groupid;
+
+void main()
+{
+	vec4 color = mainImage(f_uv * iResolution.xy);
+	if (color.a <= 0.0) discard;
+	fragment_color = color * 0.5;
+}
diff --git a/GLMakie/assets/shader/shadertoy.vert b/GLMakie/assets/shader/shadertoy.vert
new file mode 100644
index 00000000000..e928ede2a5e
--- /dev/null
+++ b/GLMakie/assets/shader/shadertoy.vert
@@ -0,0 +1,11 @@
+{{GLSL_VERSION}}
+
+in vec2 vertices;
+in vec2 uv;
+out vec2 f_uv;
+uniform mat4 projection, view, model;
+
+void main(){
+    gl_Position = projection * view * vec4(vertices, 0, 1);
+    f_uv = uv;
+}
diff --git a/GLMakie/experiments/clouds.frag b/GLMakie/experiments/clouds.frag
new file mode 100644
index 00000000000..32bc3cf6f15
--- /dev/null
+++ b/GLMakie/experiments/clouds.frag
@@ -0,0 +1,122 @@
+float noise(in vec3 x) {
+    vec3 p = floor(x);
+    vec3 f = fract(x);
+    f = f * f * (3.0 - 2.0 * f);
+    vec2 uv = (p.xy + vec2(37.0, 17.0) * p.z) + f.xy;
+    vec2 rg = texture(iChannel0, (uv + 0.5) / 256.0, -100).yx;
+    return -1.0 + 2.0 * mix(rg.x, rg.y, f.z);
+}
+
+float map5(in vec3 p) {
+    vec3 q = p - vec3(0.0, 0.1, 1.0) * iGlobalTime;
+    float f;
+    f = 0.50000 * noise(q);
+    q = q * 2.02;
+    f += 0.25000 * noise(q);
+    q = q * 2.03;
+    f += 0.12500 * noise(q);
+    q = q * 2.01;
+    f += 0.06250 * noise(q);
+    q = q * 2.02;
+    f += 0.03125 * noise(q);
+    return clamp(1.5 - p.y - 2.0 + 1.75 * f, 0.0, 1.0);
+}
+
+float map4(in vec3 p) {
+    vec3 q = p - vec3(0.0, 0.1, 1.0) * iGlobalTime;
+    float f;
+    f = 0.50000 * noise(q);
+    q = q * 2.02;
+    f += 0.25000 * noise(q);
+    q = q * 2.03;
+    f += 0.12500 * noise(q);
+    q = q * 2.01;
+    f += 0.06250 * noise(q);
+    return clamp(1.5 - p.y - 2.0 + 1.75 * f, 0.0, 1.0);
+}
+float map3(in vec3 p) {
+    vec3 q = p - vec3(0.0, 0.1, 1.0) * iGlobalTime;
+    float f;
+    f = 0.50000 * noise(q);
+    q = q * 2.02;
+    f += 0.25000 * noise(q);
+    q = q * 2.03;
+    f += 0.12500 * noise(q);
+    return clamp(1.5 - p.y - 2.0 + 1.75 * f, 0.0, 1.0);
+}
+float map2(in vec3 p) {
+    vec3 q = p - vec3(0.0, 0.1, 1.0) * iGlobalTime;
+    float f;
+    f = 0.50000 * noise(q);
+    q = q * 2.02;
+    f += 0.25000 * noise(q);
+    ;
+    return clamp(1.5 - p.y - 2.0 + 1.75 * f, 0.0, 1.0);
+}
+
+vec3 sundir = normalize(vec3(-1.0, 0.0, -1.0));
+
+vec4 integrate(in vec4 sum, in float dif, in float den, in vec3 bgcol, in float t) {
+    // lighting
+    vec3 lin = vec3(0.65, 0.68, 0.7) * 1.3 + 0.5 * vec3(0.7, 0.5, 0.3) * dif;
+    vec4 col = vec4(mix(1.15 * vec3(1.0, 0.95, 0.8), vec3(0.65), den), den);
+    col.xyz *= lin;
+    col.xyz = mix(col.xyz, bgcol, 1.0 - exp(-0.003 * t * t));
+    // front to back blending
+    col.a *= 0.4;
+    col.rgb *= col.a;
+    return sum + col * (1.0 - sum.a);
+}
+
+#define MARCH(STEPS,MAPLOD) for(int i=0; i<STEPS; i++) { vec3  pos = ro + t*rd; if( pos.y<-3.0 || pos.y>2.0 || sum.a > 0.99 ) break; float den = MAPLOD( pos ); if( den>0.01 ) { float dif =  clamp((den - MAPLOD(pos+0.3*sundir))/0.6, 0.0, 1.0 ); sum = integrate( sum, dif, den, bgcol, t ); } t += max(0.1,0.02*t); }
+
+vec4 raymarch(in vec3 ro, in vec3 rd, in vec3 bgcol) {
+    vec4 sum = vec4(0.0);
+
+    float t = 0.0;
+
+    MARCH(30, map5);
+    MARCH(30, map4);
+    MARCH(30, map3);
+    MARCH(30, map2);
+
+    return clamp(sum, 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 = normalize(cross(cu, cw));
+    return mat3(cu, cv, cw);
+}
+
+vec4 render(in vec3 ro, in vec3 rd) {
+    // background sky
+    float sun = clamp(dot(sundir, rd), 0.0, 1.0);
+    vec3 col = vec3(0.6, 0.71, 0.75) - rd.y * 0.2 * vec3(1.0, 0.5, 1.0) + 0.15 * 0.5;
+    col += 0.2 * vec3(1.0, .6, 0.1) * pow(sun, 8.0);
+
+    // clouds
+    vec4 res = raymarch(ro, rd, col);
+    col = col * (1.0 - res.w) + res.xyz;
+
+    // sun glare
+    col += 0.1 * vec3(1.0, 0.4, 0.2) * pow(sun, 3.0);
+
+    return vec4(col, 1.0);
+}
+
+vec4 mainImage(in vec2 fragCoord) {
+    vec2 p = (-iResolution.xy + 2.0 * gl_FragCoord.xy) / iResolution.y;
+    vec2 m = iMouse.xy / iResolution.xy;
+
+    // camera
+    vec3 ro = 4.0 * normalize(vec3(sin(3.0 * m.x), 0.4 * m.y, cos(3.0 * m.x)));
+    vec3 ta = vec3(0.0, -1.0, 0.0);
+    mat3 ca = setCamera(ro, ta, 0.0);
+    // ray
+    vec3 rd = ca * normalize(vec3(p.xy, 1.5));
+
+    return render(ro, rd);
+}
diff --git a/GLMakie/experiments/gpu-examples.jl b/GLMakie/experiments/gpu-examples.jl
new file mode 100644
index 00000000000..5edce7751fd
--- /dev/null
+++ b/GLMakie/experiments/gpu-examples.jl
@@ -0,0 +1,78 @@
+using CUDA
+
+const N = 1024  # Number of particles
+const dt = 0.01  # Time step
+const G = 6.67430e-11  # Gravitational constant
+const eps = 1e-3  # Softening factor to avoid singularity
+
+
+# Struct to represent a particle
+struct Particle
+    pos::CUDA.CuVector{Point3f}
+    vel::CUDA.CuVector{Vec3f}
+    mass::CUDA.CuVector{Float32}
+end
+
+# Initialize particles with random positions, velocities, and masses
+function init_particles(n)
+    pos = CUDA.fill(Point3f(0.0f0, 0.0f0, 0.0f0), n)
+    vel = CUDA.fill(Vec3f(0.0f0, 0.0f0, 0.0f0), n)
+    mass = CUDA.fill(0.0f0, n)
+
+    for i in 1:n
+        pos[i] = rand(Point3f) .* 100
+        vel[i] = (rand(Vec3f) .- 0.5) .* 10
+        mass[i] = rand(Float32) * 10 + 1
+    end
+
+    Particle(pos, vel, mass)
+end
+
+# GPU kernel to calculate forces and update positions and velocities
+function update_particles!(pos, mass, vel, dt, G, eps, N)
+    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
+    @inbounds if i <= N
+        force = Vec3f(0.0f0, 0.0f0, 0.0f0)
+        pos_i = pos[i]
+
+        for j in 1:N
+            if i != j
+                pos_j = pos[j]
+                diff = pos_j - pos_i
+                distSqr = sum(diff .* diff) + eps * eps
+                invDist = 1.0f0 / sqrt(distSqr)
+                invDist3 = invDist * invDist * invDist
+                mass_product = mass[i] * mass[j]
+                force += diff .* (G * mass_product * invDist3)
+            end
+        end
+
+        accel = force ./ mass[i]
+        vel[i] += accel .* dt
+        pos[i] += vel[i] .* dt
+    end
+    return
+end
+
+# Main simulation function
+function run_simulation(particles, dt, G, eps, N)
+    threads = 256
+    blocks = ceil(Int, N / threads)
+    @cuda threads = threads blocks = blocks update_particles!(particles, dt, G, eps, N)
+end
+
+# Initialize particles
+particles = init_particles(N)
+
+# Number of simulation steps
+steps = 100
+
+# Run the simulation
+run_simulation(particles, dt, steps, G, eps, N)
+
+# Fetch the results from the GPU
+positions = Array(pos)
+println("Final positions of particles:")
+for pos in positions
+    println(pos)
+end
diff --git a/GLMakie/experiments/monster.frag b/GLMakie/experiments/monster.frag
new file mode 100644
index 00000000000..6ef416c6c02
--- /dev/null
+++ b/GLMakie/experiments/monster.frag
@@ -0,0 +1,217 @@
+// Created by inigo quilez - iq/2013
+// License Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.
+
+vec3 hash3( float n )
+{
+    return fract(sin(vec3(n,n+1.0,n+2.0))*vec3(13.5453123,31.1459123,37.3490423));
+}
+
+vec3 noise( in float x )
+{
+    float p = floor(x);
+    float f = fract(x);
+    f = f*f*(3.0-2.0*f);
+    return mix( hash3(p+0.0), hash3(p+1.0),f);
+}
+
+
+mat4 rotationMat( in vec3 xyz )
+{
+    vec3 si = sin(xyz);
+    vec3 co = cos(xyz);
+
+    return mat4( co.y*co.z,                co.y*si.z,               -si.y,       0.0,
+                 si.x*si.y*co.z-co.x*si.z, si.x*si.y*si.z+co.x*co.z, si.x*co.y,  0.0,
+                 co.x*si.y*co.z+si.x*si.z, co.x*si.y*si.z-si.x*co.z, co.x*co.y,  0.0,
+                 0.0,                      0.0,                      0.0,        1.0 );
+}
+
+const float s = 1.1;
+
+mat4 mm;
+
+vec3 map( vec3 p )
+{
+    float k = 1.0;
+    float m = 1e10;
+    for( int i=0; i<22; i++ )
+    {
+        m = min( m, dot(p,p)/(k*k) );
+        p = (mm*vec4((abs(p)),1.0)).xyz;
+        k*= s;
+    }
+
+
+    float d = (length(p)-0.25)/k;
+
+    float h = p.z - 0.35*p.x;
+
+    return vec3( d, m, h );
+}
+
+vec3 intersect( in vec3 ro, in vec3 rd )
+{
+    float t = 0.0;
+    for( int i=0; i<100; i++ )
+    {
+        vec3 res = map( ro+rd*t );
+        if( res.x<0.0002 ) return vec3(t,res.yz);
+        t += res.x;
+        if( t>9.0 ) break;
+    }
+
+    return vec3( -1.0 );
+}
+
+vec3 calcNormal( in vec3 pos, float e )
+{
+    vec3 eps = vec3(e,0.0,0.0);
+
+    return normalize( vec3(
+           map(pos+eps.xyy).x - map(pos-eps.xyy).x,
+           map(pos+eps.yxy).x - map(pos-eps.yxy).x,
+           map(pos+eps.yyx).x - map(pos-eps.yyx).x ) );
+}
+
+float softshadow( in vec3 ro, in vec3 rd, float mint, float k )
+{
+    float res = 1.0;
+    float t = mint;
+    for( int i=0; i<32; i++ )
+    {
+        float h = map(ro + rd*t).x;
+        h = max( h, 0.0 );
+        res = min( res, k*h/t );
+        t += clamp( h, 0.001, 0.1 );
+        if( res<0.01 || t>6.0 ) break;
+    }
+    return clamp(res,0.0,1.0);
+}
+
+float calcAO( in vec3 pos, in vec3 nor )
+{
+    float totao = 0.0;
+    for( int aoi=0; aoi<16; aoi++ )
+    {
+        vec3 aopos = -1.0+2.0*hash3(float(aoi)*213.47);
+        aopos *= sign( dot(aopos,nor) );
+        aopos = pos + nor*0.01 + aopos*0.04;
+        float dd = clamp( map( aopos ).x*4.0, 0.0, 1.0 );
+        totao += dd;
+    }
+    totao /= 16.0;
+
+    return clamp( totao*totao*50.0, 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 = normalize( cross(cu,cw) );
+    return mat3( cu, cv, cw );
+}
+
+vec4 mainImage(in vec2 fragCoord )
+{
+    vec2 q = fragCoord.xy / iResolution.xy;
+    vec2 p = -1.0 + 2.0 * q;
+    p.x *= iResolution.x/iResolution.y;
+    vec2 m = vec2(0.5);
+    m = iMouse.xy/iResolution.xy;
+
+    // animation
+    float time = iGlobalTime;
+    time += 15.0*smoothstep(  15.0, 25.0, iGlobalTime );
+    time += 20.0*smoothstep(  65.0, 80.0, iGlobalTime );
+    time += 35.0*smoothstep( 105.0, 135.0, iGlobalTime );
+    time += 20.0*smoothstep( 165.0, 180.0, iGlobalTime );
+    time += 40.0*smoothstep( 220.0, 290.0, iGlobalTime );
+    time +=  5.0*smoothstep( 320.0, 330.0, iGlobalTime );
+    float time1 = (time-10.0)*1.5 - 167.0;
+    float time2 = time;
+
+    mm = rotationMat( vec3(0.4,0.1,3.4) +
+                      0.15*sin(0.1*vec3(0.40,0.30,0.61)*time1) +
+                      0.15*sin(0.1*vec3(0.11,0.53,0.48)*time1));
+    mm[0].xyz *= s;
+    mm[1].xyz *= s;
+    mm[2].xyz *= s;
+    mm[3].xyz = vec3( 0.15, 0.05, -0.07 ) + 0.05*sin(vec3(0.0,1.0,2.0) + 0.2*vec3(0.31,0.24,0.42)*time1);
+
+    // camera
+    float an = 1.0 + 0.1*time2 - 6.2*m.x;
+    float cr = 0.15*sin(0.2*time2);
+    vec3 ro = (2.4 + 0.6*smoothstep(10.0,20.0,time2))*vec3(sin(an),0.25,cos(an));
+    vec3 ta = vec3( 0.0, 0.0 + 0.13*cos(0.3*time2), 0.0 );
+    ta += 0.05*noise(  0.0 + 1.0*time );
+    ro += 0.05*noise( 11.3 + 1.0*time );
+    // camera-to-world transformation
+    mat3 ca = setCamera( ro, ta, cr );
+    // ray direction
+    vec3 rd = ca * normalize( vec3(p.xy,3.0) );
+
+    // raymarch
+    vec3 tmat = intersect(ro,rd);
+
+    // shade
+    vec3 col = vec3(0.0);
+    vec4 bg = vec4(0.7,0.5,0.7, 1.0);
+    if( tmat.z>-0.5 )
+    {
+        // geometry
+        vec3 pos = ro + tmat.x*rd;
+        vec3 nor = calcNormal(pos, 0.005);
+        vec3 sor = calcNormal(pos, 0.010);
+
+        // material
+        vec3 mate = vec3(1.0);
+        mate = mix( vec3(0.5,0.5,0.2), vec3(0.5,0.3,0.0), 0.5 + 0.5*sin(4.0+8000.0*tmat.y)  );
+        mate = mix( vec3(1.0,0.9,0.8), mate, 0.5 + 0.5*sin(4.0+20.0*tmat.z) );
+        mate.x *= 1.15;
+
+        // lighting
+        float occ = 1.1*calcAO( pos, nor );
+        occ *= 0.75 + 0.25*clamp(tmat.y*400.0,0.0,1.0);
+
+        // diffuse
+        col = vec3(0.0);
+        for( int i=0; i<32; i++ )
+        {
+            //vec3 rr = normalize(-1.0 + 2.0*texture2D( iChannel2, vec2((0.5+float(i)),0.5)/256.0,-100.0).xyz);
+            vec3 rr = normalize(-1.0 + 2.0*hash3(float(i)*123.5463));
+            rr = normalize( nor + 7.0*rr );
+            rr = rr * sign(dot(nor,rr));
+            float ds = occ;//softshadow( pos, rr, 0.01, 32.0 );
+            col += pow( bg.xyz, vec3(2.2) ) * dot(rr,nor) * ds;
+        }
+        col /= 32.0;
+
+        col *= 1.8;
+
+        // subsurface
+        col *= 1.0 + 1.0*vec3(1.0,0.6,0.1)*pow(clamp(1.0+dot(rd,sor),0.0,1.0),2.0)*vec3(1.0);
+
+        // specular
+        float fre = pow( clamp(1.0+dot(rd,nor),0.0,1.0), 5.0 );
+        vec3 ref = reflect( rd, nor );
+        float rs = softshadow( pos, ref, 0.01, 32.0 );
+        col += 1.8 * (0.04 + 12.0*fre) * occ * pow( bg.xyz, vec3(2.0) ) * rs;
+
+        col *= mate;
+    }
+    else
+    {
+        // background
+        discard;
+    }
+
+    // gamma
+    col = pow( clamp( col, 0.0, 1.0 ), vec3(0.45) );
+
+    // vigneting
+    col *= 0.5 + 0.5*pow( 16.0*q.x*q.y*(1.0-q.x)*(1.0-q.y), 0.1 );
+
+    return vec4( col, 1.0 );
+}
diff --git a/GLMakie/experiments/noise.png b/GLMakie/experiments/noise.png
new file mode 100644
index 00000000000..32a561c7f68
Binary files /dev/null and b/GLMakie/experiments/noise.png differ
diff --git a/GLMakie/experiments/raytracing.frag b/GLMakie/experiments/raytracing.frag
new file mode 100644
index 00000000000..266d1feaf7f
--- /dev/null
+++ b/GLMakie/experiments/raytracing.frag
@@ -0,0 +1,298 @@
+// Created by inigo quilez - iq/2013
+// License Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.
+
+// A list of usefull distance function to simple primitives, and an example on how to
+// do some interesting boolean operations, repetition and displacement.
+//
+// More info here: http://www.iquilezles.org/www/articles/distfunctions/distfunctions.htm
+
+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 udRoundBox(vec3 p, vec3 b, float r) {
+return length(max(abs(p) - b, 0.0)) - r;
+}
+
+float sdTorus(vec3 p, vec2 t) {
+return length(vec2(length(p.xz) - t.x, p.y)) - t.y;
+}
+
+float sdHexPrism(vec3 p, vec2 h) {
+vec3 q = abs(p);
+#if 0
+return max(q.z - h.y, max((q.x * 0.866025 + q.y * 0.5), q.y) - h.x);
+#else
+float d1 = q.z - h.y;
+float d2 = max((q.x * 0.866025 + q.y * 0.5), q.y) - h.x;
+return length(max(vec2(d1, d2), 0.0)) + min(max(d1, d2), 0.);
+#endif
+}
+
+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 sdTriPrism(vec3 p, vec2 h) {
+vec3 q = abs(p);
+#if 0
+return max(q.z - h.y, max(q.x * 0.866025 + p.y * 0.5, - p.y) - h.x * 0.5);
+#else
+float d1 = q.z - h.y;
+float d2 = max(q.x * 0.866025 + p.y * 0.5, - p.y) - h.x * 0.5;
+return length(max(vec2(d1, d2), 0.0)) + min(max(d1, d2), 0.);
+#endif
+}
+
+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));
+}
+
+float sdCone(in vec3 p, in vec3 c) {
+vec2 q = vec2(length(p.xz), p.y);
+#if 0
+return max(max(dot(q, c.xy), p.y), - p.y - c.z);
+#else
+float d1 = - p.y - c.z;
+float d2 = max(dot(q, c.xy), p.y);
+return length(max(vec2(d1, d2), 0.0)) + min(max(d1, d2), 0.);
+#endif
+}
+
+float length2(vec2 p) {
+return sqrt(p.x * p.x + p.y * p.y);
+}
+
+float length6(vec2 p) {
+p = p * p * p;
+p = p * p;
+return pow(p.x + p.y, 1.0 / 6.0);
+}
+
+float length8(vec2 p) {
+p = p * p;
+p = p * p;
+p = p * p;
+return pow(p.x + p.y, 1.0 / 8.0);
+}
+
+float sdTorus82(vec3 p, vec2 t) {
+vec2 q = vec2(length2(p.xz) - t.x, p.y);
+return length8(q) - t.y;
+}
+
+float sdTorus88(vec3 p, vec2 t) {
+vec2 q = vec2(length8(p.xz) - t.x, p.y);
+return length8(q) - t.y;
+}
+
+float sdCylinder6(vec3 p, vec2 h) {
+return max(length6(p.xz) - h.x, abs(p.y) - h.y);
+}
+
+//----------------------------------------------------------------------
+
+float opS(float d1, float d2) {
+return max(- d2, d1);
+}
+
+vec2 opU(vec2 d1, vec2 d2) {
+return (d1.x < d2.x) ? d1 : d2;
+}
+
+vec3 opRep(vec3 p, vec3 c) {
+return mod(p, c) - 0.5 * c;
+}
+
+vec3 opTwist(vec3 p) {
+float c = cos(10.0 * p.y + 10.0);
+float s = sin(10.0 * p.y + 10.0);
+mat2 m = mat2(c, - s, s, c);
+return vec3(m * p.xz, p.y);
+}
+
+//----------------------------------------------------------------------
+
+vec2 map(in vec3 pos) {
+vec2 res = opU(vec2(sdPlane(pos), 1.0), vec2(sdSphere(pos - vec3(0.0, 0.25, 0.0), 0.25), 46.9));
+res = opU(res, vec2(sdBox(pos - vec3(1.0, 0.25, 0.0), vec3(0.25)), 3.0));
+res = opU(res, vec2(udRoundBox(pos - vec3(1.0, 0.25, 1.0), vec3(0.15), 0.1), 41.0));
+res = opU(res, vec2(sdTorus(pos - vec3(0.0, 0.25, 1.0), vec2(0.20, 0.05)), 25.0));
+res = opU(res, vec2(sdCapsule(pos, vec3(- 1.3, 0.20, - 0.1), vec3(- 1.0, 0.20, 0.2), 0.1), 31.9));
+res = opU(res, vec2(sdTriPrism(pos - vec3(- 1.0, 0.25, - 1.0), vec2(0.25, 0.05)), 43.5));
+res = opU(res, vec2(sdCylinder(pos - vec3(1.0, 0.30, - 1.0), vec2(0.1, 0.2)), 8.0));
+res = opU(res, vec2(sdCone(pos - vec3(0.0, 0.50, - 1.0), vec3(0.8, 0.6, 0.3)), 55.0));
+res = opU(res, vec2(sdTorus82(pos - vec3(0.0, 0.25, 2.0), vec2(0.20, 0.05)), 50.0));
+res = opU(res, vec2(sdTorus88(pos - vec3(- 1.0, 0.25, 2.0), vec2(0.20, 0.05)), 43.0));
+res = opU(res, vec2(sdCylinder6(pos - vec3(1.0, 0.30, 2.0), vec2(0.1, 0.2)), 12.0));
+res = opU(res, vec2(sdHexPrism(pos - vec3(- 1.0, 0.20, 1.0), vec2(0.25, 0.05)), 17.0));
+
+res = opU(res, vec2(opS(udRoundBox(pos - vec3(- 2.0, 0.2, 1.0), vec3(0.15), 0.05), sdSphere(pos - vec3(- 2.0, 0.2, 1.0), 0.25)), 13.0));
+res = opU(res, vec2(opS(sdTorus82(pos - vec3(- 2.0, 0.2, 0.0), vec2(0.20, 0.1)), sdCylinder(opRep(vec3(atan(pos.x + 2.0, pos.z) / 6.2831, pos.y, 0.02 + 0.5 * length(pos - vec3(- 2.0, 0.2, 0.0))), vec3(0.05, 1.0, 0.05)), vec2(0.02, 0.6))), 51.0));
+res = opU(res, vec2(0.7 * sdSphere(pos - vec3(- 2.0, 0.25, - 1.0), 0.2) +
+    0.03 * sin(50.0 * pos.x) * sin(50.0 * pos.y) * sin(50.0 * pos.z), 65.0));
+res = opU(res, vec2(0.5 * sdTorus(opTwist(pos - vec3(- 2.0, 0.25, 2.0)), vec2(0.20, 0.05)), 46.7));
+
+return res;
+}
+
+vec2 castRay(in vec3 ro, in vec3 rd) {
+float tmin = 1.0;
+float tmax = 20.0;
+
+#if 0
+float tp1 = (0.0 - ro.y) / rd.y;
+if (tp1 > 0.0) tmax = min(tmax, tp1);
+float tp2 = (1.6 - ro.y) / rd.y;
+if (tp2 > 0.0) {
+if (ro.y > 1.6) tmin = max(tmin, tp2);
+else tmax = min(tmax, tp2);
+}
+#endif
+
+float precis = 0.002;
+float t = tmin;
+float m = - 1.0;
+for (int i = 0;
+i < 50;
+i ++) {
+vec2 res = map(ro + rd * t);
+if (res.x < precis || t > tmax) break;
+t += res.x;
+m = res.y;
+}
+
+if (t > tmax) m = - 1.0;
+return vec2(t, m);
+}
+
+float softshadow(in vec3 ro, in vec3 rd, in float mint, in float tmax) {
+float res = 1.0;
+float t = mint;
+for (int i = 0;
+i < 16;
+i ++) {
+float h = map(ro + rd * t).x;
+res = min(res, 8.0 * h / t);
+t += clamp(h, 0.02, 0.10);
+if (h < 0.001 || t > tmax) break;
+}
+return clamp(res, 0.0, 1.0);
+
+}
+
+vec3 calcNormal(in vec3 pos) {
+vec3 eps = vec3(0.001, 0.0, 0.0);
+vec3 nor = vec3(map(pos + eps.xyy).x - map(pos - eps.xyy).x, map(pos + eps.yxy).x - map(pos - eps.yxy).x, map(pos + eps.yyx).x - map(pos - eps.yyx).x);
+return normalize(nor);
+}
+
+float calcAO(in vec3 pos, in vec3 nor) {
+float occ = 0.0;
+float sca = 1.0;
+for (int i = 0;
+i < 5;
+i ++) {
+float hr = 0.01 + 0.12 * float(i) / 4.0;
+vec3 aopos = nor * hr + pos;
+float dd = map(aopos).x;
+occ += - (dd - hr) * sca;
+sca *= 0.95;
+}
+return clamp(1.0 - 3.0 * occ, 0.0, 1.0);
+}
+
+vec3 render(in vec3 ro, in vec3 rd) {
+vec3 col = vec3(0.8, 0.9, 1.0);
+vec2 res = castRay(ro, rd);
+float t = res.x;
+float m = res.y;
+if (m > - 0.5) {
+vec3 pos = ro + t * rd;
+vec3 nor = calcNormal(pos);
+vec3 ref = reflect(rd, nor);
+
+        // material
+col = 0.45 + 0.3 * sin(vec3(0.05, 0.08, 0.10) * (m - 1.0));
+
+if (m < 1.5) {
+
+float f = mod(floor(5.0 * pos.z) + floor(5.0 * pos.x), 2.0);
+col = 0.4 + 0.1 * f * vec3(1.0);
+}
+
+        // lighitng
+float occ = calcAO(pos, nor);
+vec3 lig = normalize(vec3(- 0.6, 0.7, - 0.5));
+float amb = clamp(0.5 + 0.5 * nor.y, 0.0, 1.0);
+float dif = clamp(dot(nor, lig), 0.0, 1.0);
+float bac = clamp(dot(nor, normalize(vec3(- lig.x, 0.0, - lig.z))), 0.0, 1.0) * clamp(1.0 - pos.y, 0.0, 1.0);
+float dom = smoothstep(- 0.1, 0.1, ref.y);
+float fre = pow(clamp(1.0 + dot(nor, rd), 0.0, 1.0), 2.0);
+float spe = pow(clamp(dot(ref, lig), 0.0, 1.0), 16.0);
+
+dif *= softshadow(pos, lig, 0.02, 2.5);
+dom *= softshadow(pos, ref, 0.02, 2.5);
+
+vec3 brdf = vec3(0.0);
+brdf += 1.20 * dif * vec3(1.00, 0.90, 0.60);
+brdf += 1.20 * spe * vec3(1.00, 0.90, 0.60) * dif;
+brdf += 0.30 * amb * vec3(0.50, 0.70, 1.00) * occ;
+brdf += 0.40 * dom * vec3(0.50, 0.70, 1.00) * occ;
+brdf += 0.30 * bac * vec3(0.25, 0.25, 0.25) * occ;
+brdf += 0.40 * fre * vec3(1.00, 1.00, 1.00) * occ;
+brdf += 0.02;
+col = col * brdf;
+
+col = mix(col, vec3(0.8, 0.9, 1.0), 1.0 - exp(- 0.0005 * 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 = normalize(cross(cu, cw));
+return mat3(cu, cv, cw);
+}
+
+vec4 mainImage(in vec2 fragCoord) {
+vec2 q = gl_FragCoord.xy / iResolution.xy;
+vec2 p = - 1.0 + 2.0 * q;
+p.x *= iResolution.x / iResolution.y;
+vec2 mo = iMouse.xy / iResolution.xy;
+
+float time = 15.0 + iGlobalTime;
+
+        // camera
+vec3 ro = vec3(- 0.5 + 3.2 * cos(0.1 * time + 6.0 * mo.x), 1.0 + 2.0 * mo.y, 0.5 + 3.2 * sin(0.1 * time + 6.0 * mo.x));
+vec3 ta = vec3(- 0.5, - 0.4, 0.5);
+
+        // camera-to-world transformation
+mat3 ca = setCamera(ro, ta, 0.0);
+
+        // ray direction
+vec3 rd = ca * normalize(vec3(p.xy, 2.5));
+
+        // render
+vec3 col = render(ro, rd);
+
+col = pow(col, vec3(0.4545));
+
+        return vec4(col, 1.0);
+}
diff --git a/GLMakie/experiments/shadertoy.jl b/GLMakie/experiments/shadertoy.jl
new file mode 100644
index 00000000000..1f83f20d151
--- /dev/null
+++ b/GLMakie/experiments/shadertoy.jl
@@ -0,0 +1,29 @@
+using GLMakie, GeometryBasics, LinearAlgebra, FileIO
+GLMakie.activate!(float=true, render_on_demand=false, vsync=true)
+frag = read(joinpath(@__DIR__, "clouds.frag"), String)
+img = load(joinpath(@__DIR__, "noise.png"))
+begin
+    GLMakie.closeall()
+    s=Scene()
+    shadertoy!(s,
+        Rect2f(-1, -1, 2, 2), frag;
+        uniforms=Dict{Symbol, Any}(
+            :iChannel0 => GLMakie.Sampler(img; x_repeat=:repeat, minfilter=:linear)
+        )
+    )
+    s
+end
+
+begin
+    GLMakie.closeall()
+    s = Scene()
+    shadertoy!(s, Rect2f(-1, -1, 2, 2), read(joinpath(@__DIR__, "raytracing.frag"), String))
+    scren = display(s)
+end
+
+begin
+    f, ax, pl = shadertoy(Rect2f(-1, -1, 2, 2), read(joinpath(@__DIR__, "raytracing.frag"), String))
+    foreach(n->deregister_interaction!(ax, n), keys(interactions(ax)))
+    f
+end
+shadertoy(Rect2f(-1, -1, 2, 2), read(joinpath(@__DIR__, "monster.frag"), String))
diff --git a/GLMakie/src/GLAbstraction/GLRender.jl b/GLMakie/src/GLAbstraction/GLRender.jl
index 146eeb6dd44..e1ca88dc80f 100644
--- a/GLMakie/src/GLAbstraction/GLRender.jl
+++ b/GLMakie/src/GLAbstraction/GLRender.jl
@@ -85,7 +85,7 @@ function render(renderobject::RenderObject, vertexarray=renderobject.vertexarray
                         error("Uniform tuple too long: $(length(value))")
                     end
                 catch e
-                    @warn error("uniform $key doesn't work with value $(renderobject.uniforms[key])") exception=(e, Base.catch_backtrace())
+                    @warn error("uniform $key doesn't work with value $(typeof(renderobject.uniforms[key]))") exception=(e, Base.catch_backtrace())
                 end
             end
         end
diff --git a/GLMakie/src/GLAbstraction/GLUtils.jl b/GLMakie/src/GLAbstraction/GLUtils.jl
index 4176bdc97b9..28884d301a1 100644
--- a/GLMakie/src/GLAbstraction/GLUtils.jl
+++ b/GLMakie/src/GLAbstraction/GLUtils.jl
@@ -75,7 +75,7 @@ macro gen_defaults!(dict, args)
                         doc_strings = :(doc_strings[$key_sym] = $elem)
                     elseif isa(elem, Symbol)
                         opengl_convert_target = quote
-                            if GLAbstraction.matches_target($elem, $key_name)
+                            if $(GLAbstraction).matches_target($elem, $key_name)
                                 gl_convert_targets[$key_sym] = $elem
                             end
                         end
diff --git a/GLMakie/src/drawing_primitives.jl b/GLMakie/src/drawing_primitives.jl
index f364b2e3468..bee8c21df5b 100644
--- a/GLMakie/src/drawing_primitives.jl
+++ b/GLMakie/src/drawing_primitives.jl
@@ -336,7 +336,11 @@ function cached_robj!(robj_func, screen, scene, plot::AbstractPlot)
         end
         robj = robj_func(gl_attributes)
 
-        get!(gl_attributes, :ssao, Observable(false))
+        gl_attributes[:ssao] = plot.ssao
+        gl_attributes[:transparency] = plot.transparency
+        gl_attributes[:overdraw] = plot.overdraw
+        gl_attributes[:visible] = plot.visible
+        gl_attributes[:fxaa] = plot.fxaa
         screen.cache2plot[robj.id] = plot
         return robj
     end
@@ -1018,3 +1022,53 @@ function draw_atomic(screen::Screen, scene::Scene, plot::Voxels)
         return draw_voxels(screen, tex, gl_attributes)
     end
 end
+
+function draw_atomic(screen::Screen, scene::Scene, plot::Makie.ShaderToy)
+    robj = get!(screen.cache, objectid(plot)) do
+        uniforms = copy(plot.uniforms[])
+        uniforms[:model] = plot.model
+        uniforms[:px_per_unit] = screen.px_per_unit
+        connect_camera!(plot, uniforms, scene.camera, :data)
+
+        m = plot.rect[]
+        uniforms[:faces] = indexbuffer(decompose(GLTriangleFace, m))
+        uniforms[:vertices] = GLBuffer(decompose(Point2f, m))
+        uniforms[:uv] = GLBuffer(GeometryBasics.decompose_uv(m))
+
+        templates = Dict(
+            "SHADERTOY_INPUTS" => """
+            uniform vec2 iResolution;
+            uniform vec2 iMouse;
+            uniform float iGlobalTime;
+            uniform sampler2D iChannel0;
+            uniform sampler2D iChannel1;
+            uniform sampler2D iChannel2;
+            uniform sampler2D iChannel3;
+            """,
+            "TOY_SHADER" => plot.shader[],
+            "buffers" => output_buffers(screen, to_value(false)),
+            "buffer_writes" => output_buffer_writes(screen, to_value(false))
+        )
+        @gen_defaults! uniforms begin
+            iResolution = scene.camera.resolution
+            iChannel0 = nothing => Texture
+            iChannel1 = nothing => Texture
+            iChannel2 = nothing => Texture
+            iChannel3 = nothing => Texture
+            iGlobalTime = lift(x -> Float32(x.time), scene.events.tick)
+            iMouse = lift(Vec2f, scene.events.mouseposition)
+            iGlobalTime = 0.0f0
+            shader = GLVisualizeShader(
+                screen,
+                "fragment_output.frag", "shadertoy.frag", "shadertoy.vert";
+                view=templates)
+        end
+        get!(uniforms, :ssao, Observable(false))
+        get!(uniforms, :transparency, Observable(true))
+        robj = GLMakie.assemble_shader(uniforms)
+        screen.cache2plot[robj.id] = plot
+        return robj
+    end
+    push!(screen, scene, robj)
+    return robj
+end
diff --git a/MakieCore/src/basic_plots.jl b/MakieCore/src/basic_plots.jl
index 35a0c5b0aa2..7a06ef43a49 100644
--- a/MakieCore/src/basic_plots.jl
+++ b/MakieCore/src/basic_plots.jl
@@ -777,3 +777,14 @@ or other array-like output.
     """Sets the color of the arrow head. Will copy `color` if set to `automatic`."""
     arrowcolor = automatic
 end
+
+@recipe ShaderToy (rect::Rect2f, shader::String) begin
+    """Defines code that replaces the specified keys in shaders.
+    E.g. `{{my_code}}` will be replaced by `shader_injections["my_code"]`.
+    These code injections happen before compile time and require
+    `GLMakie.closeall()` to be overwritten.
+    """
+    uniforms = Dict{Symbol,Any}()
+    mixin_generic_plot_attributes()...
+    fxaa = false
+end
diff --git a/src/Makie.jl b/src/Makie.jl
index 6d6a6452068..47f8aab23d9 100644
--- a/src/Makie.jl
+++ b/src/Makie.jl
@@ -153,6 +153,7 @@ include("camera/old_camera3d.jl")
 
 # basic recipes
 include("basic_recipes/convenience_functions.jl")
+include("basic_recipes/shaders.jl")
 include("basic_recipes/ablines.jl")
 include("basic_recipes/annotations.jl")
 include("basic_recipes/arc.jl")
@@ -381,7 +382,8 @@ include("deprecated.jl")
 export Arrows  , Heatmap  , Image  , Lines  , LineSegments  , Mesh  , MeshScatter  , Poly  , Scatter  , Surface  , Text  , Volume  , Wireframe, Voxels
 export arrows  , heatmap  , image  , lines  , linesegments  , mesh  , meshscatter  , poly  , scatter  , surface  , text  , volume  , wireframe, voxels
 export arrows! , heatmap! , image! , lines! , linesegments! , mesh! , meshscatter! , poly! , scatter! , surface! , text! , volume! , wireframe!, voxels!
-
+import MakieCore: ShaderToy, shadertoy, shadertoy!
+export ShaderToy, shadertoy, shadertoy!
 export AmbientLight, PointLight, DirectionalLight, SpotLight, EnvironmentLight, RectLight, SSAO
 export FastPixel
 
diff --git a/src/basic_recipes/shaders.jl b/src/basic_recipes/shaders.jl
new file mode 100644
index 00000000000..8712c050805
--- /dev/null
+++ b/src/basic_recipes/shaders.jl
@@ -0,0 +1,3 @@
+Makie.plot!(plot::ShaderToy) = plot
+point_iterator(p::ShaderToy) = decompose(Point3f, p.rect[])
+convert_arguments(::Type{<:ShaderToy}, args...) = args