website articles
sphere functions

Intro



If you are an engineer, you might have heard this before:


"Well, first let us consider a spherical cow..."



This is a humorous observation of the way engineers often operate. But the truth is that replacing precise problems that have only approximated, expensive or even unknown solution by approximated problems that do have an exact, cheap or known problems is sometimes a good strategy. And so, in the context of graphics, sometimes you can replace a cow, indeed, by a sphere (for example, if it's not too close to the camera, or if you are in an occlusion pass).

In fact, spheres are useful in too many ways to be ignored. I have written a few articles already on the subject of computing analytically different properties of spheres that are useful for rendering, which are scattered over this website. Some others I have not documented yet, but this page can be seen as an index to those articles and techniques, and a summary of the resulting formulas. I do not recommend to simply copy&paste them without reading the associated article first so you understand the context where they can be used, and their strengths and pitfalls.

I have used these techniques both in demos and also in production of featured films. For example, I have used a bunch of spheres to shape its occlusion and directional shadows analytically rather than resorting to shadow maps. Also, if you are rendering massive terrains made of hundreds of millions of rocks, bushes and trees, you might want to organize things in bounding spheres and check analytically (without rasterization whatsoever), how many pixels does each bounding sphere take in screen space in to do level of detail decisions. I have used spherical density to do localized color grading, and for fog volumes.

Whether you are using spheres to approximate cows, or you are making a demo or game with spheres for artistic reasons, probably some of the formulas bellow will be of use to you.



Functions



In the code below, pi means 3.141593... and eps is a small positive number close to zero.



Sphere Intersection
float sphIntersect( vec3 ro, vec3 rd, vec4 sph )
{
    vec3 oc = ro - sph.xyz;
    float b = dot( oc, rd );
    float c = dot( oc, oc ) - sph.w*sph.w;
    float h = b*b - c;
    if( h<0.0 ) return -1.0;
    h = sqrt( h );
    return -b - h;
}
http://www.iquilezles.org/blog/?p=2411

Sphere Density
float sphDensity( vec3 ro, vec3 rd, vec4 sph, float dbuffer )
{
    float ndbuffer = dbuffer/sph.w;
    vec3  rc = (ro - sph.xyz)/sph.w;
	
    float b = dot(rd,rc);
    float c = dot(rc,rc) - 1.0;
    float h = b*b - c;
    if( h<0.0 ) return 0.0;
    h = sqrt( h );
    float t1 = -b - h;
    float t2 = -b + h;

    if( t2<0.0 || t1>ndbuffer ) return 0.0;
    t1 = max( t1, 0.0 );
    t2 = min( t2, ndbuffer );

    float i1 = -(c*t1 + b*t1*t1 + t1*t1*t1/3.0);
    float i2 = -(c*t2 + b*t2*t2 + t2*t2*t2/3.0);
    return (i2-i1)*(3.0/4.0);
}
Maths / Article

Example / Code

Sphere Projection
float projectSphere( vec4 sph, in mat4 cam, float fle )
{
    vec3  o = (cam*vec4(sph.xyz,1.0)).xyz;
	
    float r2 = sph.w*sph.w;
    float z2 = o.z*o.z;	
    float l2 = dot(o,o);
    float k1 = l2-r2;
    float k2 = r2-z2;
	
    return -pi*fle*fle*r2*sqrt(abs(k1/k2))/k2;
}
Maths / Article

Example / Code

Sphere Visibility
int sphereVisibility( vec4 sA, vec4 sB, vec3 c )
{
    vec3 ac = sA.xyz - c;
    vec3 bc = sB.xyz - c;
    float ia = 1.0/length(ac);
    float ib = 1.0/length(bc);
    float k0 = dot(ac,bc)*ia*ib;
    float k1 = sA.w*ia;
    float k2 = sB.w*ib;
    float m1 = k0*k0 + k1*k1 + k2*k2;
    float m2 = 2.0*k0*k1*k2;

         if( m1 + m2 - 1.0 < 0.0 ) return 1;
    else if( m1 - m2 - 1.0 < 0.0 ) return 2;
    return 3;
}
Maths / Article

Example / Code

Sphere Ambient Occlusion
float sphOcclusion( vec3 pos, vec3 nor, vec4 sph )
{
    vec3  di = sph.xyz - pos;
    float l  = length(di);
    float nl = dot(nor,di/l);
    float h  = l/sph.w;
    float h2 = h*h;
    float k2 = 1.0 - h2*nl*nl;

    float res = max(0.0,nl)/h2;
    if( k2 > 0.0 ) // approx. for penetration
    {
        res = clamp(0.5*(nl*h+1.0)/h2,0.0,1.0);
        res = sqrt( res*res*res );
    }
    return res;
}
Maths / Article

Example / Code

Sphere Soft Shadow
float sphSoftShadow( vec3 ro, vec3 rd, vec4 sph, float k )
{
    vec3 oc = ro - sph.xyz;
    float b = dot( oc, rd );
    float c = dot( oc, oc ) - sph.w*sph.w;
    float h = b*b - c;
    return (b>0.0)?step(-eps,c):smoothstep(0.0,1.0,h*k/b);
}

float sphSoftShadow( vec3 ro, vec3 rd, vec4 sph, float k )
{
    vec3 oc = ro - sph.xyz;
    float r = sph.w*sph.w;
    float b = dot( oc, rd );
    float c = dot( oc, oc ) - r;
    float h = b*b - c;
    float d = -sph.w + sqrt( max(0.0,r-h));
    float t = -b     - sqrt( max(0.0,h) );
    return (t<0.0)?1.0:smoothstep(0.0,1.0,k*d/t);
}
Maths / Article

Example / Code

Sphere Distance
vec2 sphDistances( vec3 ro, vec3 rd, vec4 sph )
{
    vec3 oc = ro - sph.xyz;
    float b = dot( oc, rd );
    float c = dot( oc, oc ) - sph.w*sph.w;
    float h = b*b - c;
    float d = sqrt( max(0.0,sph.w*sph.w-h)) - sph.w;
    return vec2( d, -b-sqrt(max(h,0.0)) );
}
Example / Code

Sphere Motion
vec2 sphMotion( vec3 ro, vec3 rd, vec4 sph, vec3 ve, out vec3 nor )
{
    float t = -1.0;
    float s = 0.0;
    nor = vec3(0.0);
	
    vec3  rc = ro - sph.xyz;
    float A = dot(rc,rd);
    float B = dot(rc,rc) - sph.w*sph.w;
    float C = dot(ve,ve);
    float D = dot(rc,ve);
    float E = dot(rd,ve);
    float aab = A*A - B;
    float eec = E*E - C;
    float aed = A*E - D;
    float k = aed*aed - eec*aab;
		
    if( k>0.0 )
    {
        k = sqrt(k);
        float hb = (aed - k)/eec;
        float ha = (aed + k)/eec;
        float ta = max( 0.0, ha );
        float tb = min( 1.0, hb );
		
        if( ta < tb )
        {
            ta = 0.5*(ta+tb);
            float k1 = A - E*ta;
            float k2 = B + C*ta*ta - 2.0*D*ta;
            t = -k1 - sqrt( k1*k1 - k2 );
            s = 2.0*(tb - ta);
			
            vec3 pa = ro + rd*t;
            vec3 pb = sph.xyz + ta*ve;
            nor = normalize( pa - pb );
        }
    }

    return vec2(t,s);
}
Example / Code