RixShadingUtils.h

#ifndef RixShadingUtils_h
#define RixShadingUtils_h

/* $Revision: #4 $  $Date: 2016/03/07 $
# ------------------------------------------------------------------------------
#
# Copyright (c) 2012-2014 Pixar Animation Studios. All rights reserved.
#
# The information in this file (the "Software") is provided for the
# exclusive use of the software licensees of Pixar.  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 Animation Studios.
#
# 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.
#
# Pixar
# 1200 Park Ave
# Emeryville CA 94608
#
# ------------------------------------------------------------------------------
*/

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

#include "RixShading.h"

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

#include <cmath>
#include <cassert>
#include <algorithm>

// single-precision floating point constants
#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 )

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);

} // namespace RixConstants

PRMAN_INLINE RtFloat
RixDegreesToRadians(RtFloat degrees)
{
    return degrees * F_DEGTORAD;
}

PRMAN_INLINE RtFloat
RixRadiansToDegress(RtFloat 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
}

PRMAN_INLINE int
RixIsFinite(RtFloat3 const &v)
{
    return RixIsFinite(v.x) && RixIsFinite(v.y) && RixIsFinite(v.z);
}

PRMAN_INLINE RtFloat
RixMin(RtFloat x, RtFloat y)
{
    return x < y ? x : y;
}

PRMAN_INLINE RtFloat
RixMax(RtFloat x, RtFloat y)
{
    return x > y ? x : y;
}

PRMAN_INLINE RtFloat
RixClamp(RtFloat x, RtFloat min, RtFloat max)
{
    return x < min ? min : (x > max) ? max : x;
}

PRMAN_INLINE RtFloat
RixFractional(RtFloat x)
{
    return (x - floorf(x));
}

PRMAN_INLINE void
RixFresnelDielectric(RtFloat VdN, RtFloat eta, RtFloat *Kr)
{
    RtFloat rcos_i;
    if(VdN < 0.f)
        rcos_i = -VdN;
    else
        rcos_i = VdN;

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

PRMAN_INLINE void
RixFresnelDielectric(const RtVector3 &Vn, const RtNormal3 &Nn, RtFloat eta,
                     RtFloat *Kr, RtVector3 *Rn=NULL, RtVector3 *Tn=NULL)
{
    RtFloat VdN = Dot(Nn, Vn);
    RtFloat rcos_i, sign;
    if(VdN < 0.f)
    {
        rcos_i = -VdN;
        sign = -1.f;
    }
    else
    {
        rcos_i = VdN;
        sign = 1.f;
    }

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

    RtFloat etaSq = eta * eta;
    RtFloat 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
    {
        RtFloat rcos_t = std::sqrt(rcos_tSq);
        RtFloat rpar = (eta*rcos_t - rcos_i) / (eta*rcos_t + rcos_i);
        RtFloat 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;
    }
}

PRMAN_INLINE void
RixFresnelDielectric(RixShadingContext &sCtx, RtFloat eta,
                     RtFloat *Kr, RtFloat *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);
    RtInt 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];
        }
    }
}

PRMAN_INLINE void
RixFresnelConductorApprox(RtFloat VdN, RtFloat eta, RtFloat kappa, RtFloat *Kr)
{
    if(kappa > .01f)
    {
        // approximation for dielectric/conductor interface
        if(eta > 0.f)
            eta = 1.f / eta;
        kappa = 1.f / kappa;
        float cosSq = VdN*VdN;
        float cos2eta = 2.f * eta * VdN;
        RtFloat t0 = eta * eta + kappa * kappa;
        RtFloat t1 = t0 * cosSq;
        RtFloat rperp = (t0 - cos2eta + cosSq) / (t0 + cos2eta + cosSq);
        RtFloat rpar = (t1 - cos2eta + 1.f) / (t1 + cos2eta + 1.f);
        *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(RtFloat VdN, RtFloat eta, RtFloat kappa, RtFloat *Kr)
{
    if(kappa > .01f)
    {
        float cosSq = VdN * VdN,
              sinSq = 1.f-cosSq,
              sin4 = sinSq*sinSq;

        float t0 = eta*eta - kappa*kappa - sinSq;
        float aSqPlusbSq = std::sqrt(t0*t0 + 4.f*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);
    }
}

PRMAN_INLINE void
RixFresnelConductor(const RtVector3 &Vn, const RtNormal3 &Nn,
                    RtFloat eta, RtFloat kappa,
                    RtFloat *Kr, RtVector3 *Rn=NULL)
{
    RtFloat VdN = Dot(Nn, Vn);
    RixFresnelConductorApprox(VdN, eta, kappa, Kr);
    if(Rn) *Rn = 2.0f * VdN * Nn - Vn;
}

PRMAN_INLINE void
RixFresnelConductor(RtFloat 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;
}

PRMAN_INLINE void
RixFresnelConductor(RtVector3 const &Vn, RtNormal3 const &Nn,
                    RtColorRGB const &eta, RtColorRGB const &kappa,
                    RtColorRGB *Kr, RtVector3 *Rn=NULL)
{
    RtFloat VdN = Dot(Nn, Vn);
    RixFresnelConductor(VdN, eta, kappa, Kr);
    if(Rn) *Rn = 2.0f * VdN * Nn - Vn;
}

PRMAN_INLINE RtFloat
RixSchlickFresnelWeight(RtFloat NdV)
{
    if (NdV <= 0.f) return 1.f;
    if (NdV >= 1.f) return 0.f; // 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
}

PRMAN_INLINE RtVector3
RixReflect(const RtVector3 &Vn, const RtNormal3 &Nn)
{
    RtFloat VdN = Dot(Nn, Vn);
    return  2.0f * VdN * Nn - Vn; // reflection direction
}

PRMAN_INLINE RtVector3
RixReflect(RtVector3 const &Vn, RtVector3 const &Nn, RtFloat VdN)
{
    return 2.0f * VdN * Nn - Vn; // reflection direction
}

PRMAN_INLINE int
RixRefract(const RtVector3 &Vn, const RtNormal3 &Nn, RtFloat eta, RtVector3 &Tn)
{
    int err = 0;
    RtFloat VdN = Dot(Vn, Nn);
    RtFloat rcos_i, sign;
    if (VdN < 0.f)
    {
        rcos_i = -VdN;
        sign = -1.f;
    }
    else
    {
        rcos_i = VdN;
        sign = 1.f;
    }

    RtFloat etaSq = eta*eta;
    RtFloat 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
    {
        RtFloat 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,
                          RtFloat *VdotNf=NULL)
{
    RtNormal3 Nf;
    float d = Dot(Vn, Nn);
    if(d > 0.f)
        Nf = Nn;
    else
    {
        Nf = -Nn;
        d = -d;
    }
    if(VdotNf)
        *VdotNf = d;
    return Nf;
}

PRMAN_INLINE RtNormal3
RixGetBackwardFacingNormal(RtVector3 const &Vn, RtNormal3 const &Nn,
                          RtFloat *VdotNb=NULL)
{
    RtNormal3 Nb;
    float d = Dot(Vn, Nn);
    if(d < 0.f)
        Nb = Nn;
    else
    {
        Nb = -Nn;
        d = -d;
    }
    if(VdotNb)
        *VdotNb = d;
    return Nb;
}

PRMAN_INLINE void
RixCosDirectionalDistribution(const RtFloat2 &xi, const RtVector3 &n,
                              const RtVector3 &t0, const RtVector3 &t1,
                              RtVector3 &outDir, float &cosTheta)
{
    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;
}

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);
}

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;
}

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);
}

PRMAN_INLINE RtVector3
RixSphericalDistribution(const RtFloat2 &xi)
{
    RtVector3 outDir;
    outDir.z = xi.y * 2.0f - 1.0f; // cosTheta
    float sinTheta = 1.0f - outDir.z * outDir.z;
    if (sinTheta <= 0.0f)
    {
        outDir.x = 0.0f;
        outDir.y = 0.0f;
    }
    else
    {
        sinTheta = std::sqrt(sinTheta);
        float phi = xi.x * F_TWOPI;
        outDir.x = sinTheta * cosf(phi);
        outDir.y = sinTheta * sinf(phi);
    }
    return outDir;
}

PRMAN_INLINE int
RixChooseAndRemap(RtFloat &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 RtFloat
RixTraceBias(const RtPoint3 & /* org */, const RtNormal3 &N, const RtVector3 &dir,
             RtFloat biasR, RtFloat 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 RtFloat biasR,
                  const RtFloat biasT)
{
    RtFloat biasAmt = RixTraceBias(org, N, dir, biasR, biasT);

    return org + biasAmt * dir;
}


PRMAN_INLINE void
RixSinCos(RtFloat phi, RtFloat *sinphi, RtFloat *cosphi)
{
#ifdef LINUX
    sincosf(phi, sinphi, cosphi);
#else
    *sinphi = sinf(phi);
    *cosphi = cosf(phi);
#endif
}

PRMAN_INLINE RtVector3
RixSphericalDirection(RtFloat sintheta, RtFloat costheta,
                   RtFloat sinphi, RtFloat cosphi)
{
   return RtVector3(sintheta*cosphi, sintheta*sinphi, costheta);
}

PRMAN_INLINE RtVector3
RixSphericalDirection(RtFloat sintheta, RtFloat costheta, RtFloat phi)
{
    RtFloat sinphi, cosphi;
    RixSinCos(phi, &sinphi, &cosphi);
    return RixSphericalDirection(sintheta, costheta, sinphi, cosphi);
}


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);
}

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, RtFloat 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, RtFloat 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.f - 2.f*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, RtFloat 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, RtFloat 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.f - t) + x1;
        return x3;
    }
}

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

// RixBoxStep:  returns 0 when val < min and 1 otherwise
PRMAN_INLINE RtFloat
RixBoxStep(RtFloat min, RtFloat 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 RtFloat
RixLinearStep(RtFloat min, RtFloat max, RtFloat 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 RtFloat
RixSmoothStep(RtFloat min, RtFloat max, RtFloat val)
{
    if (min < max)
    {
        if (val <= min)  return 0.f;
        if (val >= max) return 1.f;
        val = (val - min)/(max - min);
    }
    else 
    if (min > max)
    {
        if (val <= max) return 1.f;
        if (val >= min)  return 0.f;
        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 RtFloat
RixGaussStep(RtFloat min, RtFloat max, RtFloat 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);
}

PRMAN_INLINE RtFloat
RixSolidAngle2Spread(RtFloat solidangle)
{
    // 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
    RtFloat omega_div_2pi = solidangle * 0.15915494f;

    // Compute cos(theta)
    RtFloat costheta = 1.0f - omega_div_2pi;
    assert(costheta > 0.0f);

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


#define RixAlloca(s)  alloca(s)

// Inputs Nf and Tn are assumed to be normalised, but are not
// necessarily orthogonal. They must not be parallel however.
PRMAN_INLINE void
RixComputeShadingBasis(RtNormal3 const &Nf, RtVector3 const &Tn,
                      RtVector3 &TX, RtVector3 &TY)
{
    TY = Cross(Nf, Tn);
    TY.Normalize();
    TX = Cross(TY, Nf);
}

PRMAN_INLINE void
RixDebugBasis(RixShadingContext *sc, RtPoint3 &o,
              RtVector3 &x, RtVector3 &y, RtVector3 &z)
{
#ifndef NDEBUG
    RtFloat 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)
{
    RtFloat const *u;
    RtFloat const *du;
    sCtx->GetBuiltinVar(RixShadingContext::k_u, &u);
    sCtx->GetBuiltinVar(RixShadingContext::k_du, &du);

    int numPts = sCtx->numPts;
    RixShadingContext::Allocator pool(sCtx);
    RtFloat* newu = (RtFloat*)pool.AllocForPattern<RtFloat>(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)
{
    RtFloat const *v;
    RtFloat const *dv;
    sCtx->GetBuiltinVar(RixShadingContext::k_v, &v);
    sCtx->GetBuiltinVar(RixShadingContext::k_dv, &dv);

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

PRMAN_INLINE void
RixBump(RixShadingContext const *sCtx, RtFloat const* disp,
        RtFloat const* dispU, RtFloat const* dispV,
        RtNormal3* Nn)
{
    int numPts = sCtx->numPts;
    RtFloat const *du;
    RtFloat const *dv;
    RtVector3 const *dPdu;
    RtVector3 const *dPdv;
    RtVector3 const *Norig;
    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);
    for(int i = 0; i < numPts; i++)
    {
        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();
    }
}

PRMAN_INLINE RtInt
RixIsMatte(RixShadingContext const &sCtx)
{
    static const char    *k_matteName = "Matte";
    static const int      k_matteLen  = sizeof(RtFloat);
    static const RtInt    k_matteDef  = 0;
    
    RtInt matte = k_matteDef;

    if (RixRenderState *state = 
        (RixRenderState *) sCtx.GetRixInterface(k_RixRenderState))
    {
        RixRenderState::Type matteType;
        RtFloat matteVal;
        RtInt matteRet;
        RtInt matteCount;
        
        matteRet = state->GetAttribute(k_matteName, &matteVal, 
                                       k_matteLen,  &matteType, &matteCount);
   
        if (matteRet != 0)
            matte = -1;
        else if (matteCount == 1)
            matte = (RtInt) matteVal;
    }
    else
    {
        matte = -1;
    }

    return matte;
}

#endif // RixShadingUtils_h