Raymarching Toolbox Thread

category: code [glöplog]

There are a lot of separate Raymarching threads around and its a lot to dig through for code. So heres a thread for collecting all the interesting stuff.

First the basics:
Add object to object
min(a,b)
Subtract object b from object a
max(a,-b)
Find intersection of 2 objects
max(a,b)

Distance to sphere
float sphere(float x, float y, float z, float r)
{
return sqrt(x * x + y * y + z * z) - r;
}

Distance to cube
float cube(float x,float y,float z,float size){
return max(max(abs(x)-size,abs(y)-size),abs(z)-size);
}

Signed distance to cube
float signedDistToBox( in vec3 p, in vec3 b )
{
const vec3 di = abs(p) - b;
return min( maxcomp(di), length(max(di,0.0)) );
}

IQ's small jump dist march
pos += 0.25 * dir * dist;

Calculate Normal
float e = …;

vec3 n = vec3(ƒ(p + {e, 0, 0}) - ƒ(p - {e, 0, 0}), ƒ(p + {0, e, 0}) - ƒ(p - {0, e, 0}), ƒ(p + {0, 0, e}) - ƒ(p - {0, 0, e}));
n = normalize(n);

Decipher of Youth Uprisings AO algorithm
for (float i = 0.; i < 5.; ++i) {
colour -= vec3(i * k - ƒ(p + n * i * k)) / pow(2., i);
}

Unk of Quite's awesome extrusion field code
for(float i=0; i<4; i++) {
float3 q=1+i*i*.18*(1+4*(1+.3*sin(t*.001))*sin(float3(5.7,6.4,7.3)*i*1.145+.3*sin(h.w*. 015)*(3+i))),g=(frac(p*q)-.5)/q;
d1=min(d1+.03,max(d1,max(abs(g.x),max(abs(g.y),abs(g.z)))-.148));
}

Full shader sources
http://pastebin.com/9PRG1PfN - Cdak by Quite
http://pastebin.com/41hKL0qa - Slisesix by RGBA

added on the 2011-02-08 23:14:47 by Mewler

Yes, the way to solve fragmentation is to add ANOTHER thread...

added on the 2011-02-08 23:16:12 by kusma

You accidentally

added on the 2011-02-08 23:22:17 by trc_wm

Rotations...

Code:


void rX(inout vec3 p, float a) {
	float c,s;vec3 q=p;
	c = cos(a); s = sin(a);
	p.y = c * q.y - s * q.z;
	p.z = s * q.y + c * q.z;
}

void rY(inout vec3 p, float a) {
	float c,s;vec3 q=p;
	c = cos(a); s = sin(a);
	p.x = c * q.x + s * q.z;
	p.z = -s * q.x + c * q.z;
}

void rZ(inout vec3 p, float a) {
	float c,s;vec3 q=p;
	c = cos(a); s = sin(a);
	p.x = c * q.x - s * q.y;
	p.y = s * q.x + c * q.y;
}
void rXCS(inout vec3 p, float c, float s) {
	vec3 q=p;
	p.y = c * q.y - s * q.z;
	p.z = s * q.y + c * q.z;
}


void rYCS(inout vec3 p, float c, float s) {
	vec3 q=p;
	p.x = c * q.x + s * q.z;
	p.z = -s * q.x + c * q.z;
}

void rZCS(inout vec3 p, float c, float s) {
	vec3 q=p;
	p.x = c * q.x - s * q.y;
	p.y = s * q.x + c * q.y;
}

modified hsv (original idea unc)

Code:



vec3 hsv(float h,float s,float v) {
	return mix(vec3(1.),clamp((abs(fract(h+vec3(3.,2.,1.)/3.)*6.-3.)-1.),0.,1.),s)*v;
}

AO (d is the delta also known as "k" in deciphers version, i the number of steps)

Code:


float ao(vec3 p, vec3 n, float d, float i) {
	float o;
	for (o=1.;i>0.;i--) {
		o-=(i*d-abs(f(p+n*i*d)))/pow(2.,i);
	}
	return o;
}

You can do fake SSS with that AO: instead of ambientOcc = ao(p,n) try sss = ao(p,d) where d is the ray direction

Some primitives...

Code:


//rT radius of the torus, rT radius of the ring iirc.
float torus(vec3 p,float rT,float rR) {
	return length(vec2(length(p.xz) - rT,p.y)) - rR;
}

//n normal of the plane
float plane(vec3 p, vec3 n) {
	return dot(p, n);
}

//rounded box from iq
float rbox(vec3 p, vec3 s, float r)
{
    return length(max(abs(p)-s+vec3(r),0.0))-r;
}

//Ugly hexagonal prism implementation, r radius, h = height * 2
float hex(vec3 p, float r,float h) {
	vec3 n = vec3(1.0,0.0,0.0);
	float a = 3.1415/3.0;
	float c = cos(a);
	float s = sin(a);
	float d = plane(p - n * r, n);
	rYCS(n, c, s);
	d = max(d, plane(p - n * r, n));
	rYCS(n, c, s);
	d = max(d, plane(p - n * r, n));
	rYCS(n, c, s);
	d = max(d, plane(p - n * r, n));
	rYCS(n, c, s);
	d = max(d, plane(p - n * r, n));
	rYCS(n, c, s);
	d = max(d, plane(p - n * r, n));
	d = max(d,-p.y-h);
	d = max(d,p.y-h);
	return d;
}

//raymarching, p - position/origin, d - direction, s - stepping multiplier (0.25-1.0), t - threshold, n maximum number of steps
//returns a vector vec3(steps/(n-1), stepped length, last stepped length) 
vec3 rm(inout vec3 p, inout vec3 d, float s, float t, float n) {
	vec3 q=p;
	float l,i;
	for (i=0.;i<1.;i+=1.0/n) {
		l=abs(f(p));
		p+=l*d*s;
		if (l<t) break;
	};
	return vec3(i,length(q-p),l);
}

added on the 2011-02-08 23:29:24 by las

Thanks las! Thats exactly where I wanted to go with this thread. :)

added on the 2011-02-08 23:52:20 by Mewler

Code:


//r1-lower radius, r2 upper radius, c = height/2
float cone(vec3 p, float r1, float r2,float c) {
	float d = length(p.xz)-mix(r1, r2, (c+p.y)/(c+c));
	d = max(d,-p.y-c);
	d = max(d,p.y-c);
	return d;
}

//signed distance to box (iq) - glsl interpretation
float sbox(vec3 p, vec3 s) {
	vec3 d = abs(p) - s;
	return min(max(d.x,max(d.y, d.z)), length(max(d, 0.0)));
}

added on the 2011-02-09 00:33:35 by las

Quote:

You can do fake SSS with that AO: instead of ambientOcc = ao(p,n) try sss = ao(p,d) where d is the ray direction

Correction: sss = 1.0 - ao(p,d); - you can also refract d first and use a far too big stepsize and you will end up with some really weird but maybe cool results :D

added on the 2011-02-09 00:57:05 by las

Link collection to other (maybe) useful RM threads and some papers.

Pouet threads
Raymarching Beginners' Thread
This might be a good place for beginners to solve problems.

So, what do distance field equations look like? And how do we solve them?
Covers texel's optimisation stuff and lots of other techniques - might be one of the most useful threads on RM.

Cdak by Quite - How?
Making of Cdak.

Problem with RayMarching using DistanceFunctions when applying Twist-Distortion.
Stepsize problems.

Raymarching idea: Invertible Surfaces
Yes it's called signed distance functions - and you can do refractions and stuff!

Useful PDFs:
Rendering Worlds With Two Triangles
iq's great RWWTT talk

Sphere tracing: a geometric method for the antialiased ray tracing of implicit surfaces
The Zeno-paper from John C. Hart - 1996

added on the 2011-02-09 04:01:38 by las

Actually damn useful thread. Thanks mewler + las :)

added on the 2011-02-09 12:12:44 by psonice

Nice thread, very useful indeed.

I've been thinking. Maybe we should move all the raymarching threads to another forum. A place where we can make as many threads as we like without anyone complaining and discuss each and any technical aspect in as much detail as we feel inclined to.

Also we'd then be totally free from pouet-style I-crap-on-your-thread intermissions.

How's that sound, raymarchers?

added on the 2011-02-10 00:01:34 by vibrator

Dull? Technical discussion is all well and good, but there's nothing like a cute kitten photo to freshen up your brain when it's getting overworked.

added on the 2011-02-10 00:08:45 by psonice

Please use this only as a toolbox. kthxbye.

Code:


//simple glsl conversion of the normal generation code for c&p...
vec3 normal(vec3 p) {
	vec3 e = vec3(e, 0.0, 0.0);
	vec3 n;
	n.x = f(p + e.xyy) - f(p - e.xyy);
	n.y = f(p + e.yxy) - f(p - e.yxy);
	n.z = f(p + e.yyx) - f(p - e.yyx);
	return normalize(n);
}

//simple cylinder
float cyl(vec3 p, float r,float c) {
	float d = length(p.xz)-r;
	d = max(d,-p.y-c);
	d = max(d,p.y-c);
	return d;
}

added on the 2011-02-10 00:20:01 by las

The magnitude of the gradient of a distance field is 1 by definition. There's really no need to normalise:

Code:

float e = …;
vec3 n = vec3
(
  ƒ( p.x+e, p.y, p.z ) - ƒ( p.x-e, p.y, p.z ),
  ƒ( p.x, p.y+e, p.z ) - ƒ( p.x, p.y-e, p.z ),
  ƒ( p.x, p.y, p.z+e ) - ƒ( p.x, p.y, p.z-e )
) / ( 2 * e );

Or maybe even:

Code:

float e = …;
vec3 n = ( vec3( ƒ( p.x+e, p.y, p.z ), ƒ( p.x, p.y+e, p.z ), ƒ( p.x, p.y, p.z+e ) ) - p ) / e;

added on the 2011-02-10 13:16:26 by doomdoom

Eh... not - p but - vec3( ƒ(p), ƒ(p), ƒ(p) ), of course. Head hurt. It's a 4-point sample.

added on the 2011-02-10 13:20:45 by doomdoom

GLSL Verision

Code:


vec3 g(vec3 p) {
	vec2 e = vec2(0.01, 0.0);
	return (vec3(f(p+e.xyy),f(p+e.yxy),f(p+e.yyx))-f(p))/e.x;
}

added on the 2011-02-10 19:07:49 by las

Generalized Distance Functions

Nice paper on shape generation.

Code:


vec3 n1 = vec3(1.000,0.000,0.000);
vec3 n2 = vec3(0.000,1.000,0.000);
vec3 n3 = vec3(0.000,0.000,1.000);
vec3 n4 = vec3(0.577,0.577,0.577);
vec3 n5 = vec3(-0.577,0.577,0.577);
vec3 n6 = vec3(0.577,-0.577,0.577);
vec3 n7 = vec3(0.577,0.577,-0.577);
vec3 n8 = vec3(0.000,0.357,0.934);
vec3 n9 = vec3(0.000,-0.357,0.934);
vec3 n10 = vec3(0.934,0.000,0.357);
vec3 n11 = vec3(-0.934,0.000,0.357);
vec3 n12 = vec3(0.357,0.934,0.000);
vec3 n13 = vec3(-0.357,0.934,0.000);
vec3 n14 = vec3(0.000,0.851,0.526);
vec3 n15 = vec3(0.000,-0.851,0.526);
vec3 n16 = vec3(0.526,0.000,0.851);
vec3 n17 = vec3(-0.526,0.000,0.851);
vec3 n18 = vec3(0.851,0.526,0.000);
vec3 n19 = vec3(-0.851,0.526,0.000);

// p as usual, e exponent (p in the paper), r radius or something like that
float octahedral(vec3 p, float e, float r) {
	float s = pow(abs(dot(p,n4)),e);
	s += pow(abs(dot(p,n5)),e);
	s += pow(abs(dot(p,n6)),e);
	s += pow(abs(dot(p,n7)),e);
	s = pow(s, 1./e);
	return s-r;
}


float dodecahedral(vec3 p, float e, float r) {
	float s = pow(abs(dot(p,n14)),e);
	s += pow(abs(dot(p,n15)),e);
	s += pow(abs(dot(p,n16)),e);
	s += pow(abs(dot(p,n17)),e);
	s += pow(abs(dot(p,n18)),e);
	s += pow(abs(dot(p,n19)),e);
	s = pow(s, 1./e);
	return s-r;
}

float icosahedral(vec3 p, float e, float r) {
	float s = pow(abs(dot(p,n4)),e);
	s += pow(abs(dot(p,n5)),e);
	s += pow(abs(dot(p,n6)),e);
	s += pow(abs(dot(p,n7)),e);
	s += pow(abs(dot(p,n8)),e);
	s += pow(abs(dot(p,n9)),e);
	s += pow(abs(dot(p,n10)),e);
	s += pow(abs(dot(p,n11)),e);
	s += pow(abs(dot(p,n12)),e);
	s += pow(abs(dot(p,n13)),e);
	s = pow(s, 1./e);
	return s-r;
}

float toctahedral(vec3 p, float e, float r) {
	float s = pow(abs(dot(p,n1)),e);
	s += pow(abs(dot(p,n2)),e);
	s += pow(abs(dot(p,n3)),e);
	s += pow(abs(dot(p,n4)),e);
	s += pow(abs(dot(p,n5)),e);
	s += pow(abs(dot(p,n6)),e);
	s += pow(abs(dot(p,n7)),e);
	s = pow(s, 1./e);
	return s-r;
}

float ticosahedral(vec3 p, float e, float r) {
	float s = pow(abs(dot(p,n4)),e);
	s += pow(abs(dot(p,n5)),e);
	s += pow(abs(dot(p,n6)),e);
	s += pow(abs(dot(p,n7)),e);
	s += pow(abs(dot(p,n8)),e);
	s += pow(abs(dot(p,n9)),e);
	s += pow(abs(dot(p,n10)),e);
	s += pow(abs(dot(p,n11)),e);
	s += pow(abs(dot(p,n12)),e);
	s += pow(abs(dot(p,n13)),e);
	s += pow(abs(dot(p,n14)),e);
	s += pow(abs(dot(p,n15)),e);
	s += pow(abs(dot(p,n16)),e);
	s += pow(abs(dot(p,n17)),e);
	s += pow(abs(dot(p,n18)),e);
	s += pow(abs(dot(p,n19)),e);
	s = pow(s, 1./e);
	return s-r;
}

For higher values of e (>~50) it seems to be kinda "unstable". I guess it's the float precision...

added on the 2011-02-10 20:28:58 by las

Who is moderating pouet.net nowadays? Please state names and church.

added on the 2011-02-10 20:58:49 by Hatikvah

truncated octahedral

added on the 2011-02-10 21:21:09 by las

Quote:

How's that sound, raymarchers?

Did that just imply that people like Kusma and I are *not* raymarchers? Know your facts, sir!!

...in all seriousness, though, dbfinteractive.com is a great forum for this kind of thing. It'd be much appreciated there; very helpful guys.

added on the 2011-02-10 23:28:12 by ferris

las: if you want "high values of e", just take the limit directly:

Code:


  s = abs(dot(p, n4));
  s = max(s, abs(dot(p, n5)));
  s = max(s, abs(dot(p, n6)));
  // ...
  return s-r;

added on the 2011-02-11 03:30:59 by ryg

also, this:

Code:


//simple cylinder
float cyl(vec3 p, float r,float c) {
	float d = length(p.xz)-r;
	d = max(d,-p.y-c);
	d = max(d,p.y-c);
	return d;
}

is just... eww! c'mon! basic algebra ftw!

Code:


float cyl(vec3 p, float r, float c) {
    return max(length(p.xz)-r, abs(p.y)-c);
}

added on the 2011-02-11 03:36:28 by ryg

That cylinder thing is of course right... *cough*
+you are of course also right about that limit stuff.
I must have been under the influence of shiny glowing things.
But sometimes one just wants to have that exponent :)
Nice, thank you ryg.

Code:


//Hard edges on the truncated octahedral - ryg style
float toctahedral(vec3 p, float r) {
	float s = abs(dot(p,n1));
	s = max(s,abs(dot(p,n2)));
	s = max(s,abs(dot(p,n3)));
	s = max(s,abs(dot(p,n4)));
	s = max(s,abs(dot(p,n5)));
	s = max(s,abs(dot(p,n6)));
	s = max(s,abs(dot(p,n7)));
	return s-r;
}

added on the 2011-02-11 08:39:40 by las

Cheap GLSL noise:

Code:


float rand(vec2 coordinate) {
    return fract(sin(dot(coordinate.xy, vec2(12.9898, 78.233))) * 43758.5453);
}

added on the 2011-02-11 10:04:53 by raer

very cool thread... thanks las and others

added on the 2011-02-11 10:07:54 by auld

Fancier: 2D/3D/4D simplex noise implementation based on the new algorithm by Ken Perlin. Too long to post.

added on the 2011-02-11 10:23:36 by raer

pouët.net

Raymarching Toolbox Thread

login