//---------------------------------------------------------------------------
// solidHypertexture_voronoi.sl
// Author: Yukinori Inagaki
// Date: May 29 2008
//
// Reference: hypertexture.sl
//        _Advanced RenderMan: Creating CGI for Motion Picture_, 
//        by Anthony A. Apodaca and Larry Gritz, Morgan Kaufmann, 1999.
//
//---------------------------------------------------------------------------
  
  
#include "noises.h"
  
float volumedensity (point Pobj;  float blend, noisefreq, noisefreq2, stepsize;
                     output point cellcenter)
{
    float density = 0, density2 = 0;
    if(blend != 1)
    {
        voronoi_f1_3d (Pobj*noisefreq, 1, density, cellcenter);
        /* Increase Contrast */
           density = pow(clamp(density,0,1), 4);
    }
    if(blend != 0)
    {
        density2 = 0.5 + 0.5 * fBm(Pobj*noisefreq2, stepsize*noisefreq2, 3, 2, 0.6);
    }
    density = mix(density, density2, blend);
    return density;
}
  
normal compute_normal (point P; float den, blend, noisefreq, noisefreq2, stepsize)
{
    float density (point p) {
    point cellcenter;
    extern float blend, noisefreq, noisefreq2, stepsize;
    return volumedensity (p, blend, noisefreq, noisefreq2, stepsize, cellcenter);
    }
  
    normal norm = normal (density(P + vector (stepsize/10, 0, 0)) - den, 
               density(P + vector (0, stepsize/10, 0)) - den, 
               density(P + vector (0, 0, stepsize/10)) - den );
    return normalize(norm);
}
  
surface
solidHypertexture_voronoi (float Kd = 1;
          float end = 1.5;
          float stepsize = 0.001;
          float blend = 0.9;
          float noisefreq = 2.0;
          float noisefreq2 = 10;
          float thresh = 0.5)
{
    Ci = Oi = 0;
    /* Do not shade the back of the sphere -- only the front! */
    if (N.I < 0) 
    {
        point  Pobj = transform ("object", P);
        vector Iobj = normalize (vtransform ("object", I));
        float t0, t1;
    
        float d = stepsize;
        
        /* Calculate a reasonable step size */
        float ss = min (stepsize, end);
        point cellcenter;
        point Psamp = Pobj;
        float last_dtau = volumedensity (Psamp, blend, noisefreq, noisefreq2, stepsize, cellcenter);
        
        color result = 0;
        normal Npsamp;
        point  lastPnt = 0;
        if (last_dtau >= thresh){
            result = diffuse(N);
            Oi = 1;}
        else{
        while (d <= end && (comp(Oi,0)<1 || comp(Oi,1)<1 || comp(Oi,2)<1))
        {
            /* Take a step and get the density */
            ss = clamp (ss, 0.0001, end);
            d += ss;
            Psamp = Pobj + d*Iobj;
            float dtau = volumedensity (Psamp, blend, noisefreq, noisefreq2, stepsize, cellcenter);
            color li = 0;
            if (dtau > thresh)
            {
                Npsamp = compute_normal(Psamp, dtau, blend, noisefreq, noisefreq2, stepsize);
                
                result = diffuse(normalize(Npsamp));
                Oi = 1;
                break;
            }
        }
        }
        Ci = result*Kd;
    }
}