The PdbEditor executes a pdb shader once for each particle contained in the pdb file(s) specified. Before executing the shader, the particle attributes have been updated to the values for the current pdb particle. Those values are accessed via the getAtt function, and can be modified via the setAtt function. For example
float radius;
getAtt("radiusPP", &radius);
retrieves the value of the particle radius, and places it into the variable radius. The names of the attributes, such as "radiusPP" in this example, are the ones contained in the Pdb file, and are generally based on the same names given by Maya or any other pdb generating package. Other attribute names commonly used are in the table below:
Attribute
Data Format
Description
radiusPP
float
radius
rgbPP
vector
color
position
vector
3D position of particle center
velocity
vector
3D velocity of particle
age
float
age of particle
id
int
particle id number
Probably to two most important attributes are the id and position. Because of that, the example editor shading routines below generally begin with the code segment
int id;
getAtt("id", &id);
vector position;
getAtt("position", &position);
The shading technology behind the pdb editor allows for a wide flexibility to manipulate particles, perform mathematical operations, compute vector mathematics, etc. The details can be found in its documentation. To some extent, many of the details can be learned just from following the examples below.
The shader function emit() performs application specific work within the editor. In the partman application, the emit() call sends the particle attributes the the renderer, which renders a sphere created according to those attributes. In other applications, the emit() call performs other functions.
A simple example editor shader using the emit() call is:
int main()
{
int id;
getAtt("id", &id);
vector position;
getAtt("position", &position);
/* Change the position just for fun */
vector displacement;
displacement[0] = 2.3;
displacement[0] = 0.6777;
displacement[0] = -45;
position = position + displacement;
setAtt("position", position);
/* Now emit the particle */
emit();
return 7;
}
This editor shader will simply shift every particle in the pdb file over by the vector amount (2.3, 0.6777, -45). As each one is shifted, emit() is called to dispose of it.
The reason why emit() is so important is because it may be called any number of times for each pdb "guide" particle. This allows us to turn one pdb particle into many "child" particles. As each child particle is created, the attributes are changed to suit the need, and emit() is called to dispose of it.
Here is a simple example of how to create child particles and emit them:
int main()
{
/* Retrieve guide particle info */
int id;
getAtt("id", &id);
srand48(id);
vector position;
getAtt("position", &position);
float radius;
getAtt("radiusPP", &radius);
/* Determine how many child particles are desired */
float childmax;
childmax = 100;
/* Reduce size of child particles
so that they fit inside guide particle */
float childradius;
childradius = 0.3 * radius / pow(childmax, 0.3333);
setAtt("radiusPP", childradius);
float child;
child = 0;
/* Loop over all children, creating new positions */
vector d;
while( child < childmax )
{
d[0] = drand48() - 0.5;
d[1] = drand48() - 0.5;
d[2] = drand48() - 0.5;
d = radius * d / sqrt( d . d );
d = d + position;
setAtt("position", d);
emit();
child = child + 1;
}
return 7;
}
The figures below show the effect of the child creation process using this code sample.
Guide particles
Child particles
In fact, the generic layout of the child creation editor can be outlined as follows:
/* Define global data structures for guide and child particles */
/* guide particle data */
int guideid;
float guideradius;
vector guideposition;
vector guidecolor;
vector guidevelocity;
float guideage;
/* child particle data */
float childcounter;
float maxchildcount;
float childradius;
vector childposition;
vector childcolor;
/* others... */
/* .... */
/* Define functions to do the work */
void getGuideParticleInfo()
{
getAtt("id", &guideid);
getAtt("radiusPP", &guideradius);
getAtt("position", &guideposition);
getAtt("rgbPP", &guidecolor);
getAtt("velocity", &guidevelocity);
getAtt("age", &guideage);
return;
}
void setChildParticleInfo()
{
/* dependent on application */
/* set child attributes that are fixed */
return;
}
void initChildParticle ()
{
childcounter = 0;
/* if using random numbers, initialize seed
to track with the guide particle */
srand48(guideid);
/* give inital values of changing attributes */
return;
}
int moreChildren ()
{
int flag;
flag = 0;
childcounter = childcounter + 1;
if(childcounter <= maxchildcount) {flag = 1;}
return flag;
}
int setChildParticle ()
{
/* Do something to distinquish this child
from the many others being created. */
return 1; /* return 1 to emit, 0 to not emit */
}
/* Operate main() in a fairly generic way */
int main()
{
/* Retrieve guide particle info of use */
getGuideParticleInfo();
/* Setup and Initialize child parameters */
setChildParticleInfo();
initChildParticle();
/* Loop over all children */
int emitdecision;
while( moreChildren() )
{
emitdecision = setChildParticle();
if(emitdecision==1){emit();}
}
return 7;
}
With this set-up, a variety of algorithms can be brought to bear to spread child particles around the volume of space based on pseudo-dynamical forms, guide particle data, statistics, and any other paradigm. The next section is devoted to providing several methods of filling the setChildParticle() function to achieve various ends. Many of these methods can be used simultaneously or in sequence.
Most of the algorithms in this section for placing child particles are based on some randomizing scheme. This is done in order to achieve complexity in the gross structure of accumulated particle collection. One of the primary issues in taking this approach is the control of the shape and size of the collection after all of the randomization has taken place.
Uniform space filling
In this approach, child particles are placed randomly within a specified volume, without history. By history, I mean that the placement of a child has no correlation with the placement of any other child. For example, uniformly filling a rectangular box has this form:
float boxheight, boxwidth, boxlength; /* values set by some criterion */
vector displacement;
int setChildParticle_UniformRectangleFill()
{
displacement[0] = boxheight * (drand48() - 0.5);
displacement[1] = boxwidth * (drand48() - 0.5);
displacement[2] = boxlength * (drand48() - 0.5);
childposition = guideposition + displacement;
setAtt("position", childposition);
return 1;
}
The figure below is an example of uniform rectangle filling.
Rectangular uniform filling
For filling a sphere, the process is only a little different:
float sphereradius; /* value set by some criterion */
vector displacement;
vector RandomUnitVector()
{
vector ruv;
float theta, phi;
theta = 3.14159265 * drand48();
phi = 2.0 * 3.14159265 * drand48();
ruv[0] = sin(theta) * cos(phi);
ruv[1] = sin(theta) * sin(phi);
ruv[2] = cos(theta);
return ruv;
}
int setChildParticle_UniformSphereFill()
{
displacement = RandomUnitVector();
displacement = sphereradius * displacement * drand48();
childposition = guideposition + displacement;
setAtt("position", childposition);
return 1;
}
The result is shown in the figure below.
Spherical uniform filling
To generate the images in this section, the setChildParticleInfo() routine has the form
void setChildParticleInfo()
{
/* dependent on application */
/* set child attributes that are fixed */
childradius = 0.2*guideradius;
setAtt("radiusPP", childradius);
maxchildcount = 100;
sphereradius = guideradius;
boxheight = guideradius;
boxwidth = 1.5 * guideradius;
boxlength = 2 * guideradius;
return;
}
Cauliflowers
A cauliflower is a hierarchy of spheres. Beginning from a guide sphere, a collection of smaller spheres are placed on its surface. If you desire to go to one more level, each of the spheres on the surface of the guide sphere can have a set of spheres placed on their surfaces. This process can continue as far as you like in a recursive way. The figure below illustrates one, two and three levels of sphere placement.
One recursion of cauliflower growth.
Two recursions of cauliflower growth.
Three recursions of cauliflower growth.
Four recursions of cauliflower growth.
Carrying this much farther would lead to a special object called Fractal Growth Pattern, which is also related to the random walk discussed in sections below.
For our algorithm, the new important parameters are:
/* cauliflower distribution */
float cauliflowerscale;
int nbcauliflowerclumps, nbcauliflowerrecursions;
The parameter cauliflowerscale controls the relative size of spheres in each recursion level. At each level, nbcauliflowerclumps spheres are placed on the surface of each sphere from the previous level. Finally, the total number of recursion levels is nbcauliflowerrecursions.
These three parameters are used to recursively create the cauliflower. The approach is:
void putbumps(vector gpos, float gradius, int nb, float scale, int rec)
{
int p, trackrec;
float theta, phi, radius;
vector pos;
/* track the number of recursions */
trackrec = rec-1;
/* loop over particles to clump onto surface of parent */
p = 0;
while(p < nb)
{
/* place the particle on the surface of its parent */
pos = RandomUnitVector();
pos = gpos + gradius * pos;
setAtt("position", pos);
/* give the particle a reduced size */
radius = gradius * scale;
setAtt("radiusPP", radius);
emit();
/* now recurse to put particles on the surface of this particle */
if(trackrec > 0){putbumps(pos, radius, nb, scale, trackrec);}
p = p + 1;
}
return;
}
int setChildParticle_Cauliflower()
{
emit(); /* emit guide particle */
/* recursively places and emits particles */
putbumps(guideposition, childradius, nbcauliflowerclumps,
cauliflowerscale, nbcauliflowerrecursions);
return 0; /* dont emit because emission occured during recursion */
}
Several new things are going on in this algorith. First, because of the recursion, all of the emit() calls take place inside setChildParticle_Cauliflower(), and its return value is zero so that no additional emit() calls are made. The second new technique is recursion. With each recursive call to putbumps(), the next level of spheres are placed, with a sphere radius that is larger or smaller than the previous level by a factor of cauliflowerscale. By setting this parameter less than one, the bumps get smaller and smaller. In the example images above, the values used were:
/* cauliflower placement */
cauliflowerscale = 0.5;
nbcauliflowerclumps = 20;
nbcauliflowerrecursions = 3;
Curves in space
Space curves are simply a parametric representation of the continuous sequence of points. We can draw the curve simply by placing particles in sequence along the length of the curve path. Positions on the curve are generated by knowing the position of an anchor point, and a parameter that traverses the lenght of the curve.
In the algorithm below, we generate a simple space curve which has constant helicity and curvature. For this, the important data to track is
/* Constant helicity and curvature space curve */
float pathlength, dpath, helicity, curvature;
vector Tangent, Normal, Binormal, T, N, B, G;
/* Space Curve Routines */
void ComputeSpaceCurveGamma()
{
float hh;
hh = 1.0/sqrt(1.0 + helicity*helicity);
G = (helicity * B - T ) * hh;
return;
}
vector ComputeSpaceCurvePosition()
{
vector Position;
float cl, sl, theta;
theta = pathlength * curvature;
sl = sin(theta);
cl = cos(theta);
float hh;
hh = 1.0/sqrt(1.0 + helicity*helicity);
Position = T * theta;
Position = Position - N * cl * hh;
Position = Position + G * (theta - sl) * hh;
Position = Position / curvature;
return Position;
}
void ComputeSpaceCurveTangent()
{
float cl, sl, theta;
theta = pathlength * curvature;
sl = sin(theta);
cl = cos(theta);
float hh;
hh = 1.0/sqrt(1.0 + helicity*helicity);
Tangent = T;
Tangent = Tangent + N * sl * hh;
Tangent = Tangent + G * (1.0-cl) * hh;
return;
}
void ComputeSpaceCurveNormal()
{
float cl, sl, theta;
theta = pathlength * curvature;
sl = sin(theta);
cl = cos(theta);
float hh;
hh = 1.0/sqrt(1.0 + helicity*helicity);
Normal = N * cl + G * sl;
return;
}
void ComputeSpaceCurveBinormal()
{
float cl, sl, theta;
theta = pathlength * curvature;
sl = sin(theta);
cl = cos(theta);
float hh;
hh = 1.0/sqrt(1.0 + helicity*helicity);
Binormal = B - ( N * sl + G * (1.0-cl) ) * helicity * hh;
return;
}
Space curves drawn from the guide particles.
Basic Random Walk
A random walk is basically just that. The basic of each successive child particle is a random step away from the position of the previous child particle. Over the space of many steps, the connected path can form a complex, fractally shape. Of course, there are many types of randomness, producing random walks with a variety of characteristics. We explore several of these in this and the next four algorithm sections.
The figure below was generated from the most basic random walk method. To accomplish this, the particles are very small (0.005 times the guide particle radius), and lots of them are used (20000 per guide particle). When even more particles are used in the next image (200,000 per guide particle, over 15 million total), you can clearly see that the random walk marches all over space.
The most basic random walk. There are 20000 child particles for each guide particle, over 1.5 million total.
Longer version of the basic random walk. There are 200000 child particles for each guide particle, over 15 million total.
The random walk process uses the following code:
vector RandomStep()
{
vector ruv;
ruv[0] = drand48() - 0.5;
ruv[1] = drand48() - 0.5;
ruv[2] = drand48() - 0.5;
return ruv;
}
int setChildParticle_BasicRandomWalk()
{
childposition = childposition + walkstep * RandomStep();
setAtt("position", childposition);
return 1;
}
The routine RandomStep() generates a vector with components lying in the range (-0.5,0.5). The current child particle is placed at a position that is randomly displaced from the previous particle, with the distance of the displacement in the range (0,sqrt(3)/2 walkstep).
There are no bounds on the size or length of random walk achievable. Because of the randomness in the direction and size of the steps however, there is a theoretical estimate of the approximate range the random walk will occupy. Based on the statistics of the random numbers being used, the root mean square step sizes is S = walkstep/sqrt(24). After maxchildcount steps, the range of the random walk should be roughly (S sqrt(maxchildcount)).
Correlated Random Walk
As might be expected, the correlated random walk is very similar to the basic random walk. The difference comes in one spot only: the random walk set is no longer independent from particle to particle. A statistical correlation is introduced in the following way:
int setChildParticle_CorrelatedRandomWalk()
{
walk = walk * mixin + walkstep * RandomStep() * mixout;
childposition = childposition + walk;
setAtt("position", childposition);
return 1;
}
The parameters mixin and mixout are mix the previous value of the random step with the new one, and so 0 < mixin < 1, and 0 < mixout < 1, with mixin + mixout = 1. The special case mixin = 0 is the uncorrelated basic random walk, while at the other extreme mixin = 1 produces purley straight lines.
The range in between uncorrelated random walk and straight line posseses a wide range of behaviors that we will try to organize in this section.
The first thing to sort out is what the impact of correlation is. The figure below shows that clearly: correlation turns the path that the particles are laid out on into a "smooth" curve with unpredictable twists and turns. As the mixin parameter approaches 1, the number of kinks and turns become fewer. The example below for example, using 50000 particles per guide particle, with short spacing between them. The mixin is 0.9999. Thats right, there are four 9's to the right of the decimal point.
A correlated random walk, with correlation mixin = 0.9999. Despite the very high correlation, there is considerable structure in each path. There are 50000 child particles for each guide particle.
To give you some idea of the importance of 0.9999 versus 0.999 for example, the image below was produced with mixin = 0.999. There is a substantial difference between the two.
A correlated random walk, with correlation mixin = 0.999. There are 50000 child particles for each guide particle.
Finally, below is the same set of particles as the other two, but with mixin = 0.99.
A correlated random walk, with correlation mixin = 0.99.
Spherical vs Cartesian walk
In the previous random walk examples, the vector step was of a "Cartesian" form, which means that the direction and length of the step were subject to randomization each step. This was accomplished by randomizing the three cartesian components of the step independently. An alternative useful method of conducting a random was is to randomize the direction of each step, but leave the magnitude of the step fixed. This is "spherical" walk. In practice, in requires only a small change in the code to accomplished. Using the random unit vector generator from the uniform spherical filling, the setChildParticle code is simply:
int setChildParticle_SphericalRandomWalk()
{
childposition = childposition + walkstep * RandomUnitVector();
setAtt("position", childposition);
return 1;
}
So, in the spherical random walk, successive steps are a distance exacly walkstep apart, whereas in the basic cartesian random walk, the distance between steps can be as little as 0 and as much as walkstep. This has the effect of leaving the spherical random walk volume less dense in the center, as shown in the example below.
A spherical (and uncorrelated) random walk.
Just as with the cartesian random walk, we can introduce correlation in the walk. This is accomplished this way:
int setChildParticle_CorrelatedSphericalRandomWalk()
{
walk = walk * mixin + walkstep * RandomUnitVector() * mixout;
childposition = childposition + walk;
setAtt("position", childposition);
return 1;
}
We have introduced the mixin and mixout parameters that same as in the basic cartesian case. An example of correlated spherical walk is below:
A spherical random walkwith correlation mixin = 0.999.
Levy Flight
A Levy Flight is a kind of random walk. The fundamental difference between Levy Flights and other random walks is that the length of the steps of a Levy Flight are chosen at random from a power law distribution. The impact of that choice is that the averge size of the walk after a number of steps is very different from other random walks. Levy Flights tend to have regions that look similar to ordinary random walks, but the regions are connects by occasional long single steps. This produces a nice clumping effect, as illustrated in the figure below.
A Levy Flight random walk. In this case, the parameter settings were LevyMu = 2.2 and LevyMin = 0.03 * guideradius. There was no correlation in the steps, and there were 2000 child particles (i.e. Levy steps) per guide particle.
The code for performing a Levy Flight is
vector RandomLevyStep()
{
vector ruv;
ruv = RandomUnitVector();
float step;
step = LevyMin * pow(drand48(), 1.0/(1.0-LevyMu));
ruv = step * ruv;
return ruv;
}
int setChildParticle_RandomLevyWalk()
{
childposition = childposition + RandomLevyStep();
setAtt("position", childposition);
return 1;
}
Directed walk
Walking the Surface of Implicit Surfaces
As discussed so far, random walks traverse space without any real restrictions on where they travel. For some purposes, it is useful to perform a random walk in order to get the complex structure, but limit the region the walk can occur it. Limiting the walk is the focus of this section.
Before proceeding, the convention chosen here for implicit surfaces is that the implicit function has a value of zero on the surface, has a positive value inside the surface, and a negative value outside the surface.
First, we need a method of computing the relevant information about the implicit surface. For this problem, we need routines that can retrieve two pieces of information:
For any point in space, what is the value of the implicit funcition.
For any point in space, what is the "normal" vector, defined as the unit vector that lines up with the gradient of the implicit function at that point, pointing outwards from the surface.
These two type of information will be contained in two function: getISFunction() and getISNormal(). For the simplest case of a spherical implicit surface, these can take the form
float getISFunction(vector r)
{
/* simple circle case */
vector test;
test = r - positionIS;
float returnvalue;
returnvalue = 1.0 - (test.test)/(radiusIS*radiusIS);
return returnvalue;
}
int getISSign(vector r)
{
/* > 0 => inside surface; < 0 => outside surface */
float ISFunctionvalue;
ISFunctionvalue = getISFunction(r);
int returnvalue;
returnvalue = 0;
if(ISFunctionvalue > 0){ returnvalue = 1; }
if(ISFunctionvalue < 0){ returnvalue = -1; }
return returnvalue;
}
vector getISNormal(vector r)
{
/* simple circle case */
vector test;
test = r - positionIS;
test = test / sqrt(test . test);
return test;
}
We have also added the routine getISSign() which uses getISFunction() to decide if a point in space is inside or outside the implicit surface. While we have built functions for spherical implicit surfaces, the process described here is applicable for any implicit surface.
The random walk illustrated in setChildParticle_CorrelatedTraverseIS() below is altered at each step to that is preferentially moves along the normal toward the surface of the Implicit Surface. At any step, the walk determines whether the current child position is inside or outside the surface, then re-orients the random step so that it moves toward the surface. There is still a random element in each step perpendicular to the normal.
int setChildParticle_CorrelatedTraverseIS()
{
float ISsign;
vector ISnormal;
ISnormal = getISNormal(childposition);
/* ISsign > 0 => inside; ISsign < 0 => outside */
ISsign = getISSign(childposition);
/* step in a random direction perpendicular to IS */
walk = RandomUnitVector();
walk = walk - (walk . ISnormal) * ISnormal;
walk = walk / sqrt(walk . walk);
walkIS = walkIS * mixin + stepIS * (ISsign * ISnormal + walk) * mixout;
childposition = childposition + walkIS;
setAtt("position", childposition);
return 1;
}
The two examples results below demonstrate that the random walk is bound to the implicit surface.
A random walk along the surface of an Implicit Surface (sphere), with mixin = 0.0.
The same random walk along the surface of an Implicit Surface as above, with mixin = 0.99.
Filling the inside (and outside) of Implicit Surfaces
It takes only a simple modification of the previous setChildParticle technique to have the random walk fill the inside of the implicit surface. The difference is that, when the particle is inside the implicit surface, the random walk is unaffected (whereas previously it was biased toward the surface). So now the routine setChildParticle_CorrelatedFillIS() looks like:
int setChildParticle_CorrelatedFillIS()
{
float ISsign;
/* ISsign > 0 => inside; ISsign < 0 => outside */
ISsign = getISSign(childposition);
/* step in a random direction perpendicular to IS */
walk = RandomUnitVector();
if(ISsign < 0) /* Outside: redirect inside */
{
vector ISnormal;
ISnormal = getISNormal(childposition);
walk = walk - (walk . ISnormal) * ISnormal;
walk = walk / sqrt(walk . walk);
walkIS = walkIS * mixin + stepIS * (ISsign * ISnormal + walk) * mixout;
}
else /* Inside: let it go */
{
walkIS = walkIS * mixin + stepIS * walk;
}
childposition = childposition + walkIS;
setAtt("position", childposition);
return 1;
}
As you can see, the difference is that now if the particle is inside, the random walk is unaltered.
The examples below demonstrate filling for various random walk correlations.
A random walk inside an Implicit Surface (sphere), with mixin = 0.0.
A random walk inside an Implicit Surface, with mixin = 0.5.
A random walk inside an Implicit Surface, with mixin = 0.9.
A random walk inside an Implicit Surface, with mixin = 0.95.
A random walk inside an Implicit Surface, with mixin = 0.99.
Time-dependence (pulses, shocks, waves, and cheesy dynamics)
The idea behind pulse is that a procedural method exists to tag positions along a random walk, and control any behavior with that tag. Because we have access to the frame number via the parameter$frame provided by the pdbEditor, the tag can be positions along the walk in a time dependent way.
An example of this is the following code:
float cosh(float arg)
{
float b;
b = exp(arg);
return 0.5*(b + 1.0/b);
}
float ComputeSpread()
{
float pl;
/* compute relative centroid of pulse based on speed and time */
pl = pathlength - pulse_spreadspeed * $frame;
/* use cosh function for a on-off behavior */
return pulse_minspread + (pulse_maxspread-pulse_minspread) / cosh(pl/pulse_spreadwavethickness);
}
int setChildParticle_Pulse()
{
/* Start out similar to a space curve */
childposition = guideposition + ComputeSpaceCurvePosition();
ComputeSpaceCurveNormal();
ComputeSpaceCurveBinormal();
ComputeSpaceCurveTangent();
pathlength = pathlength + dpath;
/* Now build a cloud around that position with a thickness */
float spread, angle;
spread = ComputeSpread();
int cloudsize, cloudcount;
cloudsize = pulse_cloudsize * (spread/pulse_minspread);
cloudcount = 0;
while(cloudcount < cloudsize)
{
displacement = RandomStep();
/* extract and scale component along tangent */
pulse_position = childposition + dpath * (displacement.Tangent)*Tangent;
/* extract and scale component perp to tangent */
displacement = displacement - (displacement.Tangent) * Tangent;
pulse_position = pulse_position + spread * displacement;
setAtt("position", pulse_position);
emit();
cloudcount = cloudcount + 1;
}
return 0;
}
What does this code do? The first few lines of setChildParticle_Pulse() just set up the walk as a SpaceCurve, exact as in the previous Curves in space section. Then a routine called ComputeSpread() is invoked. The number that comes from this is used in two ways. First, a set of "grandchild" particles are created, and the number of them depends on the spread value - the higher the spread, the more particles. Each grandchild is created as part of an ordinary random placement, but the placement is confined to the disk perpendicular to the space curve at that point. The distance that the placement extends perpendicular to the curve is controlled by the value of spread also. Effectively, spread controls the local thickness of the curve in space.
Now examine how spread is computed. The computation involves the inverse of the cosh() function, which is plotted below:
A plot of the 1/cosh() function.
So the ComputeSpread() routine returns the value pulse_minspread for all points on a walk except when pl = pathlength - pulse_spreadspeed * $frame nears a value of zero, where it smoothly changes to the value pulse_maxspread. But that near-zero point is a funciton of frame number, so over a sequence of frames, the segment of a walk with ComputeSpread() different from pulse_minspread moves farther out on the walk, giving the illusion of pulse propagation. This image below illustrates a frame of pulses in a set of space curves.
A set of particle strings with a pulse. As time progresses,
the pulse propagates along each string.
This effect can be modified so that instead of producing a limited region of a pulse, a shock-wave style effect can be produced. In this effect, once ComputeSpread() has smoothly changed from the value pulse_minspread to the value pulse_maxspread, it remain there. This is accomplished by the line
