[WIP] Shadertoy and shader plot types

GLMakie/assets/shader/shadertoy.frag

@@ -0,0 +1,16 @@
+{{GLSL_VERSION}}
+
+{{SHADERTOY_INPUTS}}
+{{TOY_SHADER}}
+
+in vec2 f_uv;
+
+uniform uint objectid;
+
+void write2framebuffer(vec4 color, uvec2 id);
+
+void main()
+{
+    vec4 color = mainImage(f_uv * iResolution.xy);
+    write2framebuffer(color, uvec2(objectid, 0));
+}

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;
+}

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);
+}

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.6743e-11  # Gravitational constant
+const eps = 1.0e-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
+
+    return 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)
+    return @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