RixIESInline.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 RixIESInline_h
#define RixIESInline_h

#include <math.h>                       // for floorf, powf, asinf, atan2f, etc
#include <algorithm>                    // for max, min
#include <cfloat>                       // for FLT_MIN
#include <string>                       // for string
#include "RixInterfaces.h"              // for RixTexture, etc
#include "RixShadingUtils.h"            // for F_PI
#include "prmanapi.h"                   // for PRMAN_INLINE
#include "RiTypesHelper.h"              // for RtColorRGB, RtFloat3, etc

PRMAN_INLINE void
RixIES::EvaluateIESProfile(RtFloat3 const& P, RtColorRGB* emission) const
{
    // if we haven't read a map...
    if (!m_iesProfile.profile)
    {
        emission->r = emission->g = emission->b = 1.f;
        return;
    }

    float vScale      = 0.5f;
    RtFloat3 localDir = P;

    // assume projection center is at shader space 0
    // Transform(P, &localDir, m_profileSpace, 1);

    if (!m_hasNegativeDeterminant)
        localDir *= -1.0f;

    float u;
    float v;
    calcLatLongParam(1, (float*)&localDir, m_iesProfile.isExr, &u, &v);

    if ((u != u) || (v != v))
    {
        // FIXME ideally we would want to report an error as opposed to silently
        // return 1
        emission->r = emission->g = emission->b = 1.f;
        return;
    }

    if (u < 0.0f || u > 1.0f || v < 0.0f || v > 1.0f)
    {
        // FIXME ideally we would want to report an error as opposed to silently
        // return 1
        emission->r = emission->g = emission->b = 1.f;
        return;
    }

    // rescale
    float coneSize     = 1.0f;
    float profileScale = 1.0f;
    if (m_coneAngle > 0.0f)
    {
        coneSize     = m_coneAngle / (float)F_PI;
        profileScale = coneSize;
    }
    if (m_iesProfile.profileScale > 0.0f)
        profileScale = m_iesProfile.profileScale;

    // on axis = 0, opposite = 1

    // map the cone to 0:1 on the map
    v = v * coneSize; // allow overflow to be dealt with by profileScale

    // profileScale 180 (1.0) = full map < 180 restricts the region

    v = v / profileScale;

    // orient against map
    v = 1.0f - v;

    v = std::max(std::min(v * vScale, 1.0f), 0.0f);

    // now exactly on axis == 1, opposite = 0

    float uf = u * static_cast<float>(m_iesProfile.nu);
    if (uf >= m_iesProfile.nu)
        uf -= static_cast<float>(m_iesProfile.nu);

    int u0 = static_cast<int>(floorf(uf));
    int u1 = u0 + 1;
    if (u1 >= m_iesProfile.nu)
        u1 = 0; // wrap

    float fu = uf - static_cast<float>(u0);

    if (u0 < 0 || u0 >= m_iesProfile.nu || u1 < 0 || u1 >= m_iesProfile.nu)
    {
        // FIXME ideally we would want to report an error as opposed to silently
        // return 1
        emission->r = emission->g = emission->b = 1.f;
        return;
    }

    float vf = v * (static_cast<float>(m_iesProfile.nv) - 1.0f);
    if (vf >= m_iesProfile.nv)
        vf = static_cast<float>(m_iesProfile.nv) - 1.0f;

    int v0 = static_cast<int>(floorf(vf));
    int v1 = v0 + 1;
    if (v1 >= m_iesProfile.nv)
        v1 = m_iesProfile.nv - 1;

    float fv = vf - static_cast<float>(v0);

    if (v0 < 0 || v0 >= m_iesProfile.nv || v1 < 0 || v1 >= m_iesProfile.nv)
    {
        // FIXME ideally we would want to report an error as opposed to silently
        // return 1
        emission->r = emission->g = emission->b = 1.f;
        return;
    }

    int v00 = u0 + v0 * m_iesProfile.nu;
    int v10 = u1 + v0 * m_iesProfile.nu;
    int v01 = u0 + v1 * m_iesProfile.nu;
    int v11 = u1 + v1 * m_iesProfile.nu;

    float rfu = 1.0f - fu;

    // clang-format off
    // Clamp to 0 as some intepolated maps have -ve emission!!
    if (m_iesProfile.isColor)
    {
        emission->r =
            (1.0f - fv) * ( rfu * m_iesProfile.profile[3 * v00 + 0] 
                           + fu * m_iesProfile.profile[3 * v10 + 0])
            +       fv  * ( rfu * m_iesProfile.profile[3 * v01 + 0] 
                           + fu * m_iesProfile.profile[3 * v11 + 0]);
        emission->r = std::max(0.0f, emission->r);

        emission->g =
            (1.0f - fv) * ( rfu * m_iesProfile.profile[3 * v00 + 1] 
                           + fu * m_iesProfile.profile[3 * v10 + 1])
            +       fv  * ( rfu * m_iesProfile.profile[3 * v01 + 1] 
                           + fu * m_iesProfile.profile[3 * v11 + 1]);
        emission->g = std::max(0.0f, emission->g);

        emission->b =
            (1.0f - fv) * ( rfu * m_iesProfile.profile[3 * v00 + 2] 
                           + fu * m_iesProfile.profile[3 * v10 + 2])
            +       fv  * ( rfu * m_iesProfile.profile[3 * v01 + 2] 
                           + fu * m_iesProfile.profile[3 * v11 + 2]);
        emission->b = std::max(0.0f, emission->b);
    }
    else
    {
        emission->r = 
            (1.0f - fv) * ( rfu * m_iesProfile.profile[v00] 
                           + fu * m_iesProfile.profile[v10])
            +       fv  * ( rfu * m_iesProfile.profile[v01] 
                           + fu * m_iesProfile.profile[v11]);

        emission->r = std::max(0.0f, emission->r);
        emission->g = emission->b = emission->r;
    }
    // clang-format on
}

PRMAN_INLINE bool
RixIES::ReadIESProfile(RixContext const& rixCtx, RtUString const iesProfile)
{
    if (m_iesProfile.profile)
        m_iesProfile.clear(); // clear existing if we already read data before

    // AcquireTexture will crash with an empty string...
    if (!iesProfile.Empty())
    {
        static std::string rtxPrefix("rtxplugin:RtxIES?filename=");

        // FIXME: this doesn't work on Windows...
        // m_iesProfile.mapName = rtxPrefix + std::string(iesProfile);

        // FIXME: The unique string construction should be made constant over
        // the render for performance reasons.
        std::string mapName = rtxPrefix + std::string(iesProfile.CStr());
        m_iesProfile.mapName = RtUString(mapName.c_str());

        RixMessages* msgs = (RixMessages*)rixCtx.GetRixInterface(k_RixMessages);
        RixTexture* tex   = (RixTexture*)rixCtx.GetRixInterface(k_RixTexture);
        // open the texture
        RixTexture::TxProperties txProps;
        if (tex->AcquireTexture(m_iesProfile.mapName, RixTexture::AtlasNone,
                                txProps))
        {
            msgs->Error(
                "RixIES emission profile: could not open %s",
                m_iesProfile.mapName.CStr());
            m_iesProfile.mapName = US_NULL;
            return false;
        }
        else
        {
            msgs->Info(
                "Building RixIES emission profile(\"%s\")\n",
                m_iesProfile.mapName.CStr());

            m_iesProfile.nu      = txProps.width;
            m_iesProfile.nv      = txProps.height;
            m_iesProfile.isExr   = false;
            m_iesProfile.isColor = (txProps.nchannels == 3);

            // read data buffer
            m_iesProfile.profile = new float[
                m_iesProfile.nu * m_iesProfile.nv * txProps.nchannels];

            RixTexture::TxParams txParams;
            txParams.nchannels = txProps.nchannels;
            // txParams.firstchannel = 0; // start at the first channel

            tex->TextureData(txProps, txParams, 0, 0, m_iesProfile.profile);
            tex->ReleaseTexture(txProps);

            int nchans = txProps.nchannels;
            int nu = m_iesProfile.nu;
            int nv = m_iesProfile.nv;
            float total = 0.f;
            // nv/2 is the extent of the vertical range of the profile, but
            // because profiles specify data at both the lower boundary and
            // upper boundary we need to include both the nv = 0 and nv/2 rows.
            for (int j = 0; j <= nv / 2; ++j)
            {
                // The area of a spherical cap is 2*pi*r^2*(1-cos(theta)), where
                // theta is the polar angle.  To get the area of a spherical
                // strip between two polar angles, we subtract the area of the
                // smaller cap from the area of the larger one.  Our scale
                // factor is then the ratio of the strip area to the total
                // surface area of the sphere, so the 2*pi*r^2 factor divided by
                // 4*pi*r^2 turns into just 0.5.
                // At the endpoints, we clamp the angle to 0 or pi (since the
                // strip can't extend past the poles of the sphere) and the
                // strip ends up being just a cap for those cases.
                float lowerAngle = (j == 0) ? 0 :  F_PI * (2 * j - 1) / nv;
                float upperAngle =
                    (j == nv / 2) ? F_PI : F_PI * (2 * j + 1) / nv;
                float scalefactor =
                    0.5f * (cosf(lowerAngle) - cosf(upperAngle));
                for (int i = 0; i < nu; ++i)
                {
                    for (int k = 0; k < nchans; ++k)
                    {
                        total += scalefactor * 
                            m_iesProfile.profile[((j * nu) + i) * nchans + k];
                    }
                }
            }
            m_normalizationFactor = nu * nchans / total;
            return true;
        }
    }
    return false;
}

PRMAN_INLINE float
RixIES::EvaluateConeAngle(const RtFloat3& P) const
{
    float res = 1.0f;
    if (m_coneAngle > 0.0f)
    {
        res = 0.0f;
        RtFloat3 localDir = P;
        if (m_hasNegativeDeterminant)
            localDir *= -1.0f;

        if (localDir.z > 0.0f)
        {
            localDir.Normalize();
            if (localDir.z >= m_cosConeAngle)
            {
                res = 1.0f;
            }
        }
    }
    return res;
}

PRMAN_INLINE float
RixIES::EvaluateConeAngle(RtPoint3 const& lightP, RtVector3& lightN,
                          RtPoint3 const& offsetP, RtVector3 const& segmentDir,
                          float& minT, float& maxT) const
{
    // TODO - this line segment vs cone algorithm comes from David Eberly. We either need
    // to give appropriate credit, or find an alternative method.
    float ts = lightN.Dot(offsetP + minT * segmentDir - lightP);
    float te = lightN.Dot(offsetP + maxT * segmentDir - lightP);
    if (ts < 0.0f && te < 0.0f) return 0.0f;
    if (ts * te < 0.0f)
    {
        float planeT = lightN.Dot(lightP - offsetP) / segmentDir.Dot(lightN);
        if (ts < 0.0f)
        {
            minT = planeT; 
        }
        else
        {
            maxT = planeT;
        }            
    }
    if (m_cosConeAngle < 1e-5f) return 1.0f;

    float vDn = segmentDir.Dot(lightN);
    RtVector3 lToP = offsetP - lightP;
    float c2 = vDn * vDn - segmentDir.Dot(segmentDir) * m_cosConeAngle * m_cosConeAngle;
    if (c2 == 0.0f) return 0.0f;
    float lpDn = lToP.Dot(lightN);
    float c1 = vDn * lpDn - lToP.Dot(segmentDir) * m_cosConeAngle * m_cosConeAngle;
    float c0 = lpDn * lpDn - lToP.Dot(lToP) * m_cosConeAngle * m_cosConeAngle;

    float delta = c1*c1 - c0*c2;

    if (delta > 0.0f)
    {
        float rootDelta = sqrtf(delta);
        float t0;
        float t1;
        if (c2 > 0.0f)
        {
            t0 = (-c1 - rootDelta) / c2;
            t1 = (-c1 + rootDelta) / c2;
        }
        else
        {
            t1 = (-c1 - rootDelta) / c2;
            t0 = (-c1 + rootDelta) / c2;
        }
        assert( t1 >= t0 );

        minT = std::max(minT, t0);
        maxT = std::min(maxT, t1);
        if (minT >= maxT)
        {
            return 0.0f;
        }
        return 1.0f;
    }
    else
    {
        return 0.0f;
    }

}

// Utility routines
PRMAN_INLINE float
RixIES::calcLongitude(float x, float y, bool isExr) const
{
    // Compensate for the fact that OpenEXR textures store the
    // data with a latitude range of [180, -180] from left to
    // right, while other formats store the data with a latitude
    // range of [0, 360].
    float s = atan2f(y, x);
    if (isExr)
    {
        s *= -0.5f / F_PI;
        if (s < 0.0f)
            s += 1.0f;
    }
    else
    {
        s *= 0.5f / F_PI;
        s += 0.5f;
    }
    return s;
}

PRMAN_INLINE void RixIES::calcLatLongParam(
    int n,           // number of points in BOP
    float* dir,      // array of reflection directions
    bool openEXR,    // Is this an EXR latlong?
    float* s,       // output s array pointer
    float* t) const // output t array pointer
{
    // this func calculates st values from direction for a lat-long texture map

    float x, y, z, len, zlen;
    for (int i = 0; i < n; i++, s++, t++, dir += 3)
    {
        /* compute texture coordinates from direction */
        /* normalize */
        x = dir[0];
        y = dir[1];
        z = dir[2];

        len = x * x + y * y + z * z;
        if (len > FLT_MIN)
            len = 1.f / sqrtf(len);

        zlen = z * len;

        /*  This clamps the value passed to asinf to avoid
        *   precision issues that might occur in the above
        *   computation. */
        if (zlen < -1.0f)
            zlen = -1.0f;
        if (zlen > 1.0f)
            zlen = 1.0f;

        *t = 0.5f - asinf(zlen) * 1.0f / F_PI;
        *s = calcLongitude(x, y, openEXR);
    }
}

#endif