// Shave and a Haircut
// (c) 2019 Epic Games
// US Patent 6720962

// ShaveHairBxdf - RIS bxdf for Shave hair.
//
// The hair bxdf models three specular transport paths:  
//   - R (reflection), 
//   - TRT (Transmission/Reflection/Transmission)
//   - TT (Transmission/Transmission)
//
// Specular component has been modified to match render images to those of
// the rsl version (i.e. Shave.sl).
// 
// Specular sampling part is using original PxrHair shader.
// 
// The diffuse component is based on:
// Goldman,
// "Fake Fur Rendering"
// SIGGRAPH 1997.
//
// Credits: Yosuke Katsura, Toneplus Animation Studios.

#include "RixShading.h" 
#include "RixBxdf.h" 
#include "RixIntegrator.h"
#include "RixRNG.h"
#include "RixShadingBuiltin.h"
#include "RixShadingUtils.h"

//#include "rx.h" // for RxNoise()

//#define ORIG_SPEC_CODE

static const unsigned char k_diffuseLobeId = 0;
static const unsigned char k_specularLobeId = 0;

static RixBXLobeSampled s_diffuseLobe;
static RixBXLobeSampled s_specularLobe;

static RixBXLobeTraits s_diffuseLobeTraits;
static RixBXLobeTraits s_specularLobeTraits;

inline RtColorRGB mix(const RtColorRGB& x, const RtColorRGB& y, const float alpha)
{
    return (x * (1.0f-alpha)) + (y * alpha);
}

inline float mix(float x, float y, float alpha)
{
    return (x * (1.0f-alpha)) + (y * alpha);
}

inline float clamp(float v, float min, float max)
{
    return (v < min) ? min : ((v < max) ? v : max);
}

inline RtColorRGB max(RtColorRGB x, RtColorRGB y)
{
    RtColorRGB max;
    max.r = (x.r > y.r) ? x.r : y.r;
    max.g = (x.g > y.g) ? x.g : y.g;
    max.b = (x.b > y.b) ? x.b : y.b;
    return max;
}

inline RtFloat luminance(const RtColorRGB& in)
{
    return (RtFloat)(0.299*in.r + 0.587*in.g + 0.114*in.b);
}

// stateless collection of functionality

//  Simple opacity and presence handler.
//
//  RenderMan calls these methods when it cannot determine presence or
//  opacity trivially and would have to do expensive shading computation to
//  determine them. So we want these to be faster than that.
//
//       us when opacity or presence services are required.
//       Owner should not instantiate us if the values convey
//       trivial opaque or present.
//     * GetPresence is invoked to when renderer wishes to skip
//       a more expensive shading computation.
//     * GetOpacity is invoked for shadows and must include presence.
//  
class SimpleOpacity : public RixOpacity 
{
public:
    SimpleOpacity(
        const RixShadingContext *ctx,
        RixBxdfFactory *bxdfFactory,
        const RtFloat *presence, 
        const RtColorRGB *transparency
    )
    : RixOpacity(ctx, bxdfFactory)
    , m_presence(presence)
    , m_transparency(transparency)
    {
    }

    //  Returns false if surface is trivially opaque.
    //  Otherwise it returns true and sets the opacity for each point
    //  being shaded.
    //
    //  Note that this method gets invoked for shadows as well.
    //
    virtual bool 
    GetOpacity(RtColorRGB *result)
    {
       bool ret = false;

       if (m_transparency)
       {
           for (int i = 0; i < shadingCtx->numPts; ++i)
           {
               result[i] = RixConstants::k_OneRGB - m_transparency[i];
               result[i].ClampAlbedo();
           }

           ret = true;
       }

       return ret;
   }

    //  Returns false if surface is trivially presence.
    //  Otherwise it returns true and sets the presence for each point
    //  being shaded.
    //
    virtual bool 
    GetPresence(RtFloat *result)
    {
        bool ret = false;

        if (m_presence)
        {
            for (int i = 0; i < shadingCtx->numPts; ++i)
            {
                result[i] = m_presence[i];

                if (!ret && (result[i] != 1.0f))
                {
                    ret = true;
                }
            }
        }

        return ret;
    }

private:
   const RtFloat *m_presence;
   const RtColorRGB *m_transparency;
};


class ShaveHairSpecular 
{
  public:

    ShaveHairSpecular(RixShadingContext const *sCtx,
                 RtFloat spec,
                 RtFloat glossiness,
                 RtColorRGB primarySpecColor, 
                 RtColorRGB secondarySpecColor, 
                 RtFloat const* index) 
                    // (to randomize hair highlights -- ignored) 
    : m_alphaR(0)
    , m_alphaTT(0)
    , m_alphaTRT(0)
    , m_betaR(0)
    , m_betaTT(0)
    , m_betaTRT(0)
    , m_glossiness(glossiness)
    , m_spec(spec)
    , m_primarySpecColor(primarySpecColor)
    , m_secondarySpecColor(secondarySpecColor)
    , m_bendTangentAmt(0.35f) // using same value as RSL version
    , m_secondarySHAVEglossMult(0.7f) // using same value as RSL version
    {
        RtVector3 const *Vn = NULL;
        RtVector3 const *Tn = NULL;
        RtInt nPts = sCtx->numPts;

#if 0
		RixMessages* m = (RixMessages*)sCtx->GetRixInterface(k_RixMessages);
		m->Info("ShaveHairSpecular: spec_color (%f, %f, %f), gloss %f", spec_color->r, spec_color->g, spec_color->b, Beta);
#endif
        RixShadingContext::Allocator pool(sCtx);
        m_B0n    = pool.AllocForBxdf<RtVector3>(nPts);
        m_B1n    = pool.AllocForBxdf<RtVector3>(nPts);
        m_thetaV = pool.AllocForBxdf<RtFloat>(nPts);
        m_phiV   = pool.AllocForBxdf<RtFloat>(nPts);
        RtVector3* orientation = pool.AllocForBxdf<RtVector3>(nPts);

        sCtx->GetBuiltinVar(RixShadingContext::k_Vn, &Vn);
        sCtx->GetBuiltinVar(RixShadingContext::k_Tn, &Tn);

        for (int i=0; i<nPts; ++i) {
            orientation[i] = Vn[i];

            // B0n perpendicular to orientation and curve direction
            m_B0n[i] = NormalizeCopy(Tn[i].Cross(orientation[i])); 
            // B1n in plane of Tn and orientation
            m_B1n[i] = NormalizeCopy(m_B0n[i].Cross(Tn[i]));  
            getThetaPhi(i, Tn[i], Vn[i], m_thetaV[i], m_phiV[i]);
        }

        // The following variables are used in sampling.
        // They came from original code 
        RtFloat Alpha = 7.5;   // highlightShift
        RtFloat Beta = 7.5;    // highlightWidth
        
        m_betaR = (F_PI/180.0f) * Beta;   // [5,10]
        m_betaTT = 0.5f * m_betaR;
        m_betaTRT = 3.0f * m_betaTT;

        m_alphaR = -(F_PI/180.0f) * Alpha;   // [-10,-5]
        m_alphaTT = -0.5f * m_alphaR;
        m_alphaTRT = -1.5f * m_alphaR;
    }
    
    static RixBXEvaluateDomain GetEvaluateDomain()
    {
        // Sample direct lighting on the sphere rather than the
        // hemisphere to ensure the TT lobe is included
        // in the direct lighting optimisation.
        return k_RixBXBoth; 
    }

    static RixBXLobeTraits GetAllLobeTraits()
    {
        return s_specularLobeTraits;
    }

    // utility convert a vector to two angles
    // Args:
    //	tn = input tangent vector at the sampling point
    //  dir = input vector to be converted to angles
    //  theta = output float azimuthal angle
    //  phi = output float altitude angle
    // * normal and binormal are from instance variables
    void getThetaPhi(int i, RtVector3 const & tn, RtVector3 const & dir,
                            RtFloat& theta, RtFloat& phi)
    {
        theta = F_PIDIV2 - acosf(dir.Dot(tn));    // [-Pi/2,Pi/2]
        phi = atan2f(dir.Dot(m_B0n[i]), dir.Dot(m_B1n[i])); // [-Pi,Pi]
    }

    // Utility  Add two angles in [-Pi,Pi] and get a result in [-Pi,Pi]
    static RtFloat addPhi(RtFloat phi1, RtFloat phi2)
    {
        RtFloat phiRet = phi1 + phi2; // [-2Pi, 2Pi]
        if(phiRet < -F_PI) phiRet += F_TWOPI;
        if(phiRet > F_PI) phiRet -= F_TWOPI;
        return phiRet;
    }


    RtColorRGB hspecular(int i,
			 const RtFloat angleTandL, 
			 const RtFloat angleTuandL,
			 const RtFloat angleTandV, 
			 const RtFloat angleTuandV)
    {
        RtFloat p = F_PI - angleTandL - angleTandV;
        RtFloat s = F_PI - angleTuandL - angleTuandV;
		
        RtFloat g = 1.0f/(0.101f-m_glossiness);
        RtFloat primarySpecPower = m_spec * powf(std::max(cosf(p), 0.0f), g) * F_INVPI; 
        g = 1.0f/(0.101f-m_secondarySHAVEglossMult*m_glossiness);
        RtFloat secondarySpecPower = m_spec * powf(std::max(cosf(s), 0.0f), g)* F_INVPI;

        return max(primarySpecPower * m_primarySpecColor,
                    secondarySpecPower * m_secondarySpecColor);
    }


    // generate a random phi in [-Pi, Pi]
    static RtFloat genPhi(RtFloat xi, RtFloat& pdf)
    {
        pdf = F_INVTWOPI; // 1/(2pi)
        return (xi - 0.5f) * F_TWOPI; //[-Pi, Pi]
    }

    // given a phi determine the probability that it would be picked
    static void evalPhi(RtFloat phi, RtFloat& pdf)
    {
        pdf = F_INVTWOPI; // 1/(2pi)
    }

    // We sample based on a Cauchy distribution from
    // "Importance Sampling for Hair Scattering", Ou et al.
    // To avoid rejecting samples we shift the single cauchy
    // distribution around based on phi to better target the
    // distributions we are drawing from.
    // Put the mapping from phi to alpha/beta in one place
    void getAlphaBeta(RtFloat phi, RtFloat& alpha, RtFloat& beta)
    {
        // distribute theta based on a (varying) gaussian
        RtFloat cosPhi2 = cosf(0.5f*phi);

        // invert so 0 means r and trt terms and so things stay near 0 longer
        RtFloat lerpTerm = 1.0f - cosPhi2;
        lerpTerm *= lerpTerm;
        lerpTerm *= lerpTerm; // ^4
        beta =  (1.0f - lerpTerm) * 0.75f * m_betaTRT + lerpTerm * m_betaTT;
        alpha = (1.0f - lerpTerm) * m_alphaR + lerpTerm * m_alphaTT;
    }

    // given a phi and a random val, generate theta
    RtFloat genTheta(int i, RtFloat xi, RtFloat phi, RtFloat& pdf)
    {
        RtFloat alpha, beta;
        getAlphaBeta(phi, alpha, beta);
        RtFloat shiftedThetaV = (0.5f * m_thetaV[i]) - alpha;
        RtFloat invThetaHMax = atanf((F_PIDIV4 + shiftedThetaV) / beta);
        RtFloat invThetaHMin = atanf((-F_PIDIV4 + shiftedThetaV) / beta);
        RtFloat thetaH = beta * tanf(xi *(invThetaHMax - invThetaHMin) +
                                     invThetaHMin);

        // half angle to angle gives a 1/2, the rest is the windowed
        // Cauchy distribution
        pdf = beta;
        pdf /= 2.0f * (invThetaHMax-invThetaHMin) * (thetaH*thetaH + beta*beta);

        RtFloat thetaNew = 2.0f * (thetaH + alpha) - m_thetaV[i];
        return clamp(thetaNew, -0.4999f*F_PI, 0.4999f*F_PI);
    }

    // given a theta and a phi, determine the conditional probability for theta
    void evalTheta(int i, RtFloat theta, RtFloat phi, RtFloat& pdf)
    {
        RtFloat alpha, beta;
        getAlphaBeta(phi, alpha, beta);
        RtFloat thetaH = 0.5f * (theta + m_thetaV[i]);
        thetaH -= alpha;
        RtFloat shiftedThetaV = (0.5f * m_thetaV[i]) - alpha;
        RtFloat invThetaHMax = atanf((F_PIDIV4 + shiftedThetaV) / beta);
        RtFloat invThetaHMin = atanf((-F_PIDIV4 + shiftedThetaV) / beta);

        // half angle to angle gives a 1/2, the rest is the windowed
        // Cauchy distribution
        pdf = beta;
        pdf /= 2.0f * (invThetaHMax-invThetaHMin) * (thetaH*thetaH + beta*beta);
    }

    // Generate specular sample
    void Generate(const RtVector3 &Vn,
                  const RtNormal3 &Nn,
                  const RtVector3 &dPdv,
                  const RtFloat2 &xi,
                  int i,
                  RtVector3 &On,
                  RtColorRGB &W,
                  RtFloat &FPdf, RtFloat &RPdf,
                  RixBXLobeSampled &lobeSampled)
    {
        // Tangents can be zero on fine hair tips. These cause invalid
        // direction vectors to be computed. We preempt this by setting the
        // material PDF to 0 which will disable further computation for the
        // sample
        RtVector3 Tn(dPdv.x, dPdv.y, dPdv.z);
        Tn.Normalize();
        if (Tn == RtFloat3(0.0f))
        {
            FPdf = RPdf = 0;
            return;
        }

        // generate phi uniformly around hair
        RtFloat phiPdf, thetaPdf;
        RtFloat phi = genPhi(xi.x, phiPdf);
        RtFloat phiL = addPhi(m_phiV[i], phi);
        RtFloat thetaL = genTheta(i, xi.y, phi, thetaPdf);

        // get sines and cosines
        RtFloat sinThetaL, cosThetaL, sinPhiL, cosPhiL;
        sinThetaL = sinf(thetaL);
        cosThetaL = cosf(thetaL);
        sinPhiL = sinf(phiL);
        cosPhiL = cosf(phiL);

        // convert into Ln
        On = (Tn * sinThetaL +
              m_B0n[i] * cosThetaL * sinPhiL +
              m_B1n[i] * cosThetaL * cosPhiL);

		RtVector3 V(Vn.x, Vn.y, Vn.z);
		V.Normalize();

		RtVector3 Tu = (1.0f-m_bendTangentAmt) * dPdv +  m_bendTangentAmt * -V;
		Tu.Normalize();

        // material response
		W = hspecular(i, acos(Tn.Dot(On)),
			    	 acos(Tu.Dot(On)), 
				 acos(Tn.Dot(V)),
				 acos(Tu.Dot(V)));

        // set the materialPdf for this sample.
        // compression near the poles (+- ctx->m_Tangent) gives 
        // 1/cos(thetaNew)
        FPdf = phiPdf * thetaPdf / cosThetaL;
        RPdf = FPdf;

        lobeSampled = s_specularLobe;
    }

    // Evaluate specular sample
    void Evaluate(const RtVector3 &Vn,
                  const RtVector3 &dPdv,
                  const RtVector3 &On,
                  int i,
                  RtColorRGB &W,
                  RtFloat &FPdf, RtFloat &RPdf)
    {
 		RtVector3 Tn(dPdv.x, dPdv.y, dPdv.z);
		Tn.Normalize();

        RtFloat thetaL, phiL;
        getThetaPhi(i, Tn, On, thetaL, phiL);
        RtFloat phi = addPhi(phiL, -m_phiV[i]);
        RtFloat phiPdf, thetaPdf;
        evalPhi(phi, phiPdf);
        evalTheta(i, thetaL, phi, thetaPdf);

        // material response

        RtVector3 V(Vn.x, Vn.y, Vn.z);
        V.Normalize();

        RtVector3 Tu = (1.0f-m_bendTangentAmt) * dPdv +  m_bendTangentAmt * -V;

        Tu.Normalize();
        W = hspecular(i,
                acos(Tn.Dot(On)),
                acos(Tu.Dot(On)),
                acos(Tn.Dot(V)),
                acos(Tu.Dot(V))) ;

        // compression near the poles (+- ctx->m_Tangent) gives 1/cos(theta)
        RtFloat sinTheta = On.Dot(Tn);
        RtFloat cosTheta = sqrtf(1.00001f - sinTheta * sinTheta);
        FPdf = phiPdf * thetaPdf / cosTheta;
        RPdf = FPdf;
    }

private:
    // PxrHair parameters
    // for pdf calculation
    RtFloat m_alphaR, m_alphaTT, m_alphaTRT;
    RtFloat m_betaR, m_betaTT, m_betaTRT;

    // shave parameters
    RtFloat m_glossiness;
    RtFloat m_spec;
    RtColorRGB m_primarySpecColor;
    RtColorRGB m_secondarySpecColor;
    RtFloat m_bendTangentAmt;
    RtFloat m_secondarySHAVEglossMult;

    // axes and angles
    RtVector3*  m_B0n;
    RtVector3*  m_B1n;
    RtFloat*    m_thetaV; 
    RtFloat*    m_phiV;
    RtFloat*    m_presence;
};

class GoldmanDiffuse
{
public:

    GoldmanDiffuse(RixShadingContext const *sCtx,
                   RtFloat diffuseGain,
                   RtFloat diffuseReflectGain,
                   RtFloat diffuseTransmitGain,
                   RtColorRGB const* diffuseRootColor,
                   RtColorRGB const* diffuseTipColor)
    {
        m_diffuseReflectGain = clamp(diffuseReflectGain, 0.f, 1.f);
        m_diffuseTransmitGain = clamp(diffuseTransmitGain, 0.f, 1.f);

        RtInt nPts = sCtx->numPts;
        RixShadingContext::Allocator pool(sCtx);
        m_diffuseRootColor = pool.AllocForBxdf<RtColorRGB>(nPts);
        m_diffuseTipColor  = pool.AllocForBxdf<RtColorRGB>(nPts);

        for (int i=0; i<nPts; ++i) {
            m_diffuseRootColor[i] = diffuseRootColor[i] * diffuseGain;
            m_diffuseTipColor[i] = diffuseTipColor[i] * diffuseGain; 
        }
  
        // here we assume hairs are parameterized by GSI-v with
        // 0 on root and 1 on tip.
        RtFloat const *v;
        sCtx->GetBuiltinVar(RixShadingContext::k_v, &v);
        RtColorRGB const *Cs;
        RixSCDetail detail;
        detail = sCtx->GetPrimVar("Cs", RtFloat3(1.0), (RtFloat3 const **)&Cs);
        bool varyingCs = (detail == k_RixSCVarying);

        m_C = pool.AllocForBxdf<RtColorRGB>(nPts);
        for (int i = 0; i < nPts; i++) 
        {
            m_C[i] = mix(m_diffuseRootColor[i], m_diffuseTipColor[i], v[i]);
            m_C[i] *= varyingCs ? Cs[i] : Cs[0];
            m_C[i] = max(m_C[i], RtColorRGB(0.0f));
        }
    }

    static RixBXEvaluateDomain GetEvaluateDomain() {
#if k_RixShadingVersion < 210
        return k_RixBXFront;
#else
		return k_RixBXReflect;
#endif
    }

    static RixBXLobeTraits GetAllLobeTraits() { 
        return s_diffuseLobeTraits;
    }

    // This is a first pass at Generate() for Goldman fake fur diffuse.
    // Stolen from PxrDiffuse.
    // To do: should generate sample according to Goldman fake fur diffuse term
    // or generate sample according to Lambertian (as now) but adjust W and 
    // pdfs for fake fur term?
    void Generate(const RtVector3 &Vn,
                  const RtNormal3 &Nn,
                  const RtVector3 &Tn,
                  const RtFloat2 &xi,
                  int index,
                  RtVector3 &On,
                  RtColorRGB &W,
                  RtFloat &FPdf, RtFloat &RPdf,
                  RixBXLobeSampled &lobeSampled)
    {
        RtVector3 Bn = Cross(Nn,Tn); // bitangent
        RtFloat cosTheta;
        RixCosDirectionalDistribution(xi, Nn, Bn, Tn, On, cosTheta);

        W = m_C[index] * cosTheta * F_INVPI;
        FPdf = cosTheta * F_INVPI;
        RPdf = fabsf(Dot(Vn, Nn)) * F_INVPI;
        if (RPdf > 0.0f) 
        {
            lobeSampled = s_diffuseLobe;
        }
        else
        {
            // rare case: at grazing angles V dot N can be zero or negative 
            lobeSampled.SetValid(false);
        }
    }

    // Evaluate Goldman fake fur diffuse term
    void Evaluate(const RtVector3 &Vn,
                  const RtNormal3 &Nn,
                  const RtVector3 &Tn,
                  const RtVector3 &On,
                  int index,
                  RtColorRGB &W,
                  RtFloat &FPdf, RtFloat &RPdf)
    {
        RtVector3 TcrossE = Tn.Cross(Vn);
        RtVector3 TcrossL = Tn.Cross(On);
        // k = cosine of the angle between the planes TxE and TxL
        float k = TcrossL.Dot(TcrossE); // (-1, 1)
        // directional attenuation factor
        float f_dir = 0.5f * ((1.0f+k) * m_diffuseReflectGain +
                              (1.0f-k) * m_diffuseTransmitGain);
        // cosine of the angle between T and L
        float TdotL = Tn.Dot(On);
        // Kajiya model is the sine between T and L:
        //  Kd * sin(T.L)
        //  sin = sqrt(1 - cos*cos)
        float kajiyaDiffuse = sqrtf(1.0f - TdotL*TdotL);

        // C contains Kd
        W = m_C[index] * f_dir * kajiyaDiffuse;

        FPdf = kajiyaDiffuse * F_INVPISQ;
        RPdf = FPdf;
    }

private:

    RtFloat m_diffuseReflectGain;
    RtFloat m_diffuseTransmitGain;
    RtColorRGB* m_diffuseRootColor;
    RtColorRGB* m_diffuseTipColor;
    RtColorRGB* m_C;
};

class ShaveHairBxdf : public RixBsdf
{
public:

    ShaveHairBxdf(RixShadingContext const* sc, RixBxdfFactory* bxf,
         GoldmanDiffuse* gdLobe,
         ShaveHairSpecular* hsLobe,
         RtNormal3 const* inputN)
        : RixBsdf(sc, bxf)
        , m_gdLobe(gdLobe)
        , m_hsLobe(hsLobe)
        , m_inputN(inputN)
    {
        m_lobesWanted = (gdLobe->GetAllLobeTraits() |
                         hsLobe->GetAllLobeTraits());
    }

    virtual RixBXEvaluateDomain GetEvaluateDomain() 
    { 
        return RixBXEvaluateDomain(ShaveHairSpecular::GetEvaluateDomain() |
                                   GoldmanDiffuse::GetEvaluateDomain());
    }

    virtual void GetAggregateLobeTraits(RixBXLobeTraits *t) 
    { 
        *t = m_lobesWanted;
    }

    //virtual bool EmitLocal(RtColorRGB* result)
    //{
    //}

    virtual void GenerateSample(RixBXTransportTrait transportTrait,
                                RixBXLobeTraits const *lobesWanted,
                                RixRNG *rng,
                                RixBXLobeSampled *lobeSampled,
                                RtVector3 *On, 
                                RixBXLobeWeights &W,
#if k_RixShadingVersion < 210
                                RtFloat *FPdf, RtFloat *RPdf
#else
								RtFloat *FPdf, RtFloat *RPdf,
								RtColorRGB* compTrans = NULL
#endif
								)
    {
        RtNormal3 const *Nn = m_inputN;
        RtVector3 const *Vn = NULL;
        RtVector3 const *dPdv = NULL;
        RtVector3 const *Tn = NULL;
        if (Nn == NULL) {
            shadingCtx->GetBuiltinVar(RixShadingContext::k_Nn, &Nn);
        }
        shadingCtx->GetBuiltinVar(RixShadingContext::k_Vn, &Vn);
        shadingCtx->GetBuiltinVar(RixShadingContext::k_Tn, &Tn);
        shadingCtx->GetBuiltinVar(RixShadingContext::k_dPdv, &dPdv);

        RtInt nPts = shadingCtx->numPts;
        RtFloat2 *xi = (RtFloat2*)alloca(sizeof(RtFloat2) * nPts);

        // Generate 2D sample points
        // needs to be out of main loop as it prevents vectorization
        rng->DrawSamples2D(nPts,xi);

        // Make any lobes that we may evaluate or write to active lobes, 
        // initialize their lobe weights to zero and fetch a pointer to the 
        // lobe weight arrays.
        // Important: AddActiveLobe() is a relatively cheap function, but it is 
        // expensive enough that it should be called here outside of the loop 
        // over the individual shading points.
        
        RtColorRGB *diffuseWgt = W.AddActiveLobe(s_diffuseLobe);
        RtColorRGB *specularWgt = W.AddActiveLobe(s_specularLobe);

#pragma ivdep
        //#pragma vector aligned
        for (int i = 0; i < nPts; i++) {
            RtFloat specChance;

            RixBXLobeTraits lobes = lobesWanted[i] & GetAllLobeTraits();
            bool doDiff = (lobes & s_diffuseLobeTraits).HasAny();
            bool doSpec = (lobes & s_specularLobeTraits).HasAny();
            bool doBoth = doDiff && doSpec;
            int lobe;

            if (doBoth) {
                specChance = 0.5; // 50% probability for each

				// pts[i].y already contains a random number between 0 and 1.
				// The call below uses that to select between the two lobes.
				// It then remaps that value to the full 0 to 1 range.
				//
				// For example, if pts[i].y is 0.3 then we'll get lobe 0
				// (specular). However, that would mean that all specular
				// samples would have y values < 0.5 while all diffuse
				// samples would have y values > 0.5. To prevent that y gets
				// remapped to 0.6.
                //
                // We could just generate a brand new random value in the
				// range 0 to 1, but this way we avoid additional calls to the
				// RNG, both for performance and to reduce the chances of
				// exhausting the random pool.
				//
                lobe = RixChooseAndRemap(xi[i].y, 1, &specChance);
            } else if (doSpec) {
                specChance = 1.0;
                lobe = 0;
            } else {
                specChance = 0.0;
                lobe = 1;
            }

            if (lobe == 1) { // choose Goldman diffuse reflection
                m_gdLobe->Generate(Vn[i], Nn[i], Tn[i], xi[i], i,
                                   On[i], diffuseWgt[i], FPdf[i], RPdf[i],
                                   lobeSampled[i]);
                FPdf[i] *= 1.0f - specChance;
                RPdf[i] *= 1.0f - specChance;
            } else { // lobe == 0: choose hair specular reflection (R/TT/TRT)
                m_hsLobe->Generate(Vn[i], Nn[i], dPdv[i], xi[i], i,
                                   On[i], specularWgt[i], FPdf[i], RPdf[i],
                                   lobeSampled[i]);
                FPdf[i] *= specChance;
                RPdf[i] *= specChance;
            }
        }
#if k_RixShadingVersion >= 210
		if (compTrans != NULL)
		{
			//	TODO
		}
#endif
    }

    //  0 chance of randomly drawing exactly the same vector.  
    virtual void EvaluateSample(RixBXTransportTrait transportTrait,
                                RixBXLobeTraits const *lobesWanted,
#if k_RixShadingVersion >= 210
								RixRNG* rng,
#endif
                                RixBXLobeTraits *lobesEvaluated,
                                const RtVector3 *On, RixBXLobeWeights &W,
                                RtFloat *FPdf, RtFloat *RPdf)
    { 
        RtNormal3 const *Nn = m_inputN;
        RtVector3 const *Vn = NULL;
        RtVector3 const *Tn = NULL;
        RtVector3 const *dPdv = NULL;

        if (Nn == NULL) {
            shadingCtx->GetBuiltinVar(RixShadingContext::k_Nn, &Nn);
        }
        shadingCtx->GetBuiltinVar(RixShadingContext::k_Vn, &Vn);
        shadingCtx->GetBuiltinVar(RixShadingContext::k_Tn, &Tn);
        shadingCtx->GetBuiltinVar(RixShadingContext::k_dPdv, &dPdv);
        RtInt nPts = shadingCtx->numPts;

        // Make any lobes that we may evaluate or write to active lobes, 
        // initialize their lobe weights to zero and fetch a pointer to the 
        // lobe weight arrays.
        // Important: AddActiveLobe() is a relatively cheap function, but it is 
        // expensive enough that it should be called here outside of the loop 
        // over the individual shading points.
        
        RtColorRGB *diffuseWgt = W.AddActiveLobe(s_diffuseLobe);
        RtColorRGB *specularWgt = W.AddActiveLobe(s_specularLobe);

#pragma ivdep
        //#pragma vector aligned
        for (int i = 0; i < nPts; i++) {
            RtFloat specChance;

            RixBXLobeTraits lobes = lobesWanted[i] & GetAllLobeTraits();
            bool doDiff = (lobes & s_diffuseLobeTraits).HasAny();
            bool doSpec = (lobes & s_specularLobeTraits).HasAny();
            bool doBoth = doDiff && doSpec;

            if (doBoth) {
                specChance = 0.5; // 50% probability for each
            } else if (doSpec) {
                specChance = 1.0;
            } else {
                specChance = 0.0;
            }

            FPdf[i] = RPdf[i] = 0.0f;
            lobesEvaluated[i].SetNone();

            if (doDiff) { // Goldman diffuse
                m_gdLobe->Evaluate(Vn[i], Nn[i], Tn[i], On[i], i,
                                   diffuseWgt[i], FPdf[i], RPdf[i]);
                FPdf[i] *= (1.0f - specChance);
                RPdf[i] *= (1.0f - specChance);
                lobesEvaluated[i] |= s_diffuseLobeTraits;
            }

            if (doSpec) { // hair specular (aka. glossy -- not dirac-specular)
                RtFloat specFPdf, specRPdf;
                m_hsLobe->Evaluate(Vn[i], dPdv[i], On[i], i,
                                   specularWgt[i], specFPdf, specRPdf);
                FPdf[i] += specChance * specFPdf;
                RPdf[i] += specChance * specRPdf;
                lobesEvaluated[i] |= s_specularLobeTraits;
            }
        }
    }

    virtual void EvaluateSamplesAtIndex(RixBXTransportTrait transportTrait,
                                        RixBXLobeTraits const &lobesWanted,
#if k_RixShadingVersion >= 210
										RixRNG* rng,
#endif
                                        RtInt index, RtInt nSamples,
                                        RixBXLobeTraits *lobesEvaluated,
                                        const RtVector3 *On,
                                        RixBXLobeWeights &W,
                                        RtFloat *FPdf, RtFloat *RPdf)
    {
        RtNormal3 const *Nn = m_inputN;
        RtVector3 const *Vn = NULL;
        RtVector3 const *dPdv = NULL;
        RtVector3 const *Tn = NULL;
        if (Nn == NULL) {
            shadingCtx->GetBuiltinVar(RixShadingContext::k_Nn, &Nn);
        }
        shadingCtx->GetBuiltinVar(RixShadingContext::k_Vn, &Vn);
        shadingCtx->GetBuiltinVar(RixShadingContext::k_Tn, &Tn);
        shadingCtx->GetBuiltinVar(RixShadingContext::k_dPdv, &dPdv);

        RixBXLobeTraits lobes = lobesWanted & GetAllLobeTraits();
        bool doDiff = (lobes & s_diffuseLobeTraits).HasAny();
        bool doSpec = (lobes & s_specularLobeTraits).HasAny();
        bool doBoth = doDiff && doSpec;
        RtFloat specChance; 

        if (doBoth) {
            specChance = 0.5; // 50% probability for each
        } else if (doSpec) {
            specChance = 1.0;
        } else {
            specChance = 0.0;
        }

        // Make any lobes that we may evaluate or write to active lobes, 
        // initialize their lobe weights to zero and fetch a pointer to the 
        // lobe weight arrays.
        // Important: AddActiveLobe() is a relatively cheap function, but it is 
        // expensive enough that it should be called here outside of the loop 
        // over the individual shading points.
        
        RtColorRGB *diffuseWgt = doDiff 
            ? W.AddActiveLobe(s_diffuseLobe) : NULL;
        RtColorRGB *specularWgt = doSpec 
            ? W.AddActiveLobe(s_specularLobe) : NULL;

#pragma ivdep
        //#pragma vector aligned
        for (int i = 0; i < nSamples; i++) {

            FPdf[i] = RPdf[i] = 0.0f;
            lobesEvaluated[i].SetNone();

            if (doDiff) { // Goldman diffuse
                m_gdLobe->Evaluate(Vn[i], Nn[i], Tn[i], On[i], i,
                                   diffuseWgt[i], FPdf[i], RPdf[i]);
                FPdf[i] *= (1.0f - specChance);
                RPdf[i] *= (1.0f - specChance);
                lobesEvaluated[i] |= s_diffuseLobeTraits;
            }

            if (doSpec) { // hair specular
                RtFloat specFPdf, specRPdf;
                m_hsLobe->Evaluate(Vn[i], dPdv[i], On[i], i,
                                   specularWgt[i], specFPdf, specRPdf);
                FPdf[i] += specChance * specFPdf;
                RPdf[i] += specChance * specRPdf;
                lobesEvaluated[i] |= s_specularLobeTraits;
            }
        }
    }

  private:
    GoldmanDiffuse* m_gdLobe;
    ShaveHairSpecular*   m_hsLobe;
    RixBXLobeTraits m_lobesWanted;
    RtNormal3 const* m_inputN;
};


// ShaveHairBxdf

class ShaveHairBxdfFactory : public RixBxdfFactory 
{
public:
    ShaveHairBxdfFactory();
    ~ShaveHairBxdfFactory();
    
    virtual int Init(RixContext &ctx, char const *pluginpath);
    RixSCParamInfo const *GetParamTable();
    virtual void Finalize(RixContext &ctx);

    virtual int CreateInstanceData(RixContext &ctx,
                                   char const *handle,
                                   RixParameterList const *plist,
                                   InstanceData *idata);

    virtual int GetInstanceHints(RtConstPointer instanceData) const;

    virtual void Synchronize(RixContext &ctx, RixSCSyncMsg syncMsg,
                             RixParameterList const *parameterList);

    virtual RixBsdf *BeginScatter(RixShadingContext const *ctx,
                                  RixBXLobeTraits const &lobesWanted,
                                  RixSCShadingMode sm,        
                                  RtConstPointer instanceData);
    virtual void EndScatter(RixBsdf *bsdf);

    virtual RixOpacity *BeginOpacity(RixShadingContext const *,
                                     RixSCShadingMode,
                                     RtConstPointer instancedata);
    virtual void EndOpacity(RixOpacity *); 

private:

    // default values for diffuse
    RtFloat    m_diffGain;
    RtFloat    m_diffReflGain;
    RtFloat    m_diffTransGain;
    RtColorRGB m_diffRootColor;
    RtColorRGB m_diffTipColor;

    RtFloat    m_index;

    // default values for presence
    RtFloat    m_presence;

	//	Default values for Shave parameters.
    RtColorRGB	m_Os;
    RtColorRGB	m_rootcolor;
    RtFloat	m_SHAVEambdiff;
    RtFloat	m_SHAVEgloss;
    RtFloat	m_SHAVEopacity;
    RtFloat	m_SHAVEselfshad;
    RtFloat	m_SHAVEspec;
    RtColorRGB	m_SHAVEspec_color;
    RtColorRGB	m_SHAVEspec_color2;
    RtColorRGB	m_tipcolor;
};

extern "C" PRMANEXPORT RixBxdfFactory *CreateRixBxdfFactory(const char *hint)
{
    return new ShaveHairBxdfFactory();
}

extern "C" PRMANEXPORT void DestroyRixBxdfFactory(RixBxdfFactory *factory)
{
    delete (ShaveHairBxdfFactory *) factory;
}

/*-----------------------------------------------------------------------*/
ShaveHairBxdfFactory::ShaveHairBxdfFactory()
{
    // diffuse
    m_diffGain = 1;
    m_diffReflGain = 1;
    m_diffTransGain = 1;
    m_diffRootColor = RtColorRGB(0.05f);
    m_diffTipColor = RtColorRGB(0.18f);

    // specular
    m_index = -1; // used to randomize hair highlights -- currently ignored

    // transmission/presence
    m_presence = 1.0;

	//	Initialize default values for Shave parameters.
    m_Os = RtColorRGB(0.0);
    m_rootcolor = RtColorRGB(1.0);
    m_SHAVEambdiff = 0.6f;
    m_SHAVEgloss = 0.07f;
    m_SHAVEopacity = 1.0f;
    m_SHAVEselfshad = 1.0f;
    m_SHAVEspec = 0.35f;
    m_SHAVEspec_color = RtColorRGB(1.0);
    m_SHAVEspec_color2 = RtColorRGB(1.0);
    m_tipcolor = RtColorRGB(1.0);
}

ShaveHairBxdfFactory::~ShaveHairBxdfFactory()
{
}

// Init
//  should be called once per RIB-instance. We look for parameter name
//  errors, and "cache" an understanding of our graph-evaluation requirements
//  in the form of allocation sizes.
int
ShaveHairBxdfFactory::Init(RixContext &ctx, char const *pluginpath)
{
    return 0;
}

// Synchronize: delivers occasional status information
// from the renderer. Parameterlist contents depend upon the SyncMsg.
// This method is optional and the default implementation ignores all 
// events.
void 
ShaveHairBxdfFactory::Synchronize(RixContext &ctx, RixSCSyncMsg syncMsg,
                             RixParameterList const *parameterList)
{
    if (syncMsg == k_RixSCRenderBegin)
    {
#if k_RixShadingVersion < 200
        s_diffuseLobe = RixBXLookupLobeByName(ctx, false, false, true, 
                                              k_diffuseLobeId, 
                                              "Diffuse");
        s_specularLobe = RixBXLookupLobeByName(ctx, false, true, true, 
                                               k_specularLobeId, 
                                               "Specular");
#else
        s_diffuseLobe = RixBXLookupLobeByName(ctx, false, false, true, false,
                                              k_diffuseLobeId, 
                                              "Diffuse");
        s_specularLobe = RixBXLookupLobeByName(ctx, false, true, true, false,
                                               k_specularLobeId, 
                                               "Specular");
#endif
        s_diffuseLobeTraits = RixBXLobeTraits(s_diffuseLobe);
        s_specularLobeTraits = RixBXLobeTraits(s_specularLobe);
    }
}

enum paramId
{
    k_diffuseGain=0,
    k_diffuseReflectGain,
    k_diffuseTransmitGain,
    k_diffuseRootColor,
    k_diffuseTipColor,
    k_transmitRootColor,
    k_transmitTipColor,
    k_index, // (ignored)
    // k_presence,
    k_inputN,
    k_inputAOV,

    // SHAVE PARAMS
    //
    k_Os,
    k_rootcolor,
    k_SHAVEambdiff,
    k_SHAVEgloss,
    k_SHAVEopacity,
    k_SHAVEselfshad,
    k_SHAVEspec,
    k_SHAVEspec_color,
    k_SHAVEspec_color2,
    k_tipcolor,

    k_numParams
};

RixSCParamInfo const *
ShaveHairBxdfFactory::GetParamTable()
{
    static RixSCParamInfo s_ptable[] = 
    {
        // diffuse inputs - const
        RixSCParamInfo("diffuseGain", k_RixSCFloat),
        RixSCParamInfo("diffuseReflectGain", k_RixSCFloat),
        RixSCParamInfo("diffuseTransmitGain", k_RixSCFloat),
        // diffuse inputs - connectable
        RixSCParamInfo("diffuseRootColor", k_RixSCColor),
        RixSCParamInfo("diffuseTipColor", k_RixSCColor),
        // specular inputs - connectable
        RixSCParamInfo("transmitRootColor", k_RixSCColor),
        RixSCParamInfo("transmitTipColor", k_RixSCColor),
        RixSCParamInfo("index", k_RixSCFloat), // (ignored)
        // transmission/presence
        //RixSCParamInfo("presence", k_RixSCFloat),
        RixSCParamInfo("inputN", k_RixSCNormal),
        RixSCParamInfo("inputAOV", k_RixSCInteger),

	//  SHAVE PARAMS
	//
        RixSCParamInfo("Os", k_RixSCColor),
        RixSCParamInfo("rootcolor", k_RixSCColor),
        RixSCParamInfo("SHAVEambdiff", k_RixSCFloat),
        RixSCParamInfo("SHAVEgloss", k_RixSCFloat),
        RixSCParamInfo("SHAVEopacity", k_RixSCFloat),
        RixSCParamInfo("SHAVEselfshad", k_RixSCFloat),
        RixSCParamInfo("SHAVEspec", k_RixSCFloat),
        RixSCParamInfo("SHAVEspec_color", k_RixSCColor),
        RixSCParamInfo("SHAVEspec_color2", k_RixSCColor),
        RixSCParamInfo("tipcolor", k_RixSCColor),

        RixSCParamInfo() // end of table
    };
    return &s_ptable[0];
}

// Finalize:
//  companion to Init, called with the expectation that any data 
//  allocated there will be released here.
void
ShaveHairBxdfFactory::Finalize(RixContext &)
{
}

RixBsdf *
ShaveHairBxdfFactory::BeginScatter(RixShadingContext const *sCtx,
                             RixBXLobeTraits const &lobesWanted,
                             RixSCShadingMode sm,
                             RtConstPointer instanceData)
{
    RixShadingContext::Allocator pool(sCtx);

    // XXX: add support for lobesWanted


	//	SHAVE PARAMS
	//
    RtFloat		const*	SHAVEambdiff;
    RtFloat		const*	SHAVEopacity;
    RtFloat		const*	SHAVEselfshad;
    RtColorRGB	const*	rootcolor;
    RtColorRGB	const*	tipcolor;

	//	Constant params.
	//
    sCtx->EvalParam(k_SHAVEambdiff, -1, &SHAVEambdiff, &m_SHAVEambdiff, false);
    sCtx->EvalParam(k_SHAVEopacity, -1, &SHAVEopacity, &m_SHAVEopacity, false);
    sCtx->EvalParam(k_SHAVEselfshad, -1, &SHAVEselfshad, &m_SHAVEselfshad, false);

	//	Varying params.
	//
	//	WARNING: OS, rootcolor and tipcolor are available by default but can
	//			 turned off in Shave Globals in Maya.
	//
    sCtx->EvalParam(k_rootcolor, -1, &rootcolor, &m_rootcolor, true);
    sCtx->EvalParam(k_tipcolor, -1, &tipcolor, &m_tipcolor, true);

    // Diffuse lobe requested
    GoldmanDiffuse* gdLobe = NULL;

    if (1) 
    {
        // Get diffuse data
        RtFloat const* diffGain;
        RtFloat const* diffReflGain;
        RtFloat const* diffTransGain;
        RtColorRGB const* diffRootColor;
        RtColorRGB const* diffTipColor;

        // constant inputs
        sCtx->EvalParam(k_diffuseGain, -1, &diffGain, &m_diffGain, false);
        sCtx->EvalParam(k_diffuseReflectGain, -1, &diffReflGain, 
                        &m_diffReflGain, false);
        sCtx->EvalParam(k_diffuseTransmitGain, -1, &diffTransGain, 
                        &m_diffTransGain, false);

        // varying inputs
        sCtx->EvalParam(k_diffuseRootColor, -1, &diffRootColor, 
                        &m_diffRootColor, true);
        sCtx->EvalParam(k_diffuseTipColor, -1, &diffTipColor, 
                        &m_diffTipColor, true);
    
        void* mem = pool.AllocForBxdf<GoldmanDiffuse>(1);

		// used Shave parameters for root and tip color
        gdLobe = new (mem) GoldmanDiffuse(sCtx, *diffGain, *diffReflGain,
                                *diffTransGain, rootcolor, tipcolor);
    }

    // Specular lobe requested
    ShaveHairSpecular* hsLobe = NULL;

    if (1) 
    {
        // Get specular data
        RtFloat const* spec;
       	RtFloat const* gloss;
       	RtColorRGB const* primarySpecColor;
       	RtColorRGB const* secondarySpecColor;
        RtFloat const* index; // ignored

        sCtx->EvalParam(k_SHAVEspec, -1, &spec, &m_SHAVEspec, false);
        sCtx->EvalParam(k_SHAVEgloss, -1, &gloss, &m_SHAVEgloss, false);
        sCtx->EvalParam(k_SHAVEspec_color, -1, &primarySpecColor, &m_SHAVEspec_color, false);
        sCtx->EvalParam(k_SHAVEspec_color2, -1, &secondarySpecColor, &m_SHAVEspec_color2, false);
        sCtx->EvalParam(k_index, -1, &index, &m_index, true);

        void *mem = pool.AllocForBxdf<ShaveHairSpecular>(1);
		hsLobe = new (mem) ShaveHairSpecular(sCtx, *spec, *gloss, 
				*primarySpecColor, *secondarySpecColor, index);
    }

    // input normal
    // only use if connected
    RixSCType type;
    RixSCConnectionInfo cinfo;
    RtNormal3 const* inputN = NULL;
    sCtx->GetParamInfo(k_inputN, &type, &cinfo);
    if (cinfo ==  k_RixSCNetworkValue) {
        sCtx->EvalParam(k_inputN, -1, &inputN, NULL, true);
    }

    // inputAOV plug
    // we only pull on the plug to trigger connected nodes.
    {
        RtInt const *aovPlug;
        RtInt foo = 0;
        sCtx->EvalParam(k_inputAOV, -1, &aovPlug, &foo, true);
    }

    void *mem = pool.AllocForBxdf<ShaveHairBxdf>(1);
    // Must use placement-new to set up the vtable properly
    ShaveHairBxdf *eval = new (mem) ShaveHairBxdf(sCtx, this, gdLobe, hsLobe, inputN);

    return eval;
}

void 
ShaveHairBxdfFactory::EndScatter(RixBsdf *)
{
    // since we know RixBsdf was placement newed and that it has
    // no special memory de-allactions duties, this is a no-op.    
}

// CreateInstanceData:
//    analyze plist to determine our response to GetOpacityHints.
int
ShaveHairBxdfFactory::CreateInstanceData(RixContext &ctx,
                                   char const *handle,
                                   RixParameterList const *plist,
                                   InstanceData *idata)
{
    RtUInt64 req = k_TriviallyOpaque | k_ComputesPresence | k_ComputesOpacity | k_OpacityCanBeCached;

    idata->data = (void *) req; // no memory allocated, overload pointer
    idata->freefunc = NULL;
    return 0;
}

int 
ShaveHairBxdfFactory::GetInstanceHints(RtConstPointer instanceData) const
{
    // our instance data is the RixBxdfFactory::InstanceHints bitfield.
    InstanceHints const &hints = (InstanceHints const&) instanceData;
    return hints;
}

RixOpacity *
ShaveHairBxdfFactory::BeginOpacity(RixShadingContext const *sCtx,
                             RixSCShadingMode shadingMode,
                             RtConstPointer instancedata)
{
    if(sCtx->scTraits.primaryHit == false) return NULL; 
    if(shadingMode != k_RixSCPresenceQuery) return NULL;
    
    // search Os
    RtColorRGB	const *oscolor;
    RixSCDetail detail;
    detail = sCtx->GetPrimVar("Os", RtFloat3(0.0), (RtFloat3 const **)&oscolor);
    
    // make it opaque if Os does not exist
    if(detail == k_RixSCInvalidDetail) return NULL;
   
    // use luminance of Os for presence
    RixShadingContext::Allocator pool(sCtx);
    RtInt nPts = sCtx->numPts;
    RtFloat *presence = NULL;
    presence = pool.AllocForBxdf<RtFloat>(nPts);
    for (int i = 0; i < nPts; i++) 
    {
		//	'presence' is the opacity of the sample with respect to the
		//	background behind it.
		//
        presence[i] = luminance(oscolor[i]); 
    }

    void *mem = pool.AllocForBxdf<SimpleOpacity>(1);
    
    RixOpacity *result = NULL;
    result = new (mem) SimpleOpacity(sCtx, this,
                                         presence, NULL /*transmitColor*/);
    return result;
}

void
ShaveHairBxdfFactory::EndOpacity(RixOpacity *)
{
    // since we know RixOpacity was placement newed and that it has
    // no special memory de-allactions duties, this is a no-op.    
}