I have always been fascinated with space phenonema. I decided to try to model a few of the beautiful objects in the night sky. The type of shaders that I wanted to create are far too complex to implement in hardware. Since I would have to use a software renderer, I decided to write my own in order to get the maximum possible speed.

Primitives

The first thing I needed were some primitive functions from which to build my shaders. I implemented some of the functions that are built into RenderMan.

noise()

I implemented Perlin's improved noise function [2]. His reference implementation in Java can be found here. I ported the code to C++. In his implementation, he chose the gradient associated with each lattice point using a series of if statements. The lattice coordinate is hashed to an integer. The lower bits of the hash code are used to choose the gradient.
float gradientMagnitude(int hash, float x, float y, float z)
{
int h = hash & 15;                     // CONVERT LO 4 BITS OF HASH CODE
float u = h<8||h==12||h==13 ? x : y,   // INTO 12 GRADIENT DIRECTIONS.
v = h<4||h==12||h==13 ? y : z;
return ((h&1) == 0 ? u : -u) + ((h&2) == 0 ? v : -v);
}
I found it to be slightly faster (and easier) to simply put all the possible gradients in a table and do a dot product.
float gradientMagnitude(int hash, float x, float y, float z)
{
static float g3[16][3] =
{
{ 1, 1, 0},{-1, 1, 0},{ 1,-1, 0},{-1,-1, 0}, // center of cube to edges
{ 1, 0, 1},{-1, 0, 1},{ 1, 0,-1},{-1, 0,-1},
{ 0, 1, 1},{ 0,-1, 1},{ 0, 1,-1},{ 0,-1,-1},
{ 1, 1, 0},{-1, 1, 0},{ 0,-1, 1},{ 0,-1,-1}  // tetrahedron
};

int h = hash & 0xf;
return x * g3[h][0] + y * g3[h][1] + z * g3[h][2];
}

I implemented a number of noise functions that take as parameters and return as values vectors of various dimensions. One way to implement a vector valued noise function is to call noise multiple times with different evaluation points. Since distant points are not correlated, neither are vector components:
Vector2 noise2( Vector2& p )
{
return Vector2( noise( P ), noise( P ) + Vector2( 1433.2, 623.3 ))
}
Another way to do this that takes less over head is to use the rest of bits in the hash code to select a second gradient for computing a second gradient magnitude.

Ken Musgrave suggests a variation on noise() that he calls vlnoise() - variable lacunarity noise [1, p. ???]. This noise function slightly perturbs the position of the evaluation point before calling noise. This gives the noise function a "swirlier" appearance. Here are two slices of fbm(), one created with noise() and the other with vlnoise().

 Truncated fbm() with noise() Truncated fbm() with vlnoise()

cellnoise()

Cellnoise associates a unique random number between 0 and 1 with each discrete unit voxel, or cell, in space. Steve Worley used a similar function for his cellular basis function [3]. My first attempt to create cellnoise() was based on the description in his paper. I hashed the cell coordinate (i,j,k) into an integer using the function he mentioned 541i+79j+31k mod 232 and used the result as the seed to a random number generator. The simplest random number generators are linear congruent generators (LCGs). An LCG is of the simple form:

yn+1 = yn * a + b mod m

The constants have to be chosen carefully in order to generate decent random numbers. I found a collection of LCGs with descriptions and analysis. I chose a generator that performs well and implemented it. I filled an image with random numbers to see what it looked like. This is what I got:

There are obviously some problems with this. It does not look random at all. It turns out that I had a subtle bug in my program. yn and a are 32 bit integers. When they are multiplied together the intermediate is 64 bits but it is truncated to 32 bits. Casting yn to an __int64 fixed the problem.

I ran into another interesting problem when playing around with LCGs. I wanted to create random numbers within a given integer range. You can do this simply taking the mod of the random number. I plotted the result of this, scaling the integer to the range 0 to 1. Here is the function:

unsigned LCG()
{
const int a = 1664525;
const int b = 1013904223;

yn = __int64(yn) * a + b;  // m = 2^32  - truncation performs the mod

return yn;
}

float getRandom()
{
return float(LCG() % 0x10000) / (0x10000 - 1);
}

Here is the result:

See the repeating pattern? It turns out that the lower bits of an number generated by an LCG are less random than the higher bits. This generator is particularly bad. It is the one suggested in Numerical Recipes in C, Second Edition, p. 284. For this reason they suggest scaling rather than using mod:

0x10000 * double(LCG()) / LCG_MAX;

After getting a decent random number generator working I tried the approach that Worley suggests. This is what I got:

It isn't random at all. The reason for this is that the hash isn't very good. Worley suggests that permutation tables would work better. But then you have only as many random number seeds as entries in the table. I decided to go with permutation tables and use a fixed table of random numbers instead of generating them. This generates cellnoise that is random enough for my purposes:

Searching Function Space - Endless Hours of Tweaking

After I had implemented the primitive functions I needed to come up with ways of combining them to get the kind of patterns that I wanted. In particular, I wanted to find something that produced billowy clouds with long, delicate structures to create gas clouds. This was a long and frustrating process of random tweaking. This screenshot gives you an idea of the number of parameters that there are to play with.

Here are some of the images I generated:

Finally, I figured out that using vector valued fbm() to perturb the evaluation point for another fbm() did just what I wanted.

Results

Starfields

Explicit creation of a star field is easy. You just plot a bunch of points according the desired distribution. With a procedural shader we need an implicit definition of the starfield. With the implicit definition it is possible to create as many stars as desired without explicitly storing them. One method that I have seen used before is to clip off the top of a high frequency noise function. This technique has several problems. The biggest problem is that the stars aren't round. It is also difficult to control the star spatial and brightness distributions. Another approach is to simply generate a random value at each pixel. If the number falls within a certain range return white, otherwise return black. This can generate a satisfactory pattern for a still image, but it will change when animated or even from run-to-run if the random number generated is seeded differently each time. The pattern can't be antialiased either because the points have no fixed spatial extent.  I have never seen a satisfactory star field shader. For this reason I wrote my own.

The key to a good starfield is to use a lot of cellnoise() to specify the density, brightnesses, and locations of the stars. Cellnoise ensures that the starfield is reproducible because the random values are linked to the evaluation points. To begin I determine which cell I am in simply by truncating the fractional part of the evaluation point. I call cellnoise() to determine how many stars are in the current cell:

starCount = cellnoise( cellP ) * maxStars;
I wrote a special version of cellnoise() that returns N random values. I call this function to return 3 * starCount random values. A triplet will be used to define the x and y position and brightness of each star. The color of the star depends on its brightness. The color is ramped from black to white until the brightness reaches the white point. Beyond the white point, the star's radius is increased and the color is maintained as white. This idea is based on the observation that in photographs of stars, the brightest stars appear to be larger than the others. If evaluation point falls within the star I return the star's color, otherwise I return black. Because the stars have finite support, stars from neighboring cells may overlap the current cell. These too must be examined. The distribution of stars and star brightnesses may be modified by biasing the cellnoise() toward 0 or 1.
starCount = bias( cellnoise( cellP ), starCountBias ) * maxStars;
Here are some images of starfields with different distributions.

 Uniform distribution - 5 stars per cell Star count  bias = 0.3.  Distribution is more uneven. Brightness bias = 0.2. Whitepoint < MaxBrightness. Variable radius stars.

I anti-alias the stars by subsampling around the evaluation point. I know ds and dt, the spatial extents of the pixel. I  determine how much of the star is covered by the pixel by checking the four points (x, y), (x + ds/2, y), (x, y + dt/2), and (x + ds/2, y + dt/2) to see if the lie within the star.

 Stars without anti-aliasing Stars with anti-aliasing. Stars look rounder.

By modulating the starCountBias you can have more explicit control over the star density. I used a noise function to modulate the density of the stars in these images:

 Stars and density function Click for a larger version

Galaxy

Here is a recipe for building a galaxy.

 Start with a base texture of scaled fbm(). Apply a warp to the texture to get a nice spiral. Add a truncated fbm() pattern to add some "clumpiness." Add a function with an exponential falloff to create the glowing center. Finally, add some stars (click to enlarge). Close up of the heart of the galaxy (click to enlarge).

I added a star glow around the stars that influences the gas clouds in front of them. You wouldn't see this in a real galaxy because if you can see individual stars, they are probably in front of the galaxy.

Code

The code is C++. The GUI was created with  Qt. The various parameters for the shader changed too often to modify the widgets by hand or in code. I wrote a generic form widget that takes a string of parameters and creates the text boxes and check boxes. The values of the of the parameters are retrieved by name. This flexible design made it easy to experiment.

References

[1] D. Ebert et al. Texturing and Modelling: A Procedural Approach, Academic Press, 1994.

[2] K. Perlin. Improving Noise. ACM Transactions on Graphics, 21 (3), pp. 681-682, 2002.

[3] S. Worley. A Cellular Texture Basis Function. Proceedings of SIGGRAPH 96, pp. 291-294, 1996.