triaccel.h

Go to the documentation of this file.

/*
    This file is part of Mitsuba, a physically based rendering system.

    Copyright (c) 2007-2014 by Wenzel Jakob and others.

    Mitsuba is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License Version 3
    as published by the Free Software Foundation.

    Mitsuba is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

#pragma once
#if !defined(__MITSUBA_RENDER_TRIACCEL_H_)
#define __MITSUBA_RENDER_TRIACCEL_H_

#include <mitsuba/render/trimesh.h>

MTS_NAMESPACE_BEGIN

// Used when a fake triangle is used to reference a non-triangle shape instance
static const uint32_t KNoTriangleFlag = 0xFFFFFFFF;

/**
 * \brief Pre-computed triangle representation based on Ingo Wald's TriAccel layout.
 *
 * Fits into three 16-byte cache lines if single precision floats are used.
 * The k parameter is also used for classification during kd-tree construction.
 * \ingroup librender
 */
struct TriAccel {
        uint32_t k;
        Float n_u;
        Float n_v;
        Float n_d;

        Float a_u;
        Float a_v;
        Float b_nu;
        Float b_nv;

        Float c_nu;
        Float c_nv;
        uint32_t shapeIndex;
        uint32_t primIndex;

        /// Construct from vertex data. Returns '1' if there was a failure
        inline int load(const Point &A, const Point &B, const Point &C);

        /// Fast ray-triangle intersection test
        FINLINE bool rayIntersect(const Ray &ray, Float mint, Float maxt,
                Float &u, Float &v, Float &t) const;
};

inline int TriAccel::load(const Point &A, const Point &B, const Point &C) {
        static const int waldModulo[4] = { 1, 2, 0, 1 };

        Vector b = C-A, c = B-A, N = cross(c, b);

        k = 0;
        /* Determine the largest projection axis */
        for (int j=0; j<3; j++) {
                if (std::abs(N[j]) > std::abs(N[k]))
                        k = j;
        }

        uint32_t u = waldModulo[k],
                v = waldModulo[k+1];
        const Float n_k = N[k],
                denom = b[u]*c[v] - b[v]*c[u];

        if (denom == 0) {
                k = 3;
                return 1;
        }

        /* Pre-compute intersection calculation constants */
        n_u   =  N[u] / n_k;
        n_v   =  N[v] / n_k;
        n_d   =  dot(Vector(A), N) / n_k;
        b_nu  =  b[u] / denom;
        b_nv  = -b[v] / denom;
        a_u   =  A[u];
        a_v   =  A[v];
        c_nu  =  c[v] / denom;
        c_nv  = -c[u] / denom;
        return 0;
}

FINLINE bool TriAccel::rayIntersect(const Ray &ray, Float mint, Float maxt,
        Float &u, Float &v, Float &t) const {

#if 0
        static const MM_ALIGN16 int waldModulo[4] = { 1, 2, 0, 1 };
        const int ku = waldModulo[k], kv = waldModulo[k+1];
        /* Get the u and v components */
        const Float o_u = ray.o[ku], o_v = ray.o[kv], o_k = ray.o[k],
                                d_u = ray.d[ku], d_v = ray.d[kv], d_k = ray.d[k];
#else
        Float o_u, o_v, o_k, d_u, d_v, d_k;
        switch (k) {
                case 0:
                        o_u = ray.o[1];
                        o_v = ray.o[2];
                        o_k = ray.o[0];
                        d_u = ray.d[1];
                        d_v = ray.d[2];
                        d_k = ray.d[0];
                        break;
                case 1:
                        o_u = ray.o[2];
                        o_v = ray.o[0];
                        o_k = ray.o[1];
                        d_u = ray.d[2];
                        d_v = ray.d[0];
                        d_k = ray.d[1];
                        break;
                case 2:
                        o_u = ray.o[0];
                        o_v = ray.o[1];
                        o_k = ray.o[2];
                        d_u = ray.d[0];
                        d_v = ray.d[1];
                        d_k = ray.d[2];
                        break;
                default:
                        return false;
        }
#endif


#if defined(MTS_DEBUG_FP)
        if (d_u * n_u + d_v * n_v + d_k == 0)
                return false;
#endif

        /* Calculate the plane intersection (Typo in the thesis?) */
        t = (n_d - o_u*n_u - o_v*n_v - o_k) /
                (d_u * n_u + d_v * n_v + d_k);

        if (t < mint || t > maxt)
                return false;

        /* Calculate the projected plane intersection point */
        const Float hu = o_u + t * d_u - a_u;
        const Float hv = o_v + t * d_v - a_v;

        /* In barycentric coordinates */
        u = hv * b_nu + hu * b_nv;
        v = hu * c_nu + hv * c_nv;
        return u >= 0 && v >= 0 && u+v <= 1.0f;
}

MTS_NAMESPACE_END

#endif /* __MITSUBA_RENDER_TRIACCEL_H_ */