RixShadingUtils.h

Go to the documentation of this file.

/*
# ------------------------------------------------------------------------------
#
# Copyright (c) 2020 Pixar. All rights reserved.
#
# The information in this file (the "Software") is provided for the exclusive
# use of the software licensees of Pixar ("Licensees").  Licensees have the
# right to incorporate the Software into other products for use by other
# authorized software licensees of Pixar, without fee. Except as expressly
# permitted herein, the Software may not be disclosed to third parties, copied
# or duplicated in any form, in whole or in part, without the prior written
# permission of Pixar.
#
# The copyright notices in the Software and this entire statement, including the
# above license grant, this restriction and the following disclaimer, must be
# included in all copies of the Software, in whole or in part, and all permitted
# derivative works of the Software, unless such copies or derivative works are
# solely in the form of machine-executable object code generated by a source
# language processor.
#
# PIXAR DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING ALL
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL PIXAR BE
# LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
# WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION
# OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
# CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.  IN NO CASE WILL
# PIXAR'S TOTAL LIABILITY FOR ALL DAMAGES ARISING OUT OF OR IN CONNECTION WITH
# THE USE OR PERFORMANCE OF THIS SOFTWARE EXCEED $50.
#
# Pixar
# 1200 Park Ave
# Emeryville CA 94608
#
# ------------------------------------------------------------------------------
*/

#ifndef RixShadingUtils_h
#define RixShadingUtils_h

// RixShadingUtils:
//
// Utility/portability functions for RixShading plugins.
//

#include <algorithm>                    // for max
#include <cassert>                      // for assert
#include <cmath>                        // for sqrt, cosf, sinf, floorf, etc
#include <cstddef>                      // for NULL
#include "RixInterfaces.h"              // for RixRenderState, etc
#include "RixShading.h"                 // for RixShadingContext, etc
#include "prmanapi.h"                   // for PRMAN_INLINE
#include "RiTypesHelper.h"              // for RtVector3, RtColorRGB, etc
#include "RixPredefinedStrings.hpp"     // for k_N primvar

#if defined(WIN32)
#include <malloc.h>
#include <intrin.h>
#else
#include <alloca.h>
#endif


// single-precision floating point constants
// clang-format off
#define F_PI            (3.14159265f)
#define F_TWOPI         (6.283185307f)
#define F_FOURPI        (12.56637061f)
#define F_INVPI         (0.318309886f)
#define F_INVPISQ       (0.1013211836f)
#define F_INVTWOPI      (0.15915494309f)
#define F_INVFOURPI     (0.0795774715f)
#define F_PIDIV2        (1.57079632679f)
#define F_PIDIV4        (0.785398163397f)
#define F_DEGTORAD      (0.017453292520f)
#define F_RADTODEG      (57.2957795131f)
#define F_INVLN2        (1.44269504089f)
#define F_INVLN4        (0.72134752f)
#define F_LOG2E         (1.44269504088896f)
#define F_SQRT2         (1.41421356237f)
#define F_MAXDIST       (1.0e10f)
#define F_SQRTMIN       (1.0842021721e-19f)
#define F_MINDIVISOR    (1.0e-6f)
// clang-format on

namespace RixConstants
{

// some additional constants

static const RtColorRGB k_ZeroRGB(0.0f);
static const RtColorRGB k_OneRGB(1.0f);
static const RtFloat3 k_ZeroF3(0.0f);
static const RtFloat3 k_OneF3(1.0f);
static const RtFloat2 k_ZeroF2(0.0f);
static const RtFloat2 k_OneF2(1.0f);
static const float k_ZeroF(0.0f);
static const float k_OneF(1.0f);
static const float k_Infinity(1e+38f);
static const RtMatrix4x4 k_IdentityMatrix(1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1);

} // namespace RixConstants


//
// Returns the zero-based index of the first (least significant) bit set in
// a 64-bit integer.  Compiles to a single quick instruction via intrinsics.
//
// NOTE: The result is undefined if the input is zero!
//
PRMAN_INLINE unsigned int RixFindFirstSetBit(unsigned long long v)
{
#if defined(WIN32)
    unsigned long r = 0;
    _BitScanForward64( &r, v );
    return r;
#else // Linux, OS X, gcc or clang
    return __builtin_ctzll( v );
#endif
}


// Converts degrees to radians
PRMAN_INLINE float
RixDegreesToRadians(float degrees)
{
    return degrees * F_DEGTORAD;
}


// Converts radians to degrees
PRMAN_INLINE float
RixRadiansToDegrees(float rads)
{
    return rads * F_RADTODEG;
}


template <typename T>
PRMAN_INLINE T
RixSignum(T x)
{
    if (x > T(0))
        return T(1);
    if (x < T(0))
        return T(-1);
    return T(0);
}


//
// RixIsFinite is helpful in diagnosing bugs that produce NaNs
// and other exotic breads.
//
PRMAN_INLINE int
RixIsFinite(float x)
{
#if defined(WIN32)
    return (x == x) && (x != std::numeric_limits<float>::infinity()) &&
           (x != -std::numeric_limits<float>::infinity());
#else
    return std::isfinite(x);
#endif
}


//
// Returns true of all three components of a point/vector/normal are
// finite
//
PRMAN_INLINE int
RixIsFinite(RtFloat3 const& v)
{
    return RixIsFinite(v.x) && RixIsFinite(v.y) && RixIsFinite(v.z);
}


//
// Returns the minimum of two floats
//
PRMAN_INLINE float
RixMin(float x, float y)
{
    return x < y ? x : y;
}


//
// Returns the minimum of two integers
//
PRMAN_INLINE int
RixMin(int x, int y)
{
    return x < y ? x : y;
}


//
// Returns the maximum of two floats
//
PRMAN_INLINE float
RixMax(float x, float y)
{
    return x > y ? x : y;
}


//
// Returns the maximum of two integers
//
PRMAN_INLINE int
RixMax(int x, int y)
{
    return x > y ? x : y;
}


//
// Clamps float x in the interval [min,max]
//
PRMAN_INLINE float
RixClamp(float x, float min, float max)
{
    return x < min ? min : (x > max) ? max : x;
}


//
// Clamps integer x in the interval [min,max]
//
PRMAN_INLINE int
RixClamp(int x, int min, int max)
{
    return x < min ? min : (x > max) ? max : x;
}


//
// Returns the fractional part of float x
//
PRMAN_INLINE float
RixFractional(float x)
{
    return (x - floorf(x));
}


//
// RixMod is a straight port of RSL's mod() shadeop.
// It is NOT equivalent to fmodf() from math.h.
//
PRMAN_INLINE float
RixMod(float x, float y)
{
    float n = float(int(x/y));
    float result = x - n*y;
    if (result < 0.0f) result += y;
    return result;
}


//
// RixFresnelDielectric
//
// Calculate the reflection coefficient based on the following formula:
//      Kr = (rpar * rpar + rper * rper) / 2.0
//
// where:
//      rcos_i = abs(N.V)  -- with V == -I
//      rcos_t = sqrt( 1.0 - eta * eta * ( 1.0 - rcos_i * rcos_i ) )
//      rpar = ( eta * rcos_t - rcos_i ) / ( eta * rcos_t + rcos_i )
//      rper = ( eta * rcos_i - rcos_t ) / ( eta * rcos_i + rcos_t )
//
// and eta is the relative index of refraction, the ratio of
// refraction indices  ior_i / ior_t
//
// Subscript 'i' denotes the side the ray came from (incident direction,
// ie. the direction V points into); subscript 't' denotes the other
// (refraction/transmission) side.
//
// (The refraction (transmission) coefficient could be computed as:
//      Kt = (1.0 - Kr) * rcos_i / rcos_t
// but usually we just set Kt = 1 - Kr.)
//
// The reflected and transmitted (unit) vectors are optionally calculated
// using the following formula:
//      R = I - 2.0 * (N.I) * N
//        = -V + 2.0 * (N.V) * N
//        = 2.0 * (N.V) * N - V;
//      T = eta * I + (eta * rcos_i - rcos_t) * N
//        = (eta * rcos_i - rcos_t) * N - eta * V;
//
// All vectors have unit length.  Input Vn and results Rn and Tn all point
// away from the surface.  Kr is between 0 and 1.
//
PRMAN_INLINE void
RixFresnelDielectric(float VdN, float eta, float* Kr)
{
    float rcos_i;
    if (VdN < 0.0f)
        rcos_i = -VdN;
    else
        rcos_i = VdN;

    float etaSq    = eta * eta;
    float rcos_tSq = 1.0f - etaSq * (1.0f - rcos_i * rcos_i);
    if (rcos_tSq <= 0.0f)
        *Kr = 1.0f; // total internal reflection
    else
    {
        float rcos_t = std::sqrt(rcos_tSq);
        float rpar   = (eta * rcos_t - rcos_i) / (eta * rcos_t + rcos_i);
        float rper   = (eta * rcos_i - rcos_t) / (eta * rcos_i + rcos_t);
        *Kr          = 0.5f * (rpar * rpar + rper * rper);
    }
}


//
// Computes reflection and refraction directions Rn and Tn.  Both are unit vectors.
//
PRMAN_INLINE void
RixFresnelDielectric(
    const RtVector3& Vn,
    const RtNormal3& Nn,
    float eta,
    float* Kr,
    RtVector3* Rn = NULL,
    RtVector3* Tn = NULL)
{
    float VdN = Dot(Nn, Vn);
    float rcos_i, sign;

    if (VdN < 0.0f)
    {
        rcos_i = -VdN;
        sign   = -1.0f;
    }
    else
    {
        rcos_i = VdN;
        sign   = 1.0f;
    }

    // reflection direction, formula works if Nn and Vn are on
    // opposite or same sides
    if (Rn)
        *Rn = 2.0f * VdN * Nn - Vn;

    float etaSq    = eta * eta;
    float rcos_tSq = 1.0f - etaSq * (1.0f - rcos_i * rcos_i);
    if (rcos_tSq <= 0.0f)
    {
        *Kr = 1.0f; // total internal reflection
        if (Tn)
            *Tn = RtVector3(0.0f); // no refraction (transmission) dir
    }
    else
    {
        float rcos_t = std::sqrt(rcos_tSq);
        float rpar   = (eta * rcos_t - rcos_i) / (eta * rcos_t + rcos_i);
        float rper   = (eta * rcos_i - rcos_t) / (eta * rcos_i + rcos_t);
        *Kr = 0.5f * (rpar * rpar + rper * rper);
        if (Tn)
            *Tn = sign * (eta * rcos_i - rcos_t) * Nn - eta * Vn;
    }
}


//
// RixFresnelDielectric:
// a varying version of RixFresnelDielectric.
//
PRMAN_INLINE void RixFresnelDielectric(
    RixShadingContext& sCtx,
    float eta,
    float* Kr,
    float* Kt, // reflect and refract coeffs
    RtVector3* R,
    RtVector3* T) // reflect and refract dirs
{
    RtNormal3 const* Nn;
    RtVector3 const* Vn;
    sCtx.GetBuiltinVar(RixShadingContext::k_Nn, &Nn);
    sCtx.GetBuiltinVar(RixShadingContext::k_Vn, &Vn);
    int nPts = sCtx.numPts;
    if (R && T)
    {
        for (int i = 0; i < nPts; ++i)
        {
            RixFresnelDielectric(Vn[i], Nn[i], eta, Kr + i, R + i, T + i);
            Kt[i] = 1.0f - Kr[i];
        }
    }
    else
    {
        for (int i = 0; i < nPts; ++i)
        {
            RixFresnelDielectric(Vn[i], Nn[i], eta, Kr + i);
            Kt[i] = 1.0f - Kr[i];
        }
    }
}


//
// RixFresnelConductor:
// calculate the reflection coefficient for a conductor (metal)
// incoming eta, kappa are assumed ior_i / ior_t (for consistency
// with pre-existing dielectric convention). Internally eta and
// kappa are ior_t / ior_i.
//
PRMAN_INLINE void
RixFresnelConductorApprox(float VdN, float eta, float kappa, float* Kr)
{
    if (kappa > 0.01f)
    {
        // approximation for dielectric/conductor interface
        if (eta > 0.0f)
            eta = 1.0f / eta;

        kappa         = 1.0f / kappa;
        float cosSq   = VdN * VdN;
        float cos2eta = 2.0f * eta * VdN;
        float t0      = eta * eta + kappa * kappa;
        float t1      = t0 * cosSq;
        float rperp   = (t0 - cos2eta + cosSq) / (t0 + cos2eta + cosSq);
        float rpar    = (t1 - cos2eta + 1.0f) / (t1 + cos2eta + 1.0f);
        *Kr           = 0.5f * (rpar + rperp);
    }
    else
    {
        // when kappa is small, we're dielectric
        RixFresnelDielectric(VdN, eta, Kr);
    }
}


//
// RixFresnelConductor:
// more expensive computation of fresnel conductor, usually the
// approximation is quite good (except for small kappa).
// See: Moeller's Optics.
//
PRMAN_INLINE void
RixFresnelConductor(float VdN, float eta, float kappa, float* Kr)
{
    if (kappa > .01f)
    {
        float cosSq = VdN * VdN, sinSq = 1.0f - cosSq, sin4 = sinSq * sinSq;

        float t0         = eta * eta - kappa * kappa - sinSq;
        float aSqPlusbSq = std::sqrt(t0 * t0 + 4.0f * kappa * kappa * eta * eta);
        float a          = std::sqrt(0.5f * (aSqPlusbSq + t0));

        float t1 = aSqPlusbSq + cosSq;
        float t2 = 2.0f * a * VdN;

        float rperpSq = (t1 - t2) / (t1 + t2);

        float t3 = aSqPlusbSq * cosSq + sin4;
        float t4 = t2 * sinSq;

        float rparSq = rperpSq * (t3 - t4) / (t3 + t4);

        *Kr = 0.5f * (rperpSq + rparSq);
    }
    else
    {
        // when kappa is small, we're dielectric
        RixFresnelDielectric(VdN, eta, Kr);
    }
}


//
// RixFresnelConductor:
// a version of RixFresnelConductor with the same interface
// as RixFresnelDielectric.
//
PRMAN_INLINE void
RixFresnelConductor(
    const RtVector3& Vn,
    const RtNormal3& Nn,
    float eta,
    float kappa,
    float* Kr,
    RtVector3* Rn = NULL)
{
    float VdN = Dot(Nn, Vn);
    RixFresnelConductorApprox(VdN, eta, kappa, Kr);
    if (Rn)
        *Rn = 2.0f * VdN * Nn - Vn;
}


PRMAN_INLINE void
RixFresnelConductor(
    float VdN,
    RtColorRGB const& eta,
    RtColorRGB const& kappa,
    RtColorRGB* Kr)
{
    RtColorRGB kr;
    RixFresnelConductorApprox(VdN, eta.r, kappa.r, &kr.r);

    // Make sure we compute other bands only when necessary.
    // if eta and kappa are monochromatic, we only compute one band and the cost
    // is the same as a RixFresnelDielectric.
    if (eta.g == eta.r && kappa.g == kappa.r)
        kr.g = kr.r;
    else
        RixFresnelConductorApprox(VdN, eta.g, kappa.g, &kr.g);
    if (eta.b == eta.r && kappa.b == kappa.r)
        kr.b = kr.r;
    else if (eta.b == eta.g && kappa.b == kappa.g)
        kr.b = kr.g;
    else
        RixFresnelConductorApprox(VdN, eta.b, kappa.b, &kr.b);

    *Kr = kr;
}


//
// RixFresnelConductor:
// compute the reflection coefficient for three spectral bands (RGB).
// Similar interface as RixFresnelDielectric.
//
PRMAN_INLINE void
RixFresnelConductor(
    RtVector3 const& Vn,
    RtNormal3 const& Nn,
    RtColorRGB const& eta,
    RtColorRGB const& kappa,
    RtColorRGB* Kr,
    RtVector3* Rn = NULL)
{
    float VdN = Dot(Nn, Vn);
    RixFresnelConductor(VdN, eta, kappa, Kr);
    if (Rn)
        *Rn = 2.0f * VdN * Nn - Vn;
}


//
// Schlick's approximate Fresnel (only the weighting function).
// 0 at normal incidence angle.  1 at grazing angle.
// F(cosTheta, rper) = rper + (1-rper) * RixSchlickFresnelWeight
// where rper = (ior1-ior2)^2 / (ior1+ior2)
//
PRMAN_INLINE float
RixSchlickFresnelWeight(float NdV)
{
    if (NdV <= 0.0f)
        return 1.0f;

    if (NdV >= 1.0f)
        return 0.0f;  // dot prod can be slightly larger than 1.0

    float f  = 1.0f - NdV;
    float f2 = f * f;
    return f2 * f2 * f;  // std::pow(f, 5); -- result between 0.0 and 1.0
}


//
// Reflection direction: standalone computation, useful when Fresnel weights
// aren't needed
//
PRMAN_INLINE RtVector3
RixReflect(const RtVector3& Vn, const RtNormal3& Nn)
{
    float VdN = Dot(Nn, Vn);
    return 2.0f * VdN * Nn - Vn;  // reflection direction
}


//
// Reflection direction: an alternate interface, faster when VdN is already
// available
//
PRMAN_INLINE RtVector3
RixReflect(RtVector3 const& Vn, RtVector3 const& Nn, float VdN)
{
    return 2.0f * VdN * Nn - Vn; // reflection direction
}


//
// Refraction direction: standalone computation using Snell's law.
// All vectors have unit length.  Input Vn and result Tn all point
// away from the surface.  Return value is 1 and Tn = (0,0,0) if total
// internal reflection occured.
//
PRMAN_INLINE int
RixRefract(const RtVector3& Vn, const RtNormal3& Nn, float eta, RtVector3& Tn)
{
    int err = 0;
    float VdN = Dot(Vn, Nn);
    float rcos_i, sign;

    if (VdN < 0.0f)
    {
        rcos_i = -VdN;
        sign   = -1.0f;
    }
    else
    {
        rcos_i = VdN;
        sign   = 1.0f;
    }

    float etaSq    = eta * eta;
    float rcos_tSq = 1.0f - etaSq * (1.0f - rcos_i * rcos_i);
    if (rcos_tSq <= 0.0f)
    {
        Tn  = RtVector3(0.0f);
        err = 1;  // total internal reflection: no refraction
    }
    else
    {
        float rcos_t = std::sqrt(rcos_tSq);
        Tn = sign * (eta * rcos_i - rcos_t) * Nn - eta * Vn;  // refraction dir
        Tn.Normalize();
    }
    return err;
}


PRMAN_INLINE RtNormal3
RixGetForwardFacingNormal(
    RtVector3 const& Vn,
    RtNormal3 const& Nn,
    float* VdotNf = NULL)
{
    RtNormal3 Nf;
    float d = Dot(Vn, Nn);

    if (d > 0.0f)
        Nf = Nn;
    else
    {
        Nf = -Nn;
        d  = -d;
    }

    if (VdotNf)
        *VdotNf = d;

    return Nf;
}


PRMAN_INLINE RtNormal3
RixGetBackwardFacingNormal(
    RtVector3 const& Vn,
    RtNormal3 const& Nn,
    float* VdotNb = NULL)
{
    RtNormal3 Nb;
    float d = Dot(Vn, Nn);

    if (d < 0.0f)
        Nb = Nn;
    else
    {
        Nb = -Nn;
        d  = -d;
    }

    if (VdotNb)
        *VdotNb = d;

    return Nb;
}


#define ALMOSTZERO 0.005f
//
// Given two random numbers, generate a vector from a cosine distribution
// on the upper hemisphere.
// The normal n and two tangent vectors t0 and t1 are assumed to form an
// orthonormal basis: |n| = |t0| = |t1| = 1; n.t1 = n.t2 = t1.t2 = 0.
// The output outDir has unit length.
//
PRMAN_INLINE void
RixCosDirectionalDistribution(
    const RtFloat2& xi,
    const RtVector3& n,
    const RtVector3& t0,
    const RtVector3& t1,
    RtVector3& outDir,
    float& cosTheta)
{
    // In debug mode, assertion will help us identify unexpected inputs
    // where n, t0, t1 do not form an orthonormal basis.
    // Since this func is in the API, any bxdf or plugins that use this func
    // will get the assertion, which probably indicates a bug in the caller.
    assert(n.IsUnitLength());  // close to 1 enough
    assert(t0.IsUnitLength());  // close to 1 enough
    assert(t1.IsUnitLength());  // close to 1 enough
    assert(n.AbsDot(t0) < ALMOSTZERO);  // orthogonal enough
    assert(t0.AbsDot(t1) < ALMOSTZERO);  // orthogonal enough
    assert(n.AbsDot(t1) < ALMOSTZERO);  // orthogonal enough

    float e1 = xi.x * F_TWOPI;
    float z  = std::sqrt(xi.y);
    float r  = std::sqrt(std::max(0.0f, 1.0f - xi.y));  // = sqrt(1 - z^2)
    float x  = r * cosf(e1);
    float y = r * sinf(e1);
    if (z < 1.e-12f)
        z = 1.e-12f;

    outDir   = x * t0 + y * t1 + z * n;
    cosTheta = z;

    assert(outDir.IsUnitLength()); // close to 1 enough
}


//
// Given two random numbers, generate a vector from a cosine distribution
// on the upper hemisphere.  The normal n is assumed to have unit length.
// The output outDir has unit length.
//
PRMAN_INLINE void
RixCosDirectionalDistribution(
    const RtFloat2& xi,
    const RtVector3& n,
    RtVector3& outDir,
    float& cosTheta)
{
    RtVector3 t0, t1;
    n.CreateOrthonormalBasis(t0, t1);
    RixCosDirectionalDistribution(xi, n, t0, t1, outDir, cosTheta);
}


//
// Given two random numbers, generate a vector from a uniform distribution
// on the upper hemisphere.
// The normal n and two tangent vectors t0 and t1 are assumed to form an
// orthonormal basis: |n| = |t0| = |t1| = 1; n.t1 = n.t2 = t1.t2 = 0.
// The output outDir has unit length.
//
PRMAN_INLINE void
RixUniformDirectionalDistribution(
    const RtFloat2& xi,
    const RtVector3& n,
    const RtVector3& t0,
    const RtVector3& t1,
    RtVector3& outDir,
    float& cosTheta)
{
    float e1 = xi.x * F_TWOPI;
    float z  = xi.y;
    float r  = std::sqrt(std::max(0.0f, 1.0f - z * z));
    float x  = r * cosf(e1);
    float y  = r * sinf(e1);
    if (z < 1.e-12f)
        z = 1.e-12f;

    outDir   = x * t0 + y * t1 + z * n;
    cosTheta = z;
}


//
// Given two random numbers, generate a vector from a uniform distribution
// on the upper hemisphere.  The normal n is assumed to have unit length.
// The output outDir has unit length.
//
PRMAN_INLINE void
RixUniformDirectionalDistribution(
    const RtFloat2& xi,
    const RtVector3& n,
    RtVector3& outDir,
    float& cosTheta)
{
    RtVector3 t0, t1;
    n.CreateOrthonormalBasis(t0, t1);
    RixUniformDirectionalDistribution(xi, n, t0, t1, outDir, cosTheta);
}


//
// Given two random numbers and a cone angle, generate a vector from a
// uniform distribution over the cone angle around a normal.
// The normal n and two tangent vectors t0 and t1 are assumed to form an
// orthonormal basis: |n| = |t0| = |t1| = 1; n.t1 = n.t2 = t1.t2 = 0.
// The output outDir has unit length.
//
PRMAN_INLINE void
RixUniformConeDistribution(const RtFloat2& xi, const float coneAngle,
                           const RtVector3& n, const RtVector3& t0,
                           const RtVector3& t1, RtVector3& outDir,
                           float& cosTheta)
{
    float e1 = xi.x * F_TWOPI;
    float cosAngle = cosf(coneAngle);
    float z = xi.y * (1.0f - cosAngle) + cosAngle;
    float r = std::sqrt(std::max(0.0f, 1.0f - z * z));
    float x = r * cosf(e1);
    float y = r * sinf(e1);
    if (z < 1.e-12f)
        z = 1.e-12f;

    outDir   = x * t0 + y * t1 + z * n;
    cosTheta = z;
}


//
// Given two random numbers and a cone angle, generate a vector from a
// uniform distribution over the cone angle around a normal.  The
// normal n is assumed to have unit length.  The output outDir has
// unit length.
//
PRMAN_INLINE void
RixUniformConeDistribution(
    const RtFloat2& xi, const float coneAngle,
    const RtVector3& n,
    RtVector3& outDir,
    float& cosTheta)
{
    RtVector3 t0, t1;
    n.CreateOrthonormalBasis(t0, t1);
    RixUniformConeDistribution(xi, coneAngle, n, t0, t1, outDir, cosTheta);
}


//
// Given two random numbers, generate a vector from a uniform distribution
// on the entire sphere.
//
PRMAN_INLINE RtVector3
RixSphericalDistribution(const RtFloat2& xi)
{
    RtVector3 outDir;
    outDir.z = 2.0f * xi.y - 1.0f;  // cosTheta
    float sinTheta = 1.0f - outDir.z * outDir.z;  // actually sinTheta^2 here

    if (sinTheta > 0.0f)
    {
        sinTheta  = std::sqrt(sinTheta);
        float phi = xi.x * F_TWOPI;
        outDir.x  = sinTheta * cosf(phi);
        outDir.y  = sinTheta * sinf(phi);
    }
    else  // rare case
    {
        outDir.x = 0.0f;
        outDir.y = 0.0f;
    }

    return outDir;
}


//
// Randomly choose between numThresholds+1 scattering lobes using xi,
// and remap xi back to [0,1).
// Making a random selection this way and stretching the random numbers
// preserves stratification properties of a well-distributed set.
//
PRMAN_INLINE int
RixChooseAndRemap(float& xi, int numThresholds, float* thresholds)
{
    // below lowest threshold?
    if (xi < thresholds[0])
    {
        xi /= thresholds[0];
        return 0; // chose lobe 0
    }
    // between thresholds?
    for (int i = 1; i < numThresholds; i++)
    {
        if (thresholds[i - 1] <= xi && xi < thresholds[i])
        {
            xi = (xi - thresholds[i - 1]) / (thresholds[i] - thresholds[i - 1]);
            return i; // chose lobe i
        }
    }
    // must be above highest threshold
    float lastThres = thresholds[numThresholds - 1];
    xi = (xi - lastThres) / (1.0f - lastThres);
    return numThresholds; // chose last lobe
}


//
// Compute ray-tracing bias amount
//
PRMAN_INLINE float
RixTraceBias(
    const RtPoint3& /* org */,
    const RtNormal3& N,
    const RtVector3& dir,
    float biasR,
    float biasT)
{
    return (N.Dot(dir) < 0.0f) ? biasT : biasR;
}


//
// Compute ray-tracing bias amount and apply it to the ray origin
//
PRMAN_INLINE RtPoint3
RixApplyTraceBias(
    const RtPoint3& org,
    const RtNormal3& N,
    const RtVector3& dir,
    const float biasR,
    const float biasT)
{
    float biasAmt = RixTraceBias(org, N, dir, biasR, biasT);

    // Currently we bias in the ray direction, but in some cases we might
    // want to bias along the surface normal (which is also provided as an
    // input to this routine). Note that in the case of a ray in the opposite
    // direction to the normal, we always want to bias in the ray direction
    return org + biasAmt * dir;
}


//
// Compute sin and cos of the angle phi at the same time.
// This is more efficient than calling sin and cos separately.
//
PRMAN_INLINE void
RixSinCos(float phi, float* sinPhi, float* cosPhi)
{
#ifdef LINUX
    sincosf(phi, sinPhi, cosPhi);
#elif (defined(OSX) && defined(__clang__))
    __sincosf(phi, sinPhi, cosPhi);
#else
    *sinPhi = sinf(phi);
    *cosPhi = cosf(phi);
#endif
}


//
// Return a unit vector specified by angles theta and phi.
//
PRMAN_INLINE RtVector3
RixSphericalDirection(
    float sinTheta,
    float cosTheta,
    float sinPhi,
    float cosPhi)
{
    return RtVector3(sinTheta * cosPhi, sinTheta * sinPhi, cosTheta);
}


//
// Return a unit vector specified by angles theta and phi.
//
PRMAN_INLINE RtVector3
RixSphericalDirection(float sinTheta, float cosTheta, float phi)
{
    float sinPhi, cosPhi;
    RixSinCos(phi, &sinPhi, &cosPhi);
    return RixSphericalDirection(sinTheta, cosTheta, sinPhi, cosPhi);
}


//
// RixChangeBasisFrom:
// transform a vector defined relative to the basis vectors
// to the embedding space (local to global).
//
PRMAN_INLINE RtVector3
RixChangeBasisFrom(
    RtVector3 const& in,
    RtVector3 const& X,
    RtVector3 const& Y,
    RtVector3 const& Z)
{
    RtVector3 result;
    result.x = in.x * X.x + in.y * Y.x + in.z * Z.x;
    result.y = in.x * X.y + in.y * Y.y + in.z * Z.y;
    result.z = in.x * X.z + in.y * Y.z + in.z * Z.z;
    return result;
}


PRMAN_INLINE RtVector3
RixChangeBasis(
    RtVector3 const& in,
    RtVector3 const& X,
    RtVector3 const& Y,
    RtVector3 const& Z)
{
    return RixChangeBasisFrom(in, X, Y, Z);
}


//
// RixChangeBasisTo:
// transform a vector defined in the embedding space to be relative
// to the provided basis vectors.  (global to local)
//
PRMAN_INLINE RtVector3
RixChangeBasisTo(
    RtVector3 const& in,
    RtVector3 const& X,
    RtVector3 const& Y,
    RtVector3 const& Z)
{
    return RtVector3(Dot(in, X), Dot(in, Y), Dot(in, Z));
}


//
// Templated mix:  T must define multiply and add (*, + operators)
//
template <typename T>
PRMAN_INLINE T
RixMix(const T& v0, const T& v1, float m)
{
    return (1.0f - m) * v0 + m * v1;
}


//
// Templated smooth (hermite) mix:  T must define *, +, plus scalar +, -)
// This is also known as the "ease" function.
//
template <typename T>
PRMAN_INLINE T
RixSmoothMix(const T& x1, const T& x2, float t)
{
    if (t <= 0.0f)
        return x1;
    else if (t >= 1.0f)
        return x2;
    else
    {
        T x3;
        T dx = x2 - x1;
        x3   = x1 + dx * t * t * (3.0f - 2.0f * t);
        return x3;
    }
}


//
// RixEaseInMix:
// Templated ease-in:  T must define *, +, plus scalar +, -)
// interpolate from x1 to x2, with tangent at x1 of 0
//
template <typename T>
PRMAN_INLINE T
RixEaseInMix(const T& x1, const T& x2, float t)
{
    if (t <= 0.0f)
        return x1;
    else if (t >= 1.0f)
        return x2;
    else
    {
        T x3;
        T dx = x2 - x1;
        x3   = x1 + dx * t * t;
        return x3;
    }
}


//
// RixEaseOutMix:
// Templated ease-in:  T must define *, +, plus scalar +, -)
// interpolate from x1 to x2, with tangent at x2 of 0
//
template <typename T>
PRMAN_INLINE T
RixEaseOutMix(const T& x1, const T& x2, float t)
{
    if (t <= 0.0f)
        return x1;
    else if (t >= 1.0f)
        return x2;
    else
    {
        T x3;
        T dx = x2 - x1;
        x3   = t * dx * (2.0f - t) + x1;
        return x3;
    }
}


// -------------------
// Threshold functions
// -------------------


//
// RixBoxStep:  returns 0 when val < min and 1 otherwise
//
PRMAN_INLINE float
RixBoxStep(float min, float val)
{
    return val < min ? 0.0f : 1.0f;
}


//
// RixLinearStep:  returns 0 when val < min, 1 when val > max and
// a linearly interpolated value between 0 and 1 for val between min and max.
//
PRMAN_INLINE float
RixLinearStep(float min, float max, float val)
{
    if (min < max)
    {
        return val <= min
            ? 0.0f
            : (val >= max ? 1.0f : (val - min) / (max - min));
    }
    else if (min > max)
    {
        return 1.0f - (
            val <= max
                ? 0.0f
                : (val >= min ? 1.0f : (val - max) / (min - max)));
    }
    else
        return RixBoxStep(min, val);
}


//
// RixSmoothStep: returns 0 when val < min, 1 when val > max and
// a smooth Hermite interpolated value between 0 and 1 for val between
// min and max. The slope at both ends of the cubic curve is zero.
//
PRMAN_INLINE float
RixSmoothStep(float min, float max, float val)
{
    if (min < max)
    {
        if (val <= min)
            return 0.0f;
        if (val >= max)
            return 1.0f;
        val = (val - min) / (max - min);
    }
    else if (min > max)
    {
        if (val <= max)
            return 1.0f;
        if (val >= min)
            return 0.0f;
        val = 1.0f - (val - max) / (min - max);
    }
    else
        return RixBoxStep(min, val);

    return val * val * (3.0f - 2.0f * val);
}


//
// RixGaussStep:  returns 0 when val < min, 1 when val > max and
// an exponential curve between 0 and 1 for val between min and max.
// Strictly speaking, this isn't a gaussian curve.
//
PRMAN_INLINE float
RixGaussStep(float min, float max, float val)
{
    if (min < max)
    {
        if (val < min)
            return 0.0f;
        if (val >= max)
            return 1.0f;
        val = 1.0f - (val - min) / (max - min);
    }
    else if (min > max)
    {
        if (val <= max)
            return 1.0f;
        if (val > min)
            return 0.0f;
        val = (val - max) / (min - max);
    }
    else
        return RixBoxStep(min, val);

    return std::pow(2.0f, -8.0f * val * val);
}


//
// Convert a solid angle to a corresponding ray spread
//
PRMAN_INLINE float
RixSolidAngle2Spread(float solidangle)
{
    // To convert from angle to solid angle we use
    // cone-endcap formulation, assuming radius of 1
    //
    // omega = 2.pi.(1-cos(theta))
    //
    // therefore
    //
    // cos(theta) = 1 - omega/(2.pi)
    //
    // we want tan(theta)
    //
    // tan(theta) = + sqrt(1-cos^2(theta)) / cos(theta)
    //  or          - sqrt(1-cos^2(theta)) / cos(theta)

    // Clamp spread to 1.  The solid angle corresponding to spread 1 is
    // 2 pi (1 - cos(pi/4)) ~ 1.840302
    if (solidangle > 1.840302f)
        return 1.0f;

    // Divide solidangle by 2*pi
    float omega_div_2pi = solidangle * 0.15915494f;

    // Compute cos(theta)
    float cosTheta = 1.0f - omega_div_2pi;
    assert(cosTheta > 0.0f);

    // Return tan(theta)
    return std::sqrt(1.0f - cosTheta * cosTheta) / cosTheta;
}


//
// RixAlloca() is a thin veneer atop alloca(). Callers should take
// care to not allocate beyond the platform-specific limits
// on stack size.  To avoid this concern, the ShadingContext
// offers a heap-based memory pool. We implement this as a macro
// instead of an inlined function since it's difficult to enforce
// the requisite inlining in a portable fashion.  The value of
// this veneer is merely that the windows header contortions are
// centralized in this utility header.
//
#define RixAlloca(s) alloca(s)


//
// Inputs Nf and Tn are assumed to be normalised, but are not
// necessarily orthogonal. Parallel Nf and Tn are dealt with
// by ignoring Tn.
//
PRMAN_INLINE void
RixComputeShadingBasis(
    RtNormal3 const& Nf,
    RtVector3 const& Tn,
    RtVector3& TX,
    RtVector3& TY)
{
    // In debug mode, assertion will help us identify unnormalized
    // Nf and Tn inputs.
    // Since this function is in the API, any bxdf or plugins that this
    // will get the assertion, which probably indicates a bug in the caller.
    assert(Nf.IsUnitLength()); // length close to 1
    assert(Tn.IsUnitLength()); // length close to 1

    if (Nf.AbsDot(Tn) < 0.995f)
    {
        TY = Cross(Nf, Tn);
        TY.Normalize();
        TX = Cross(TY, Nf);
    }
    else
    {
        // Under rare circumstances, the input vectors
        // can be parallel. In this case, we generate
        // a basis using only the normal.
        Nf.CreateOrthonormalBasis(TX, TY);
    }
}


PRMAN_INLINE void
RixDebugBasis(
    RixShadingContext* sc,
    RtPoint3& o,
    RtVector3& x,
    RtVector3& y,
    RtVector3& z)
{
    PIXAR_ARGUSED(sc);
    PIXAR_ARGUSED(o);
    PIXAR_ARGUSED(x);
    PIXAR_ARGUSED(y);
    PIXAR_ARGUSED(z);

#ifndef NDEBUG
    float orthonormalCheck = z.Dot(Cross(x, y));
    if (orthonormalCheck < 0.98f)
    {
        const RtColorRGB r = RtColorRGB(1, 0, 0);
        const RtColorRGB g = RtColorRGB(0, 1, 0);
        const RtColorRGB b = RtColorRGB(0, 0, 1);
        RixGeoDebugger* gdb =
            (RixGeoDebugger*)sc->GetRixInterface(k_RixGeoDebugger);
        gdb->EmitPointNormal(o, x, r);
        gdb->EmitPointNormal(o, y, g);
        gdb->EmitPointNormal(o, z, b);
    }
#endif
}


//
// Modify the shading context so that u is shifted
// by one differential.
//
PRMAN_INLINE void
RixShiftCtxInU(RixShadingContext* sCtx)
{
    float const* u;
    float const* du;
    sCtx->GetBuiltinVar(RixShadingContext::k_u, &u);
    sCtx->GetBuiltinVar(RixShadingContext::k_du, &du);

    int numPts = sCtx->numPts;
    RixShadingContext::Allocator pool(sCtx);
    float* newu = (float*)pool.AllocForPattern<float>(numPts);
    for (int i = 0; i < numPts; i++)
    {
        newu[i] = u[i] + du[i];
    }
    sCtx->SetBuiltinVar(RixShadingContext::k_u, newu);
}


//
// Modify the shading context so that v is shifted
// by one differential.
//
PRMAN_INLINE void
RixShiftCtxInV(RixShadingContext* sCtx)
{
    float const* v;
    float const* dv;
    sCtx->GetBuiltinVar(RixShadingContext::k_v, &v);
    sCtx->GetBuiltinVar(RixShadingContext::k_dv, &dv);

    int numPts = sCtx->numPts;
    RixShadingContext::Allocator pool(sCtx);
    float* newv = (float*)pool.AllocForPattern<float>(numPts);
    for (int i = 0; i < numPts; i++)
    {
        newv[i] = v[i] + dv[i];
    }
    sCtx->SetBuiltinVar(RixShadingContext::k_v, newv);
}


//
// Classic Blinn formulation of bumped normal from a height field.
// The three displacements represent the height at the original point,
// at a point du units in the dPdu direction, and at a point dv units
// in the dPdv direction.
//
PRMAN_INLINE void
RixBump(
    RixShadingContext const* sCtx,
    float const* disp,
    float const* dispU,
    float const* dispV,
    RtNormal3* Nn)   // result: bumped normal (normalized)
{
    int numPts = sCtx->numPts;
    float const* du;
    float const* dv;
    RtVector3 const* dPdu;
    RtVector3 const* dPdv;
    RtVector3 const* Norig;
    RtNormal3 const* Ngn;
    RtNormal3 const* Naon;
    sCtx->GetBuiltinVar(RixShadingContext::k_du, &du);
    sCtx->GetBuiltinVar(RixShadingContext::k_dv, &dv);
    sCtx->GetBuiltinVar(RixShadingContext::k_dPdu, &dPdu);
    sCtx->GetBuiltinVar(RixShadingContext::k_dPdv, &dPdv);
    sCtx->GetBuiltinVar(RixShadingContext::k_Nn, &Norig);
    sCtx->GetBuiltinVar(RixShadingContext::k_Ngn, &Ngn);
    sCtx->GetBuiltinVar(RixShadingContext::k_Naon, &Naon);

    // Get primvar N, ie. interpolated vertex normal, if present (not normalized).
    RtNormal3 fill(0,0,0);
    RtNormal3 const *Nuser;
    RixSCDetail nDetail = sCtx->GetPrimVar(Rix::k_N, fill, &Nuser);
    bool hasUserNormals = (nDetail != k_RixSCInvalidDetail);

    for (int i = 0; i < numPts; i++)
    {
        // Only skip computations if all disp values are zero
        if(disp[i] == 0.0f && dispU[i] == 0.0f && dispV[i] == 0.0f)
        {
            Nn[i] = Norig[i];
        }
        else
        {
            // We currently assume the normal does not change
            // significantly over the neighborhood of P + du and P + dv.
            // Should this assumption be invalid a further call to
            // evaluate the normal at each location may be warranted.
            RtPoint3 P0 = Norig[i] * disp[i];
            RtPoint3 P1 = du[i] * dPdu[i] + Norig[i] * dispU[i];
            RtPoint3 P2 = dv[i] * dPdv[i] + Norig[i] * dispV[i];
            Nn[i] = Cross(P1 - P0, P2 - P0);
            Nn[i].Normalize();

            // In some cases, Cross(dPdu, dPdv) is *not* the same as Nn, nor Ngn, mostly because of
            // the camera-to-object transform determinant.
            RtNormal3 computedNgn = Cross(du[i] * dPdu[i], dv[i] * dPdv[i]);
            computedNgn.Normalize();

            float const dotNg = Dot(computedNgn, Ngn[i]);

            if (hasUserNormals)
            {
                // Compensate for faceted renders of coarse polygon meshes (large faces)

                RtNormal3 Nuser_i = Nuser[i];
                RtNormal3 Naon_i = Naon[i];

                // Flip Naon if needed for calculation of delta
                if (dotNg < 0.0f)
                {
                    Naon_i = -Naon_i;
                }

                // Interpolated user-provided vertex normals aren't usually normalized,
                // so normalize now
                Nuser_i.Normalize();

                // Compute difference between user-provided smooth normal and polymesh
                // face normal
                RtNormal3 delta = Nuser_i - Naon_i;

                // Use vertex normals to smooth out faceted renders of coarse polygon meshes
                if (delta != RtNormal3(0.0f))
                {
                   Nn[i] += delta;
                   Nn[i].Normalize();
                }
            }
            else
            {
                // Compensate for faceted renders of coarsely tessellated geometry (large
                // micropolygons)

                // Flip Nn and computedNgn if needed for calculation of delta
                if (dotNg < 0.0f)
                {
                    Nn[i] = -Nn[i];
                    computedNgn = -computedNgn;
                }

                RtNormal3 const delta = Norig[i] - computedNgn;

                // Smooth out coarsely tessellated polygon meshes by adding difference
                // between Nn and Ngn
                if (delta != RtNormal3(0.0f))
                {
                    Nn[i] += delta;
                    Nn[i].Normalize();
                }
            }

        }
    }
}


//
// This routine determines whether the shading points in the provided
// shading context belong to a matte object. Zero is returned if not,
// while a positive value indicates a matte object. A negative value is
// returned on error.
//
PRMAN_INLINE int
RixIsMatte(RixShadingContext const& sCtx)
{
    static const RtUString k_matteName("Ri:Matte");
    static const int k_matteLen    = sizeof(float);
    static const int k_matteDef    = 0;

    int matte = k_matteDef;

    if (RixRenderState* state =
            (RixRenderState*)sCtx.GetRixInterface(k_RixRenderState))
    {
        RixRenderState::Type matteType;
        float matteVal;
        int matteRet;
        int matteCount;

        matteRet = state->GetAttribute(
            k_matteName, &matteVal, k_matteLen, &matteType, &matteCount);

        if (matteRet != 0)
            matte = -1;
        else if (matteCount == 1)
            matte = (int)matteVal;
    }
    else
    {
        matte = -1;
    }

    return matte;
}


PRMAN_INLINE int
RixIsHoldout(RixShadingContext const& sCtx)
{
    static const RtUString k_holdoutName("trace:holdout");
    static const int k_holdoutLen    = sizeof(float);
    static const int k_holdoutDef    = 0;
    int holdout = k_holdoutDef;
    if (RixRenderState* state =
            (RixRenderState*)sCtx.GetRixInterface(k_RixRenderState))
    {
        RixRenderState::Type holdoutType;
        float holdoutVal;
        int holdoutCount;
        int holdoutRet = state->GetAttribute(
            k_holdoutName, &holdoutVal, k_holdoutLen, &holdoutType, &holdoutCount);
        if (holdoutRet)
            holdout = -1;
        else
            holdout = (int)holdoutVal;
    }
    return holdout;
}


PRMAN_INLINE unsigned long
RixImod(const unsigned long a, const unsigned long b)
{
    unsigned long n = a/b;
    return (a - n*b);
}


//
// Compute a hashed float value based on an input string.
// djb2
// http://www.cse.yorku.ca/~oz/hash.html
// This algorithm (k=33) was first reported by dan bernstein many years ago in
// comp.lang.c. Another version of this algorithm (now favored by Bernstein) uses
// xor: hash(i) = hash(i - 1) * 33 ^ str[i]; the magic of number 33 (why it works
// better than many other constants, prime or not) has never been adequately
// explained.
//
PRMAN_INLINE float
RixHash(const char *bufPtr, unsigned long range=65535)
{
    unsigned long hash = 5381;
    int c;
    while ((c = int(*bufPtr++)))
        hash = ((hash << 5) + hash) ^ c;  // hash * 33 + c

    return float(RixImod(hash,range));
}


//
// This function mixes and clamps the color result between a and b to avoid issues
// with floating point precision.
//
PRMAN_INLINE RtColorRGB
RixLerpRGB(RtColorRGB const &a, RtColorRGB const &b, float const mix)
{
    RtColorRGB result = b + mix * (a - b);
    result.r = a.r < b.r ? RixClamp(result.r, a.r, b.r) : RixClamp(result.r, b.r, a.r);
    result.g = a.g < b.g ? RixClamp(result.g, a.g, b.g) : RixClamp(result.g, b.g, a.g);
    result.b = a.b < b.b ? RixClamp(result.b, a.b, b.b) : RixClamp(result.b, b.b, a.b);
    return result;
}


//
// This function adjusts the normals when the geometric normals are facing away
// from the camera.
//
PRMAN_INLINE void
RixAdjustNormal(float const amount, RtVector3 const &Vn, RtVector3 const &Ngn,
                RtNormal3 &Nn)
{
    if (amount == 0.0f)
        return;

    if (Dot(Vn, Ngn) >= 0.0f)
    {
        float VdotN = Dot(Vn, Nn);
        if (VdotN <= 0.0f)
        {
            // If the input normal faces away from camera, nudge it back
            // pad an extra 1% towards V
            Nn -= amount * 1.01f * VdotN * Vn;
            Nn.Normalize();
         }
    }
}

PRMAN_INLINE RtFloat
RixBumpShadowingFactor(RtNormal3 shadingNormal, RtNormal3 bumpNormal, RtVector3 Ln)
{
    // Code from the original paper:
    // "A Microfacet-Based Shadowing Function to Solve the Bump Terminator Problem".
    // Estevez, et al 2019.
    float cos_d = RixMin(std::abs(shadingNormal.Dot(bumpNormal)), 1.0f);
    float tan2_d = (1.0f - cos_d * cos_d) / (cos_d * cos_d);
    float alpha2 = RixClamp(0.125f * tan2_d , 0.0f , 1.0f);

    float cos_i = RixMax(std::abs(shadingNormal.Dot(Ln)), 1e-6f);
    float tan2_i_square = (cos_i * cos_i);
    float tan2_i = (1.0f - cos_i * cos_i) / tan2_i_square;
    if (tan2_i_square < 1e-6f)
    {
        // Division by ZERO undifined;
        tan2_i = 1.0f;
    }

    return (2.0f / (1.0f + sqrtf(1.0f + alpha2 * tan2_i)));
}

// Procedural may have a "user:procprimid" attribute. This is a scene-level
// per-instance hash that may be combined with a procedural-level hash to build
// a globally unique one and create random variations.
PRMAN_INLINE RtFloat
RixProcPrimInstanceHash(RixRenderState *rstate)
{
    float instanceId = 0;
    RixRenderState::Type type;
    int count;
    rstate->GetAttribute(RtUString("user:procprimid"), &instanceId, sizeof(RtFloat), &type, &count);
    return instanceId;
}

//
// Convert ray spread, ray radius, and surface curvature to an equivalent
// (squared) surface roughness.
//
PRMAN_INLINE float
RixRayCurvatureToRoughnessSq(float multiplier, float curvature,
                             float incidentRayRadius, float incidentRaySpread)
{
    float sinTheta = std::max(fabsf(incidentRaySpread), fabsf(incidentRayRadius*curvature));
    sinTheta *= multiplier;

    if (sinTheta == 0.0f) return 0.0f;

    sinTheta = std::min(1.0f, sinTheta);
    float cosTheta = sqrtf(std::max(0.0f, 1.0f-sinTheta*sinTheta));
    // Precision issue, use limit:
    // sin(0.01) ~ 0.0099998 so we could even go a little larger to be safer
    float theta = (sinTheta < 0.01f) ? sinTheta : asinf(sinTheta);

    float sigmaSq = (theta*theta-2.0f) + 2.0f*cosTheta*(theta/sinTheta);
    
    // sigmaSq can be 0.0 or slightly negative due to precision: clamp at 0.0
    return std::max(0.0f, 2.0f*sigmaSq);
}


//
// Combine material roughness with a (squared) roughness computed from the ray
// derivatives and local surface curvature.  Returns a new, larger roughness.
// The colloquial term for this is "roughness mollification"; it is sometimes 
// also called "path regularization".
// (Note that for many bxdfs, the material roughness may already have been squared
// for a more uniform parameter space -- don't let that confuse you.)
//
PRMAN_INLINE float
RixMollifyRoughness(float roughness, float rayCurvatureRoughnessSq)
{
    return sqrtf(roughness*roughness + rayCurvatureRoughnessSq);
}


#endif // RixShadingUtils_h