//---------------------------------------------------------------------------
// solidHypertexture_voronoi_occ.sl
// Author: Yukinori Inagaki
// Date: May 29 2008
//
// Description: This shader uses point cloud to calculate occlusion.
//
// Reference: hypertexture.sl
//        _Advanced RenderMan: Creating CGI for Motion Picture_, 
//        by Anthony A. Apodaca and Larry Gritz, Morgan Kaufmann, 1999.
//
//---------------------------------------------------------------------------
  
  
  
#include "noises.h"
//#include "normals.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);
}
  
float pointbasedocclusion (point pp; normal nn; string filename;
                        string hitsides; float maxdist, falloff, falloffmode,
                        samplebase, clampocclusion, maxsolidangle;) 
{
    //normal Ns = shadingnormal(nn);        
    float occ;
  
    occ = occlusion(pp, nn, 0, "pointbased", 1, "filename", filename,
                      "hitsides", hitsides,
                      "maxdist", maxdist, "falloff", falloff,
                      "falloffmode", falloffmode,
                      "samplebase", samplebase, "clamp", clampocclusion,
                      "maxsolidangle", maxsolidangle);
    occ *= 20;
  
    return 1 - occ;
}
  
surface
solidHypertexture_voronoi_occ (float Kd = 1;
          float end = 3;
          float stepsize = 0.1;
          float blend = 0.9;
          float noisefreq = 2.0;
          float noisefreq2 = 10;
          float thresh = 0.5;
          string bakedfile = "";
          string hitsides = "front";
          float maxdist = 1e15;
          float falloff = 0;
          float falloffmode = 0;
          float samplebase = 0;
          float clampocclusion = 0;
          float maxsolidangle = 0.05;)
{
    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;
        color white = 1;
        normal Npsamp;
        point  lastPnt = 0;
        if (last_dtau >= thresh){
            result = pointbasedocclusion (Psamp, N, bakedfile, hitsides,
                                    maxdist, falloff, falloffmode, samplebase,
                                    clampocclusion, maxsolidangle)* white;
            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, radius, stepsize);
                    
                    result = pointbasedocclusion (Psamp, Npsamp, bakedfile, hitsides,
                                    maxdist, falloff, falloffmode, samplebase,
                                    clampocclusion, maxsolidangle)* white;
                    Oi = 1;
                    break;
                }
            }
        }
        Ci = result*Kd;
    }
}