Intro
This is a collection of ray-surface intersectors for some common primitive types. You can use them to implement a pathtracer, UI clicking, sound propagation, collision detection or anything that requires analytic raycasting. I have derived these myself, and while I have used them professionally, I haven't put them to the most strict of the tests. If you notice precision errors in any of them you can choose to improve these implementations (at some performance cost) or improve how your renderer generally handles precision problems. When possible, I've tried to get the number of operations to a minimum (for example the capped cone and capsule intersectors only check against one of the caps, even though these primitives have two).
Note that all the primitives listed here come with a link to a realtime online demo in Shadertoy. In fact, this public playlist contains all the shadertoy examples: https://www.shadertoy.com/playlist/l3dXRf, so don't miss that out.
Primitives
In all the examples below, ro is the ray origin, and rd is the ray direction. The functions always return the distance to the closest intersection, and sometimes some additional information (second intersection, normal or barycentric coordinates). Most primitives are centered at the origin (you'll need to transform the ray origin and direction accordingly), but some of the functions above accept an arbitrary primitive location and orientation when such implementation is faster than performing a ray transformation. Most of these functions only work for a ray with origin outside of the primitive. Fixing that shouldn't be difficult.
// sphere of size ra centered at point ce
vec2 sphIntersect( in vec3 ro, in vec3 rd, in vec3 ce, float ra )
{
vec3 oc = ro - ce;
float b = dot( oc, rd );
float c = dot( oc, oc ) - ra*ra;
float h = b*b - c;
if( h<0.0 ) return vec2(-1.0); // no intersection
h = sqrt( h );
return vec2( -b-h, -b+h );
}
Alternative method computes h (the squared distance from the closest ray point to the sphere, qc below) with a projection rather than by using Pythagoras' theorem. This is less precision hungry because we don't generate large numbers (in comparison to the size of the sphere) since we don't square triangle edges (b*b):
// sphere of size ra centered at point ce
vec2 sphIntersect( in vec3 ro, in vec3 rd, in vec3 ce, float ra )
{
vec3 oc = ro - ce;
float b = dot( oc, rd );
vec3 qc = oc - b*rd;
float h = ra*ra - dot( qc, qc );
if( h<0.0 ) return vec2(-1.0); // no intersection
h = sqrt( h );
return vec2( -b-h, -b+h );
}
// axis aligned box centered at the origin, with size boxSize
vec2 boxIntersection( in vec3 ro, in vec3 rd, vec3 boxSize, out vec3 outNormal )
{
vec3 m = 1.0/rd; // can precompute if traversing a set of aligned boxes
vec3 n = m*ro; // can precompute if traversing a set of aligned boxes
vec3 k = abs(m)*boxSize;
vec3 t1 = -n - k;
vec3 t2 = -n + k;
float tN = max( max( t1.x, t1.y ), t1.z );
float tF = min( min( t2.x, t2.y ), t2.z );
if( tN>tF || tF<0.0) return vec2(-1.0); // no intersection
outNormal = (tN>0.0) ? step(vec3(tN),t1)) : // ro ouside the box
step(t2,vec3(tF))); // ro inside the box
outNormal *= -sign(rd);
return vec2( tN, tF );
}
// axis aligned box centered at the origin, with dimensions "size" and extruded by "rad"
float roundedboxIntersect( in vec3 ro, in vec3 rd, in vec3 size, in float rad )
{
// bounding box
vec3 m = 1.0/rd;
vec3 n = m*ro;
vec3 k = abs(m)*(size+rad);
vec3 t1 = -n - k;
vec3 t2 = -n + k;
float tN = max( max( t1.x, t1.y ), t1.z );
float tF = min( min( t2.x, t2.y ), t2.z );
if( tN>tF || tF<0.0) return -1.0;
float t = tN;
// convert to first octant
vec3 pos = ro+t*rd;
vec3 s = sign(pos);
ro *= s;
rd *= s;
pos *= s;
// faces
pos -= size;
pos = max( pos.xyz, pos.yzx );
if( min(min(pos.x,pos.y),pos.z) < 0.0 ) return t;
// some precomputation
vec3 oc = ro - size;
vec3 dd = rd*rd;
vec3 oo = oc*oc;
vec3 od = oc*rd;
float ra2 = rad*rad;
t = 1e20;
// corner
{
float b = od.x + od.y + od.z;
float c = oo.x + oo.y + oo.z - ra2;
float h = b*b - c;
if( h>0.0 ) t = -b-sqrt(h);
}
// edge X
{
float a = dd.y + dd.z;
float b = od.y + od.z;
float c = oo.y + oo.z - ra2;
float h = b*b - a*c;
if( h>0.0 )
{
h = (-b-sqrt(h))/a;
if( h>0.0 && h<t && abs(ro.x+rd.x*h)<size.x ) t = h;
}
}
// edge Y
{
float a = dd.z + dd.x;
float b = od.z + od.x;
float c = oo.z + oo.x - ra2;
float h = b*b - a*c;
if( h>0.0 )
{
h = (-b-sqrt(h))/a;
if( h>0.0 && h<t && abs(ro.y+rd.y*h)<size.y ) t = h;
}
}
// edge Z
{
float a = dd.x + dd.y;
float b = od.x + od.y;
float c = oo.x + oo.y - ra2;
float h = b*b - a*c;
if( h>0.0 )
{
h = (-b-sqrt(h))/a;
if( h>0.0 && h<t && abs(ro.z+rd.z*h)<size.z ) t = h;
}
}
if( t>1e19 ) t=-1.0;
return t;
}
// normal of a rounded box
vec3 roundedboxNormal( in vec3 pos, in vec3 siz, in float rad )
{
return sign(pos)*normalize(max(abs(pos)-siz,0.0));
}
Plane
// plane degined by p (p.xyz must be normalized)
float plaIntersect( in vec3 ro, in vec3 rd, in vec4 p )
{
return -(dot(ro,p.xyz)+p.w)/dot(rd,p.xyz);
}
// disk center c, normal n, radius r
float diskIntersect( in vec3 ro, in vec3 rd, vec3 c, vec3 n, float r )
{
vec3 o = ro - c;
float t = -dot(n,o)/dot(rd,n);
vec3 q = o + rd*t;
return (dot(q,q)<r*r) ? t : -1.0;
}
// returns t and normal
vec4 iHexPrism( in vec3 ro, in vec3 rd, in float ra, in float he )
{
const float ks3 = 0.866025;
// normals
const vec3 n1 = vec3( 1.0,0.0,0.0);
const vec3 n2 = vec3( 0.5,0.0,ks3);
const vec3 n3 = vec3(-0.5,0.0,ks3);
const vec3 n4 = vec3( 0.0,1.0,0.0);
// slabs intersections
vec3 t1 = vec3((vec2(ra,-ra)-dot(ro,n1))/dot(rd,n1), 1.0);
vec3 t2 = vec3((vec2(ra,-ra)-dot(ro,n2))/dot(rd,n2), 1.0);
vec3 t3 = vec3((vec2(ra,-ra)-dot(ro,n3))/dot(rd,n3), 1.0);
vec3 t4 = vec3((vec2(he,-he)-dot(ro,n4))/dot(rd,n4), 1.0);
// inetsection selection
if( t1.y<t1.x ) t1=vec3(t1.yx,-1.0);
if( t2.y<t2.x ) t2=vec3(t2.yx,-1.0);
if( t3.y<t3.x ) t3=vec3(t3.yx,-1.0);
if( t4.y<t4.x ) t4=vec3(t4.yx,-1.0);
vec4 tN=vec4(t1.x,t1.z*n1);
if( t2.x>tN.x ) tN=vec4(t2.x,t2.z*n2);
if( t3.x>tN.x ) tN=vec4(t3.x,t3.z*n3);
if( t4.x>tN.x ) tN=vec4(t4.x,t4.z*n4);
float tF = min(min(t1.y,t2.y),min(t3.y,t4.y));
// no intersection
if( tN.x>tF || tF<0.0) return vec4(-1.0);
return tN; // return tF too for exit point
}
// returns t and normal
vec4 iWedge( in vec3 ro, in vec3 rd, in vec3 s )
{
// intersect plane box
vec3 m = 1.0/rd;
vec3 z = vec3(rd.x>=0.0?1.0:-1.0, rd.y>=0.0?1.0:-1.0, rd.z>=0.0?1.0:-1.0);
vec3 k = s*z;
vec3 t1 = (-ro - k)*m;
vec3 t2 = (-ro + k)*m;
float tn = max(max(t1.x, t1.y), t1.z);
float tf = min(min(t2.x, t2.y), t2.z);
if( tn>tf ) return vec4(-1.0);
// boolean with plane
float k1 = s.y*ro.x - s.x*ro.y;
float k2 = s.x*rd.y - s.y*rd.x;
float tp = k1/k2;
if( k1>tn*k2 )
return vec4(tn,-step(tn,t1)*z); // box
if( tp>tn && tp<tf )
return vec4(tp,normalize(vec3(-s.y,s.x,0.0))); // plane
return vec4(-1.0);
}
// capsule defined by extremes pa and pb, and radious ra
// Note that only ONE of the two spherical caps is checked for intersections,
// which is a nice optimization
float capIntersect( in vec3 ro, in vec3 rd, in vec3 pa, in vec3 pb, in float ra )
{
vec3 ba = pb - pa;
vec3 oa = ro - pa;
float baba = dot(ba,ba);
float bard = dot(ba,rd);
float baoa = dot(ba,oa);
float rdoa = dot(rd,oa);
float oaoa = dot(oa,oa);
float a = baba - bard*bard;
float b = baba*rdoa - baoa*bard;
float c = baba*oaoa - baoa*baoa - ra*ra*baba;
float h = b*b - a*c;
if( h >= 0.0 )
{
float t = (-b-sqrt(h))/a;
float y = baoa + t*bard;
// body
if( y>0.0 && y<baba ) return t;
// caps
vec3 oc = (y <= 0.0) ? oa : ro - pb;
b = dot(rd,oc);
c = dot(oc,oc) - ra*ra;
h = b*b - c;
if( h>0.0 ) return -b - sqrt(h);
}
return -1.0;
}
// cylinder defined by extremes a and b, and radious ra
vec4 cylIntersect( in vec3 ro, in vec3 rd, in vec3 a, in vec3 b, float ra )
{
vec3 ba = b - a;
vec3 oc = ro - a;
float baba = dot(ba,ba);
float bard = dot(ba,rd);
float baoc = dot(ba,oc);
float k2 = baba - bard*bard;
float k1 = baba*dot(oc,rd) - baoc*bard;
float k0 = baba*dot(oc,oc) - baoc*baoc - ra*ra*baba;
float h = k1*k1 - k2*k0;
if( h<0.0 ) return vec4(-1.0);//no intersection
h = sqrt(h);
float t = (-k1-h)/k2;
// body
float y = baoc + t*bard;
if( y>0.0 && y<baba ) return vec4( t, (oc+t*rd - ba*y/baba)/ra );
// caps
t = ( ((y<0.0) ? 0.0 : baba) - baoc)/bard;
if( abs(k1+k2*t)<h )
{
return vec4( t, ba*sign(y)/sqrt(baba) );
}
return vec4(-1.0);//no intersection
}
// normal at point p of cylinder (a,b,ra), see above
vec3 cylNormal( in vec3 p, in vec3 a, in vec3 b, float ra )
{
vec3 pa = p - a;
vec3 ba = b - a;
float baba = dot(ba,ba);
float paba = dot(pa,ba);
float h = dot(pa,ba)/baba;
return (pa - ba*h)/ra;
}
Cylinder
// infinite cylinder defined by a base point cb, a normalized axis ca and a radious cr
vec2 cylIntersect( in vec3 ro, in vec3 rd, in vec3 cb, in vec3 ca, float cr )
{
vec3 oc = ro - cb;
float card = dot(ca,rd);
float caoc = dot(ca,oc);
float a = 1.0 - card*card;
float b = dot( oc, rd) - caoc*card;
float c = dot( oc, oc) - caoc*caoc - cr*cr;
float h = b*b - a*c;
if( h<0.0 ) return vec2(-1.0); //no intersection
h = sqrt(h);
return vec2(-b-h,-b+h)/a;
}
// cone defined by extremes pa and pb, and radious ra and rb
Only one square root and one division is emplyed in the worst case. dot2(v) is dot(v,v)
vec4 coneIntersect( in vec3 ro, in vec3 rd, in vec3 pa, in vec3 pb, in float ra, in float rb )
{
vec3 ba = pb - pa;
vec3 oa = ro - pa;
vec3 ob = ro - pb;
float m0 = dot(ba,ba);
float m1 = dot(oa,ba);
float m2 = dot(rd,ba);
float m3 = dot(rd,oa);
float m5 = dot(oa,oa);
float m9 = dot(ob,ba);
// caps
if( m1<0.0 )
{
if( dot2(oa*m2-rd*m1)<(ra*ra*m2*m2) ) // delayed division
return vec4(-m1/m2,-ba*inversesqrt(m0));
}
else if( m9>0.0 )
{
float t = -m9/m2; // NOT delayed division
if( dot2(ob+rd*t)<(rb*rb) )
return vec4(t,ba*inversesqrt(m0));
}
// body
float rr = ra - rb;
float hy = m0 + rr*rr;
float k2 = m0*m0 - m2*m2*hy;
float k1 = m0*m0*m3 - m1*m2*hy + m0*ra*(rr*m2*1.0 );
float k0 = m0*m0*m5 - m1*m1*hy + m0*ra*(rr*m1*2.0 - m0*ra);
float h = k1*k1 - k2*k0;
if( h<0.0 ) return vec4(-1.0); //no intersection
float t = (-k1-sqrt(h))/k2;
float y = m1 + t*m2;
if( y<0.0 || y>m0 ) return vec4(-1.0); //no intersection
return vec4(t, normalize(m0*(m0*(oa+t*rd)+rr*ba*ra)-ba*hy*y));
}
// cone defined by extremes pa and pb, and radious ra and rb.
vec4 iRoundedCone( in vec3 ro, in vec3 rd, in vec3 pa, in vec3 pb, in float ra, in float rb )
{
vec3 ba = pb - pa;
vec3 oa = ro - pa;
vec3 ob = ro - pb;
float rr = ra - rb;
float m0 = dot(ba,ba);
float m1 = dot(ba,oa);
float m2 = dot(ba,rd);
float m3 = dot(rd,oa);
float m5 = dot(oa,oa);
float m6 = dot(ob,rd);
float m7 = dot(ob,ob);
// body
float d2 = m0-rr*rr;
float k2 = d2 - m2*m2;
float k1 = d2*m3 - m1*m2 + m2*rr*ra;
float k0 = d2*m5 - m1*m1 + m1*rr*ra*2.0 - m0*ra*ra;
float h = k1*k1 - k0*k2;
if( h<0.0) return vec4(-1.0);
float t = (-sqrt(h)-k1)/k2;
//if( t<0.0 ) return vec4(-1.0);
float y = m1 - ra*rr + t*m2;
if( y>0.0 && y<d2 ) return vec4(t, normalize(d2*(oa+t*rd)-ba*y));
// caps
float h1 = m3*m3 - m5 + ra*ra;
float h2 = m6*m6 - m7 + rb*rb;
if( max(h1,h2)<0.0 ) return vec4(-1.0);
vec4 r = vec4(1e20);
if( h1>0.0 )
{
t = -m3 - sqrt( h1 );
r = vec4( t, (oa+t*rd)/ra );
}
if( h2>0.0 )
{
t = -m6 - sqrt( h2 );
if( t<r.x )
r = vec4( t, (ob+t*rd)/rb );
}
return r;
}
// ellipsoid centered at the origin with radii ra
vec2 eliIntersect( in vec3 ro, in vec3 rd, in vec3 ra )
{
vec3 ocn = ro/ra;
vec3 rdn = rd/ra;
float a = dot( rdn, rdn );
float b = dot( ocn, rdn );
float c = dot( ocn, ocn );
float h = b*b - a*(c-1.0);
if( h<0.0 ) return vec2(-1.0); //no intersection
h = sqrt(h);
return vec2(-b-h,-b+h)/a;
}
vec3 eliNormal( in vec3 pos, in vec3 ra )
{
return normalize( pos/(ra*ra) );
}
// triangle degined by vertices v0, v1 and v2
vec3 triIntersect( in vec3 ro, in vec3 rd, in vec3 v0, in vec3 v1, in vec3 v2 )
{
vec3 v1v0 = v1 - v0;
vec3 v2v0 = v2 - v0;
vec3 rov0 = ro - v0;
vec3 n = cross( v1v0, v2v0 );
vec3 q = cross( rov0, rd );
float d = 1.0/dot( rd, n );
float u = d*dot( -q, v2v0 );
float v = d*dot( q, v1v0 );
float t = d*dot( -n, rov0 );
if( u<0.0 || v<0.0 || (u+v)>1.0 ) t = -1.0;
return vec3( t, u, v );
}
// ellipse defined by its center c and its two radii u and v
vec3 ellIntersect( in vec3 ro, in vec3 rd, vec3 c, vec3 u, vec3 v )
{
vec3 q = ro - c;
vec3 n = cross(u,v);
float t = -dot(n,q)/dot(rd,n);
float r = dot(u,q + rd*t);
float s = dot(v,q + rd*t);
if( r*r+s*s>1.0 ) return vec3(-1.0); //no intersection
return vec3(t,s,r);
}
vec3 ellNormal( in vec2 u, vec2 v )
{
return normalize( cross(u,v) );
}
float torIntersect( in vec3 ro, in vec3 rd, in vec2 tor )
{
float po = 1.0;
float Ra2 = tor.x*tor.x;
float ra2 = tor.y*tor.y;
float m = dot(ro,ro);
float n = dot(ro,rd);
float k = (m + Ra2 - ra2)/2.0;
float k3 = n;
float k2 = n*n - Ra2*dot(rd.xy,rd.xy) + k;
float k1 = n*k - Ra2*dot(rd.xy,ro.xy);
float k0 = k*k - Ra2*dot(ro.xy,ro.xy);
if( abs(k3*(k3*k3-k2)+k1) < 0.01 )
{
po = -1.0;
float tmp=k1; k1=k3; k3=tmp;
k0 = 1.0/k0;
k1 = k1*k0;
k2 = k2*k0;
k3 = k3*k0;
}
float c2 = k2*2.0 - 3.0*k3*k3;
float c1 = k3*(k3*k3-k2)+k1;
float c0 = k3*(k3*(c2+2.0*k2)-8.0*k1)+4.0*k0;
c2 /= 3.0;
c1 *= 2.0;
c0 /= 3.0;
float Q = c2*c2 + c0;
float R = c2*c2*c2 - 3.0*c2*c0 + c1*c1;
float h = R*R - Q*Q*Q;
if( h>=0.0 )
{
h = sqrt(h);
float v = sign(R+h)*pow(abs(R+h),1.0/3.0); // cube root
float u = sign(R-h)*pow(abs(R-h),1.0/3.0); // cube root
vec2 s = vec2( (v+u)+4.0*c2, (v-u)*sqrt(3.0));
float y = sqrt(0.5*(length(s)+s.x));
float x = 0.5*s.y/y;
float r = 2.0*c1/(x*x+y*y);
float t1 = x - r - k3; t1 = (po<0.0)?2.0/t1:t1;
float t2 = -x - r - k3; t2 = (po<0.0)?2.0/t2:t2;
float t = 1e20;
if( t1>0.0 ) t=t1;
if( t2>0.0 ) t=min(t,t2);
return t;
}
float sQ = sqrt(Q);
float w = sQ*cos( acos(-R/(sQ*Q)) / 3.0 );
float d2 = -(w+c2); if( d2<0.0 ) return -1.0;
float d1 = sqrt(d2);
float h1 = sqrt(w - 2.0*c2 + c1/d1);
float h2 = sqrt(w - 2.0*c2 - c1/d1);
float t1 = -d1 - h1 - k3; t1 = (po<0.0)?2.0/t1:t1;
float t2 = -d1 + h1 - k3; t2 = (po<0.0)?2.0/t2:t2;
float t3 = d1 - h2 - k3; t3 = (po<0.0)?2.0/t3:t3;
float t4 = d1 + h2 - k3; t4 = (po<0.0)?2.0/t4:t4;
float t = 1e20;
if( t1>0.0 ) t=t1;
if( t2>0.0 ) t=min(t,t2);
if( t3>0.0 ) t=min(t,t3);
if( t4>0.0 ) t=min(t,t4);
return t;
}
vec3 torNormal( in vec3 pos, vec2 tor )
{
return normalize( pos*(dot(pos,pos)-tor.y*tor.y - tor.x*tor.x*vec3(1.0,1.0,-1.0)));
}
float sph4Intersect( in vec3 ro, in vec3 rd, in float ra )
{
float r2 = ra*ra;
vec3 d2 = rd*rd; vec3 d3 = d2*rd;
vec3 o2 = ro*ro; vec3 o3 = o2*ro;
float ka = 1.0/dot(d2,d2);
float k3 = ka* dot(ro,d3);
float k2 = ka* dot(o2,d2);
float k1 = ka* dot(o3,rd);
float k0 = ka*(dot(o2,o2) - r2*r2);
float c2 = k2 - k3*k3;
float c1 = k1 + 2.0*k3*k3*k3 - 3.0*k3*k2;
float c0 = k0 - 3.0*k3*k3*k3*k3 + 6.0*k3*k3*k2 - 4.0*k3*k1;
float p = c2*c2 + c0/3.0;
float q = c2*c2*c2 - c2*c0 + c1*c1;
float h = q*q - p*p*p;
if( h<0.0 ) return -1.0; //no intersection
float sh = sqrt(h);
float s = sign(q+sh)*pow(abs(q+sh),1.0/3.0); // cuberoot
float t = sign(q-sh)*pow(abs(q-sh),1.0/3.0); // cuberoot
vec2 w = vec2( s+t,s-t );
vec2 v = vec2( w.x+c2*4.0, w.y*sqrt(3.0) )*0.5;
float r = length(v);
return -abs(v.y)/sqrt(r+v.x) - c1/r - k3;
}
vec3 sph4Normal( in vec3 pos )
{
return normalize( pos*pos*pos );
}
float gouIntersect( in vec3 ro, in vec3 rd, in float ka, float kb )
{
float po = 1.0;
vec3 rd2 = rd*rd; vec3 rd3 = rd2*rd;
vec3 ro2 = ro*ro; vec3 ro3 = ro2*ro;
float k4 = dot(rd2,rd2);
float k3 = dot(ro ,rd3);
float k2 = dot(ro2,rd2) - kb/6.0;
float k1 = dot(ro3,rd ) - kb*dot(rd,ro)/2.0;
float k0 = dot(ro2,ro2) + ka - kb*dot(ro,ro);
k3 /= k4;
k2 /= k4;
k1 /= k4;
k0 /= k4;
float c2 = k2 - k3*(k3);
float c1 = k1 + k3*(2.0*k3*k3-3.0*k2);
float c0 = k0 + k3*(k3*(c2+k2)*3.0-4.0*k1);
if( abs(c1) < 0.1*abs(c2) )
{
po = -1.0;
float tmp=k1; k1=k3; k3=tmp;
k0 = 1.0/k0;
k1 = k1*k0;
k2 = k2*k0;
k3 = k3*k0;
c2 = k2 - k3*(k3);
c1 = k1 + k3*(2.0*k3*k3-3.0*k2);
c0 = k0 + k3*(k3*(c2+k2)*3.0-4.0*k1);
}
c0 /= 3.0;
float Q = c2*c2 + c0;
float R = c2*c2*c2 - 3.0*c0*c2 + c1*c1;
float h = R*R - Q*Q*Q;
if( h>0.0 ) // 2 intersections
{
h = sqrt(h);
float s = sign(R+h)*pow(abs(R+h),1.0/3.0); // cube root
float u = sign(R-h)*pow(abs(R-h),1.0/3.0); // cube root
float x = s+u+4.0*c2;
float y = s-u;
float ks = x*x + y*y*3.0;
float k = sqrt(ks);
float t = -0.5*po*abs(y)*sqrt(6.0/(k+x)) - 2.0*c1*(k+x)/(ks+x*k) - k3;
return (po<0.0)?1.0/t:t;
}
// 4 intersections
float sQ = sqrt(Q);
float w = sQ*cos(acos(-R/(sQ*Q))/3.0);
float d2 = -w - c2;
if( d2<0.0 ) return -1.0; //no intersection
float d1 = sqrt(d2);
float h1 = sqrt(w - 2.0*c2 + c1/d1);
float h2 = sqrt(w - 2.0*c2 - c1/d1);
float t1 = -d1 - h1 - k3; t1 = (po<0.0)?1.0/t1:t1;
float t2 = -d1 + h1 - k3; t2 = (po<0.0)?1.0/t2:t2;
float t3 = d1 - h2 - k3; t3 = (po<0.0)?1.0/t3:t3;
float t4 = d1 + h2 - k3; t4 = (po<0.0)?1.0/t4:t4;
float t = 1e20;
if( t1>0.0 ) t=t1;
if( t2>0.0 ) t=min(t,t2);
if( t3>0.0 ) t=min(t,t3);
if( t4>0.0 ) t=min(t,t4);
return t;
}
vec3 gouNormal( in vec3 pos, float ka, float kb )
{
return normalize( 4.0*pos*pos*pos - 2.0*pos*kb );
}