Over the past two weeks I have been diving into the wonderful and mostly just confusing math behind rendering the weirdest objects in our universe, Black Holes. Today I will take you behind the curtain of this crazy world of shaders that I found myself in and my hope is that you'll learn some things too, just as I did.
First I would like to say that I know that It has been a while since I have posted last, and I'm sorry about that but it has not been a quick process to switch from 3 years of coding in scratch to line code, but I did it and am pretty good at it as you can probably see. Also for this entire thing I will be using the free browser website shaderToy to run this code. But enough with all that lets get into what we've all been waiting for, math and coding.
Rendering a simple Scene
Alright, the first thing that we need to understand to do this kind of thing is Ray Marching. This is the process of firing out a ray from every pixel on the screen and calculating what color it is. The way that I normally accomplish this is using SDFs or signed distance functions to tell you how far a point is from the surface of a shape. For example the SDF for a circle is
float sdSphere( vec2 p, float r )
{
return length( p ) - r;
}
where p is the points position and r is the radius, also for this article I will be expecting that you have a somewhat idea of what the shaderToy functions do. You can of course expand this idea to 3D by making p a vec3.
With this and a couple other SDFs which people like Inigo Quilez have worked out a bunch of them we can make a scene that looks like this :
with the code looking like this :
mat2 rot2D( float angle )
{
float c = cos( angle );
float s = sin( angle );
return mat2( -c, s, s, c );
}
float sdBox( vec3 p, vec3 b )
{
vec3 q = abs( p ) - b;
return length( max( q,0.0 ) ) + min( max( q.x,max( q.y,q.z ) ),0.0 );
}
float sdSphere( vec3 p, float r )
{
return length( p ) - r;
}
float map( vec3 p )
{
vec3 q = p;
q -= vec3( 2.0, 0.0, 0.0 );
float sphere = sdSphere( q, 1.0 );
q = p;
q -= vec3( 1.0, 0.0, 2.0 );
float box = sdBox( q, vec3( 0.75 ) );
float sphere2 = sdSphere( p, 2.0 );
return min( sphere, min( box, sphere2 ) );
}
void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
vec2 uv = ( fragCoord * 2.0 - iResolution.xy ) / iResolution.y;
vec2 m = ( iMouse.xy * 2.0 - iResolution.xy ) / iResolution.y * sqrt( 2.0 ) * 2.0;
vec3 col = vec3( 0.0 );
vec3 ro = vec3( 0.0, 0.0, -5.0 );
vec3 rd = normalize( vec3( uv, 1.0 ) );
ro.yz *= rot2D( m.y );
rd.yz *= rot2D( m.y );
ro.xz *= rot2D( m.x );
rd.xz *= rot2D( m.x );
float t = 0.0;
for( int i = 0; i < 80; i ++ )
{
vec3 p = ro + rd * t;
float d = map( p );
t += d;
if( d < 0.001 || t > 100.0 )break;
}
col = vec3( t * 0.2 );
fragColor = vec4( col, 0.0 );
}
Now before you click off, this is really much simpler than it looks, q is just the displaced position of the individual objects and rot2D is just the 2D rotation matrix in code, other than that its just all what we talked about. Now we have our simple scene made using a sphere traced ray marcher which is what this approach is called.
Warping Space
Alright now that we have our simple engine that renders very simple scenes, we are going to switch to an entirely new approach, fixed step ray marching, which is where instead of moving the farthest distance the point is to the surface of an object, we take well, fixed steps, its pretty self explanatory. We did this because it is very hard to simulate gravity when moving different distances each step, you just have to trust me on this.
The revised code for this is :
const float StepSize = 0.05;
mat2 rot2D( float angle )
{
float c = cos( angle );
float s = sin( angle );
return mat2( -c, s, s, c );
}
float sdBox( vec3 p, vec3 b )
{
vec3 q = abs( p ) - b;
return length( max( q,0.0 ) ) + min( max( q.x,max( q.y,q.z ) ),0.0 );
}
float sdSphere( vec3 p, float r )
{
return length( p ) - r;
}
float map( vec3 p )
{
vec3 q = p;
q -= vec3( 2.0, 0.0, 0.0 );
float sphere = sdSphere( q, 1.0 );
q = p;
q -= vec3( 1.0, 0.0, 2.0 );
float box = sdBox( q, vec3( 0.75 ) );
float sphere2 = sdSphere( p, 2.0 );
return min( sphere, min( box, sphere2 ) );
}
void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
vec2 uv = ( fragCoord * 2.0 - iResolution.xy ) / iResolution.y;
vec2 m = ( iMouse.xy * 2.0 - iResolution.xy ) / iResolution.y * sqrt( 2.0 ) * 2.0;
vec3 col = vec3( 1.0 );
vec3 ro = vec3( 0.0, 0.0, -5.0 );
vec3 rd = normalize( vec3( uv, 1.0 ) );
ro.yz *= rot2D( m.y );
rd.yz *= rot2D( m.y );
ro.xz *= rot2D( m.x );
rd.xz *= rot2D( m.x );
vec3 p = ro;
for( int i = 0; i < 200; i ++ )
{
p += rd * StepSize;
float d = map( p );
if( d < 0.001 )
{
col = vec3( 0.0 );
break;
}
}
fragColor = vec4( col, 0.0 );
}
which looks like :
Honestly a bit simpler now. now when we alter the rd or ray direction, the position of the ray won't rotate around the origin or camera position but the last position of the ray. Now before we warp space and time there are a few constants that we need to make :
const G = 1.0;
const c = 10.0;
const vec3 bhPos = vec3( 0.0 );
const float bhMass = 5.0;
const float bhRS = sqrt( 2.0 * G * bhMass / ( c * c ) );
Where G and c are the gravitational constant and the speed of light, bhPos and Mass are pretty self explanatory but bhRS is the Schwarzschild Radius of the black hole using the real physics equation. Now to see the black hole I will remove everything from the map function and make it return 1.0 and I will add a new function called mapBH which looks like this :
float mapBH( vec3 p )
{
vec3 q = p - bhPos;
float bh = sdSphere( q, bhRS );
return bh;
}
I also change ro to be vec3( 0.0, 0.0, -150.0 ); and const float StepSize = 0.2;, if you see nothing, you're on the right track. The reason that there is nothing is that from this distance, such a small radius just won't show up, like for a while even if you turn up the mass to like 500.0. But now we change that, now we actually warp light to our will.
To do this we need to make a new function, curve, this will curve our light around the black hole using Euler's approximation of light. This is the function :
vec3 curve( vec3 dir, vec3 p )
{
vec3 pos = p - bhPos;
float distToBH = length( pos );
vec3 newDir = normalize(
dir - ( pos * StepSize / pow( distToBH, 3.0 ) * bhMass )
);
return newDir;
}
Where pos is the position of the light ray relative to the black holes position and distToBH is what it sounds like, then using this we can approximate the curvature of light around a black hole, but before we run this, there are a few things I would like to add.
- A noise functions for the rays that miss the Black hole witch the two functions for that look like :
float hash3( vec3 p ) {
p = fract( p * 0.1031 );
p += dot( p, p.zyx + 31.32 );
return fract( ( p.x + p.y ) * p.z );
}
float noise3D( vec3 x ) {
vec3 i = floor( x );
vec3 f = fract( x );
vec3 u = f * f * f * ( f * ( f * 6.0 - 15.0 ) + 10.0 );
float a = hash3( i + vec3( 0.0, 0.0, 0.0 ) );
float b = hash3( i + vec3( 1.0, 0.0, 0.0 ) );
float c = hash3( i + vec3( 0.0, 1.0, 0.0 ) );
float d = hash3( i + vec3( 1.0, 1.0, 0.0 ) );
float e = hash3( i + vec3( 0.0, 0.0, 1.0 ) );
float fVal = hash3( i + vec3( 1.0, 0.0, 1.0 ) );
float g = hash3( i + vec3( 0.0, 1.0, 1.0 ) );
float h = hash3( i + vec3( 1.0, 1.0, 1.0 ) );
return mix( mix( mix( a, b, u.x ), mix( c, d, u.x ), u.y ),
mix( mix( e, fVal, u.x ), mix( g, h, u.x ), u.y ), u.z );
}
- a border where things are visible to better improve performance when running the simulation.
- increasing the StepSize to make the black hole more visible.
this all gives the final effect of :
and I also tried some with a simple planet and sun :
and the final code for this section being :
const float StepSize = 0.5;
const float G = 1.0;
const float c = 10.0;
const vec3 bhPos = vec3( 0.0 );
const float bhMass = 5.0;
const float bhRS = sqrt( 2.0 * G * bhMass / ( c * c ) );
mat2 rot2D( float angle )
{
float c = cos( angle );
float s = sin( angle );
return mat2( -c, s, s, c );
}
float hash3( vec3 p ) {
p = fract( p * 0.1031 );
p += dot( p, p.zyx + 31.32 );
return fract( ( p.x + p.y ) * p.z );
}
float noise3D( vec3 x ) {
vec3 i = floor( x );
vec3 f = fract( x );
vec3 u = f * f * f * ( f * ( f * 6.0 - 15.0 ) + 10.0 );
float a = hash3( i + vec3( 0.0, 0.0, 0.0 ) );
float b = hash3( i + vec3( 1.0, 0.0, 0.0 ) );
float c = hash3( i + vec3( 0.0, 1.0, 0.0 ) );
float d = hash3( i + vec3( 1.0, 1.0, 0.0 ) );
float e = hash3( i + vec3( 0.0, 0.0, 1.0 ) );
float fVal = hash3( i + vec3( 1.0, 0.0, 1.0 ) );
float g = hash3( i + vec3( 0.0, 1.0, 1.0 ) );
float h = hash3( i + vec3( 1.0, 1.0, 1.0 ) );
return mix( mix( mix( a, b, u.x ), mix( c, d, u.x ), u.y ),
mix( mix( e, fVal, u.x ), mix( g, h, u.x ), u.y ), u.z );
}
float sdBox( vec3 p, vec3 b )
{
vec3 q = abs( p ) - b;
return length( max( q,0.0 ) ) + min( max( q.x,max( q.y,q.z ) ),0.0 );
}
float sdSphere( vec3 p, float r )
{
return length( p ) - r;
}
float map( vec3 p )
{
return 1.0;
}
float mapBH( vec3 p )
{
vec3 q = p - bhPos;
float bh = sdSphere( q, bhRS );
return bh;
}
vec3 curve( vec3 dir, vec3 p )
{
vec3 pos = p - bhPos;
float distToBH = length( pos );
vec3 newDir = normalize(
dir - ( pos * StepSize / pow( distToBH, 3.0 ) * bhMass )
);
return newDir;
}
vec3 rayTrace( vec2 uv )
{
vec2 m = ( iMouse.xy * 2.0 - iResolution.xy ) / iResolution.y * sqrt( 2.0 ) * 2.0;
vec3 col = vec3( 1.0 );
vec3 ro = vec3( 0.0, 0.0, -150.0 );
vec3 rd = normalize( vec3( uv, 1.0 ) );
ro.yz *= rot2D( 0.0 );
rd.yz *= rot2D( 0.0 );
ro.xz *= rot2D( iTime / 5.0 );
rd.xz *= rot2D( iTime / 5.0 );
vec3 p = ro;
bool collided = false;
for( int i = 0; i < 1000; i ++ )
{
p += rd * StepSize;
float d = map( p );
if( d < 0.001 )
{
col = vec3( 0.2, 0.1, 0.9 );
collided = true;
break;
}
d = mapBH( p );
if( d < 0.001 )
{
col = vec3( 0.0 );
collided = true;
break;
}
rd = curve( rd, p );
}
if( collided == false )
{
col = vec3( pow( noise3D( p / 7.0 ), 10.0 ) * vec3( 0.9, 0.8, 0.6 ) );
}
return col;
}
void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
vec2 uv = ( fragCoord * 2.0 - iResolution.xy ) / iResolution.y;
vec3 col;
if( abs( uv.x ) < 0.9 && abs( uv.y ) < 0.9 )
{
col = rayTrace( uv );
} else {
col = vec3( 0.1 );
}
fragColor = vec4( col, 0.0 );
}
Accretion Disk
A black hole's accretion disk is mostly made out of dust and gas, which just so happens to be the thing that fixed step ray marching is mostly used for, volumetric rendering ( or things like clouds ).
I will admit this was weirdly hard to create in shaderToy but in the end I got to the point where it works for the black hole's Accretion Disk. So what I am going to do to learn this concept, and what I recommend you to do, is go all the way back to the start of this article and create a simple test scene with just a cube and set the color it returns to vec3( 0.39, 1.0, 0.63 );.
Now comes the fun part, I created a float right under ro and rd called density that is set to 0.0;. Now if we collided with anything instead of breaking I will just increase the density by a constant I will call cloudDensity ( keep the col being set to the color you choose ). Then at the end of the loop I will use col *= density; to set the color to be equal to the level of density that that ray went through. But the spice this up a little more I will use the noise3D function and add that to the point multiplied by the cloud density to add some variety and use p * 5.0; for the input position for the noise.
now we have a very simple cube of volumetric noise that increases in density as the ray travels along it. But now is the hard part, at least to me, lighting. For this I just made the accretion disk itself because that is how it works in real life as the gas and dust is moving so fast that It heats up to thousands of degrees kelvin. The short answer to what I did is the illumination is based on the density that the ray goes through.
To go into more detail the equation that I used for this is 'illumination += sqrt( 0.006 / pow( length( q / ( accretionDiskSize ) * 3.0 ), 2.0 ) * 10.0 );' or I = I + sqrt( 0.006 / ( distTo0( q / ADS ) * 3 ) )^2 * 10. Do with this equation what you will.
Then I just used the sdf :
float sdDiskWithHole( vec3 p, float r, float t, float h )
{
float d = length( p.xz ) - r;
float ringDist = abs( d ) - t * 0.5;
vec2 w = vec2( ringDist, abs( p.y ) - h * 0.5 );
return min( max( w.x, w.y ), 0.0 ) + length( max( w, 0.0 ) );
}
And then added back in the black hole and boom, I have the cover of this article. This is the final code with some optimizations that I did not talk about :
const float StepSize = 0.1;
const float Absorption = 0.5;
const float G = 1.0;
const float c = 10.0;
const vec3 bhPos = vec3( 0.0 );
const float bhMass = 2.0;
const float bhRS = sqrt( 2.0 * G * bhMass / ( c * c ) );
const float accretionDiskSize = 15.0;
mat2 rot2D( float angle )
{
float c = cos( angle );
float s = sin( angle );
return mat2( -c, s, s, c );
}
float sdDiskWithHole( vec3 p, float r, float t, float h )
{
float d = length( p.xz ) - r;
float ringDist = abs( d ) - t * 0.5;
vec2 w = vec2( ringDist, abs( p.y ) - h * 0.5 );
return min( max( w.x, w.y ), 0.0 ) + length( max( w, 0.0 ) );
}
float sdSphere( vec3 p, float r )
{
return length( p ) - r;
}
float rand( float seed ) {
return fract( sin( seed * 12.9898 ) * 43758.5453123 );
}
float hash3( vec3 p ) {
p = fract( p * 0.1031 );
p += dot( p, p.zyx + 31.32 );
return fract( ( p.x + p.y ) * p.z );
}
float noise3D( vec3 x ) {
vec3 i = floor( x );
vec3 f = fract( x );
vec3 u = f * f * f * ( f * ( f * 6.0 - 15.0 ) + 10.0 );
float a = hash3( i + vec3( 0.0, 0.0, 0.0 ) );
float b = hash3( i + vec3( 1.0, 0.0, 0.0 ) );
float c = hash3( i + vec3( 0.0, 1.0, 0.0 ) );
float d = hash3( i + vec3( 1.0, 1.0, 0.0 ) );
float e = hash3( i + vec3( 0.0, 0.0, 1.0 ) );
float fVal = hash3( i + vec3( 1.0, 0.0, 1.0 ) );
float g = hash3( i + vec3( 0.0, 1.0, 1.0 ) );
float h = hash3( i + vec3( 1.0, 1.0, 1.0 ) );
return mix( mix( mix( a, b, u.x ), mix( c, d, u.x ), u.y ),
mix( mix( e, fVal, u.x ), mix( g, h, u.x ), u.y ), u.z );
}
float mapGas( vec3 p )
{
vec3 q = p;
q.xy *= rot2D( 0.2 );
float sphere = sdDiskWithHole( q, bhRS * 3.0 + accretionDiskSize / 2.0 + bhMass * 2.0, accretionDiskSize, 1.0 );
return sphere;
}
float mapBH( vec3 p )
{
float bh = sdSphere( p, bhRS );
return bh;
}
vec3 curve( vec3 rd, vec3 p )
{
vec3 pos = p - bhPos;
float distToBH = length( pos );
vec3 newDir = normalize(
rd - ( pos * StepSize / pow( distToBH, 3.0 ) * bhMass )
);
return newDir;
}
vec3 rayTrace( vec2 uv )
{
vec3 col = vec3( 0.0 );
vec3 ro = vec3( 0.0, 0.0, -60.0 );
vec3 rd = normalize( vec3( uv, 5.0 ) );
ro.yz *= rot2D( 0.3 );
rd.yz *= rot2D( 0.3 );
vec3 p = ro;
float d;
float density = 0.0;
float illumination = 0.0;
for( int i = 0; i < int( length( ro ) * 2.0 / StepSize ); i ++ )
{
p += rd * StepSize;
d = mapGas( p );
if( d < 0.0 )
{
col = vec3( 1.0, 0.3, 0.0 );
vec3 q = p;
q.xz *= rot2D( iTime );
vec3 lightDir = normalize( -q );
illumination += sqrt( 0.006 / pow( length( q / ( accretionDiskSize ) * 3.0 ), 2.0 ) * 10.0 );
q = p;
q.xz *= rot2D( iTime / 20.0 );
density += noise3D( vec3( length( q ) ) + q / 5.0 ) * StepSize * Absorption;
}
d = mapBH( p );
if( d < 0.0 )
{
col = vec3( 0.0 );
break;
}
rd = curve( rd, p );
}
col += density;
col *= illumination * 3.0 * StepSize;
return col;
}
void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
vec2 uv = ( fragCoord * 2.0 - iResolution.xy ) / iResolution.y;
vec3 col;
if( abs( uv.x ) < 1.2 && abs( uv.y ) < 0.9 )
{
col = rayTrace( uv );
} else {
col = vec3( 0.02 );
}
fragColor = vec4( col, 1.0 );
}
sry the disk is small. But thank you very very much for reading and have a great rest of your day.
GitHub - https://github.com/noarmsstuidos/Black-Hole-Rendering
Top comments (1)
It's been a while, good to be back.