//---------------------------------------------------------------------------
// solidHypertexture_voronoi_bakeOcc.sl
// Author: Yukinori Inagaki
// Date: May 29 2008
//
// Description: This shader bakes point cloud for point-based occlusion.
//
// 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, sphererad, 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_bakeOcc (float Kd = 1;
          float end = 0.5;
          float stepsize = 0.1;
          float blend = 0.9;
          float noisefreq = 2.0;
          float noisefreq2 = 10;
          float thresh = 0.5;
          string bakefile = "")
{
    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;
        float a;
        if (last_dtau >= thresh)
        {
            result = diffuse(N);
            Oi = 1;
            if(bakefile != "")
            {
                a = area(Psamp, "dicing");
                bake3d(bakefile, "_area", P, N, "interpolate", 1, "_area", a);
            }
        }
        
        point last_Psamp = Psamp;
        while (d <= end)
        {
            /* 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);
            if (dtau > thresh - 0.04 && dtau < thresh + 0.04)
            {
                Npsamp = compute_normal(Psamp, dtau, blend, noisefreq, noisefreq2, stepsize);
                
                result = diffuse(normalize(Npsamp));
                Oi = 1;
                if(bakefile == "")
                    break;
                else
                {
                    a = area(Psamp, "dicing");
                    point pp = point "world" (Psamp[0],Psamp[1],Psamp[2]);
                    bake3d(bakefile, "_area", pp, Npsamp, "interpolate", 1, "_area", a);
                }
            }
            last_dtau = dtau;
            last_Psamp = Psamp;
        }
        Ci = result*Kd;    
    }
}