tebVoronoi.sl
/* 
 * Copyright (c) 1999 Thomas E. Burge.  All rights reserved.  
 * 
 *
 * FILE:  tebVoronoi.sl
 *
 *
 * DESCRIPTION:  
 *
 *    References:
 *          [WORL96]  Worley, Steven, A Cellular Texture Basis Function, 
 *                    SIGGRAPH 96 Conference Proceedings, Addison Wesley, 
 *                    1996, pp. 291-294.
 *
 */
#define DEBUG_


/* Note in the following macro, SwapP and SwapF can conflict with the
 * variable names defined in the scope that the Swap() macro was called 
 * from.  This problem showed up with the PRMan 3.8c compiler.  Also 
 * if Swap() is setup as a function, the parameters would all be of 
 * output type.  Array elements (featurepoints[2] for example) do not
 * work as output parameters.  Inshort Swap is a macro here because
 * the data was not getting to the array elements that needed to be 
 * swapped.
 */
#define Swap( p1, p2, f1, f2 ) \
{                              \
   point  SwapP;               \
   float  SwapF;               \
                               \
   SwapP  = p1;                \
   SwapF  = f1;                \
   p1 = p2;                    \
   f1 = f2;                    \
   p2 = SwapP;                 \
   f2 = SwapF;                 \
}


void FindFeaturePoints( float n; point P; float i; float j; float k; 
/* NOTE PRMan can take the following without the size in the array:
 *    output point featurepoints[];
 *    output float distances[];
 */
                        output point featurepoints[9];
                        output float distances[9]; )
{
   float  seed, npoints, dist;
   point  p;
   float  s;
   float  c;


   s = 0;
   seed = 541*i + 79*j + mod(31*k,999983);
   npoints = floor( 5 * cellnoise(s,seed) ) + 4;
   npoints *= 3; /* Three coordinates (x,y,z) per point. */
   for ( s=1; s < npoints; s+=3 ) {
      p = point( i+cellnoise(s,seed),
		 j+cellnoise(s+1,seed),
		 k+cellnoise(s+2,seed) );
      dist = distance( P, p );

      /* I wanted two for loops with a break in the first, but found
       * that the PRMan 3.8 compiler makes the break act as a return.
       */
      for ( c=0; c < n; c+=1 ) {
	 if ( dist <= distances[c] )
	    Swap( p, featurepoints[c], dist, distances[c] );
      }
   }   
}


#define DIST 0.5
surface tebVoronoi( uniform float Ka = 1.0; 
		    uniform float Kd = .5; 
		    /* cellbasistype can be assigned one of the following:
		     *   "voronoi", "medialaxis", "F1", "F2", "F3", "F4", 
		     *   or "edgedvoronoi". 
		     * voronoi      == F1
		     * medialaxis   == F2 -- has a medial-axis diagram "look"
		     * edgedvoronoi == F2 - F1
		     */
		    uniform string cellbasistype = "voronoi";
		    )
{
   uniform float  nF;
   float   x,y,z;
   float   i,j,k;
   float   xn,yn,zn;
   float   distances[9] = {1e5,1e5,1e5,1e5,1e5,1e5,1e5,1e5,1e5};
   point   featurepoints[9] = {0,0,0,0,0,0,0,0,0};
   normal  Nf;
   color   Ct;
   

   if ( cellbasistype=="voronoi" || cellbasistype=="F1" )
      nF = 1;
   else if ( cellbasistype=="medialaxis" || cellbasistype=="F2" 
	     || cellbasistype=="edgedvoronoi")
      nF = 2;
   else if ( cellbasistype=="F3" )
      nF = 3;
   else if ( cellbasistype=="F4" )
      nF = 4;

   x = xcomp(P);
   y = ycomp(P);
   z = zcomp(P);
   i = floor(x);
   j = floor(y);
   k = floor(z);

   xn = x - i;
   yn = y - j;
   zn = z - k;

   FindFeaturePoints( nF, P, i,   j,   k,   featurepoints,distances );

   /* Checking bounds gives a huge performance boost, but if nF gets
    *    greater than 4 when there are 9 points per cell the DIST bounds
    *    needs to be lowered. 
    */
#if 1
   if ( zn <= DIST )
   {
      FindFeaturePoints( nF, P, i,   j,   k-1, featurepoints,distances );
      if ( yn <= DIST )
      {
	 FindFeaturePoints( nF, P, i,   j-1, k-1, featurepoints,distances );
	 FindFeaturePoints( nF, P, i,   j-1, k,   featurepoints,distances );
      }
      else
      {
	 FindFeaturePoints( nF, P, i,   j+1, k-1, featurepoints,distances );
	 FindFeaturePoints( nF, P, i,   j+1, k,   featurepoints,distances );
      }
      if ( xn <= DIST )
      {
	 FindFeaturePoints( nF, P, i-1, j,   k-1, featurepoints,distances );
	 FindFeaturePoints( nF, P, i-1, j,   k,   featurepoints,distances );
	 if ( yn <= DIST )
	 {
	    FindFeaturePoints( nF, P, i-1, j-1, k-1, featurepoints,distances );
	    FindFeaturePoints( nF, P, i-1, j-1, k,   featurepoints,distances );
	 }
	 else
	 {
	    FindFeaturePoints( nF, P, i-1, j+1, k-1, featurepoints,distances );
	    FindFeaturePoints( nF, P, i-1, j+1, k,   featurepoints,distances );
	 }
      }
      else
      {
	 FindFeaturePoints( nF, P, i+1, j,   k-1, featurepoints,distances );
	 FindFeaturePoints( nF, P, i+1, j,   k,   featurepoints,distances );
	 if ( yn <= DIST )
	 {
	    FindFeaturePoints( nF, P, i+1, j-1, k-1, featurepoints,distances );
	    FindFeaturePoints( nF, P, i+1, j-1, k,   featurepoints,distances );
	 }
	 else
	 {
	    FindFeaturePoints( nF, P, i+1, j+1, k-1, featurepoints,distances );
	    FindFeaturePoints( nF, P, i+1, j+1, k,   featurepoints,distances );
	 }
      }
   }
   else
   {
      FindFeaturePoints( nF, P, i,   j,   k+1, featurepoints,distances );
      if ( yn <= DIST )
      {
	 FindFeaturePoints( nF, P, i,   j-1, k+1, featurepoints,distances );
	 FindFeaturePoints( nF, P, i,   j-1, k,   featurepoints,distances );
      }
      else
      {
	 FindFeaturePoints( nF, P, i,   j+1, k+1, featurepoints,distances );
	 FindFeaturePoints( nF, P, i,   j+1, k,   featurepoints,distances );
      }
      if ( xn <= DIST )
      {
	 FindFeaturePoints( nF, P, i-1, j,   k+1, featurepoints,distances );
	 FindFeaturePoints( nF, P, i-1, j,   k,   featurepoints,distances );
	 if ( yn <= DIST )
	 {
	    FindFeaturePoints( nF, P, i-1, j-1, k+1, featurepoints,distances );
	    FindFeaturePoints( nF, P, i-1, j-1, k,   featurepoints,distances );
	 }
	 else
	 {
	    FindFeaturePoints( nF, P, i-1, j+1, k+1, featurepoints,distances );
	    FindFeaturePoints( nF, P, i-1, j+1, k,   featurepoints,distances );
	 }
      }
      else
      {
	 FindFeaturePoints( nF, P, i+1, j,   k+1, featurepoints,distances );
	 FindFeaturePoints( nF, P, i+1, j,   k,   featurepoints,distances );
	 if ( yn <= DIST )
	 {
	    FindFeaturePoints( nF, P, i+1, j-1, k+1, featurepoints,distances );
	    FindFeaturePoints( nF, P, i+1, j-1, k,   featurepoints,distances );
	 }
	 else
	 {
	    FindFeaturePoints( nF, P, i+1, j+1, k+1, featurepoints,distances );
	    FindFeaturePoints( nF, P, i+1, j+1, k,   featurepoints,distances );
	 }
      }
   }
#else
   FindFeaturePoints( nF, P, i,   j-1, k-1, featurepoints, distances );
   FindFeaturePoints( nF, P, i,   j,   k-1, featurepoints, distances );
   FindFeaturePoints( nF, P, i,   j+1, k-1, featurepoints, distances );
   FindFeaturePoints( nF, P, i-1, j-1, k-1, featurepoints, distances );
   FindFeaturePoints( nF, P, i-1, j,   k-1, featurepoints, distances );
   FindFeaturePoints( nF, P, i-1, j+1, k-1, featurepoints, distances );
   FindFeaturePoints( nF, P, i+1, j-1, k-1, featurepoints, distances );
   FindFeaturePoints( nF, P, i+1, j,   k-1, featurepoints, distances );
   FindFeaturePoints( nF, P, i+1, j+1, k-1, featurepoints, distances );

   FindFeaturePoints( nF, P, i,   j-1, k,   featurepoints, distances );
   FindFeaturePoints( nF, P, i,   j+1, k,   featurepoints, distances );
   FindFeaturePoints( nF, P, i-1, j-1, k,   featurepoints, distances );
   FindFeaturePoints( nF, P, i-1, j,   k,   featurepoints, distances );
   FindFeaturePoints( nF, P, i-1, j+1, k,   featurepoints, distances );
   FindFeaturePoints( nF, P, i+1, j-1, k,   featurepoints, distances );
   FindFeaturePoints( nF, P, i+1, j,   k,   featurepoints, distances );
   FindFeaturePoints( nF, P, i+1, j+1, k,   featurepoints, distances );

   FindFeaturePoints( nF, P, i,   j-1, k+1, featurepoints, distances );
   FindFeaturePoints( nF, P, i,   j,   k+1, featurepoints, distances );
   FindFeaturePoints( nF, P, i,   j+1, k+1, featurepoints, distances );
   FindFeaturePoints( nF, P, i-1, j-1, k+1, featurepoints, distances );
   FindFeaturePoints( nF, P, i-1, j,   k+1, featurepoints, distances );
   FindFeaturePoints( nF, P, i-1, j+1, k+1, featurepoints, distances );
   FindFeaturePoints( nF, P, i+1, j-1, k+1, featurepoints, distances );
   FindFeaturePoints( nF, P, i+1, j,   k+1, featurepoints, distances );
   FindFeaturePoints( nF, P, i+1, j+1, k+1, featurepoints, distances );
#endif

#if DEBUG
   printf("P:%p\n",P);
   for (i=0;i < N;i+=1)
      printf("d[%f]:%f p:%p\n",i,distances[i],featurepoints[i]);
#endif
   /* Note setting Ct to equal distances[0] gives a voronoi texmap. 
    * Setting Ct to distances[1] gives a medial axis texmap. 
    * Using (distances[1]-distances[0]) for Ct gives a voronoi texmap
    * with a thickness given to the voronoi diagram's edges.
    */
   if ( cellbasistype=="voronoi" || cellbasistype=="F1" )
      Ct = distances[0];
   else if ( cellbasistype=="medialaxis" || cellbasistype=="F2" )
      Ct = distances[1];
   else if ( cellbasistype=="F3" )
      Ct = distances[2];
   else if ( cellbasistype=="F4" )
      Ct = distances[3];
   else if ( cellbasistype=="edgedvoronoi" )
      Ct = distances[1] - distances[0];
   else
      Ct = 0.0;

   Nf = faceforward( normalize(N), I );
   
   Ci = Os * ( Cs * Ct * (Ka*ambient()+Kd*diffuse(Nf)) );
   Oi = Os;
}