Glint BRDF: Discrete Microfacet Implementation for Nori

glints.cpp

#include <nori/bsdf.h>
#include <nori/frame.h>
#include <nori/warp.h>
#include <nori/quadtree.h>

using namespace std;

NORI_NAMESPACE_BEGIN

    class Glints : public BSDF {
    public:
        Glints(const PropertyList &propList) {
            /* RMS surface roughness */
            m_alpha = propList.getFloat("alpha", 0.1f);

            /* Interior IOR (default: BK7 borosilicate optical glass) */
            m_intIOR = propList.getFloat("intIOR", 1.5046f);

            /* Exterior IOR (default: air) */
            m_extIOR = propList.getFloat("extIOR", 1.000277f);

            /* Albedo of the diffuse base material (a.k.a "kd") */
            m_kd = propList.getColor("kd", Color3f(0.5f));

            /* To ensure energy conservation, we must scale the
               specular component by 1-kd.

               While that is not a particularly realistic model of what
               happens in reality, this will greatly simplify the
               implementation. Please see the course staff if you're
               interested in implementing a more realistic version
               of this BRDF. */
            m_ks = 1 - m_kd.maxCoeff();
        }

        float calculateG(Vector3f wv, Vector3f wh) const {
            // since the vectors are normalised, wv dot n is simply cos theta
            // where theta is the angle between wv and n
            float c = wv.dot(wh)/Frame::cosTheta(wv);
            if (c > 0) {
                float b = 1.0f/(m_alpha*Frame::tanTheta(wv));
                if (b < 1.6) {
                    return (3.535*b + 2.181*b*b)/(1 + 2.276*b + 2.577*b*b);
                }
                else {
                    return 1.0f;
                }
            }

            return .0f;

        }

        /// Evaluate the BRDF for the given pair of directions
        Color3f evalGlint(const BSDFQueryRecord &bRec, const Point2f footprint[4]) const {
            if (bRec.measure != ESolidAngle || Frame::cosTheta(bRec.wi) <= 0 || Frame::cosTheta(bRec.wo) <= 0) {
                return Color3f(.0f);
            }

            Vector3f wh = (bRec.wi + bRec.wo).normalized();

            float D = computeDiscreteD(bRec, footprint);
            float F = fresnel(wh.dot(bRec.wi), m_extIOR, m_intIOR);
            float G = calculateG(bRec.wi, wh) * calculateG(bRec.wo, wh);

            return D * F * G / (4 * Frame::cosTheta(bRec.wi) * Frame::cosTheta(bRec.wo) * Frame::cosTheta(wh));
        }

        /// Evaluate the BRDF for the given pair of directions
        Color3f eval(const BSDFQueryRecord &bRec) const {
            if (bRec.measure != ESolidAngle || Frame::cosTheta(bRec.wi) <= 0 || Frame::cosTheta(bRec.wo) <= 0) {
                return Color3f(.0f);
            }

            Vector3f wh = (bRec.wi + bRec.wo).normalized();

            float D = Warp::squareToBeckmannPdf(wh, m_alpha);
            float F = fresnel(wh.dot(bRec.wi), m_extIOR, m_intIOR);
            float G = calculateG(bRec.wi, wh) * calculateG(bRec.wo, wh);

            return m_kd/M_PI + m_ks * D * F * G / (4 * Frame::cosTheta(bRec.wi) * Frame::cosTheta(bRec.wo) * Frame::cosTheta(wh));
        }

        /// Evaluate the sampling density of \ref sample() wrt. solid angles
        float pdf(const BSDFQueryRecord &bRec) const {
            if (bRec.measure != ESolidAngle || Frame::cosTheta(bRec.wi) <= 0 || Frame::cosTheta(bRec.wo) <= 0) {
                return .0f;
            }

            Vector3f wh = (bRec.wi + bRec.wo).normalized();
            // need to consider distance blending for D
            float D = Warp::squareToBeckmannPdf(wh, m_alpha);
            float Jh = 1.0f/(4*wh.dot(bRec.wo));

            return m_ks * D * Jh + (1-m_ks) * Frame::cosTheta(bRec.wo)/M_PI;
        }

        /// Sample the ƒ
        Color3f sample(BSDFQueryRecord &bRec, const Point2f &_sample) const {
            // Note: Once you have implemented the part that computes the scattered
            // direction, the last part of this function should simply return the
            // BRDF value divided by the solid angle density and multiplied by the
            // cosine factor from the reflection equation, i.e.
            // return eval(bRec) * Frame::cosTheta(bRec.wo) / pdf(bRec);
            bRec.measure = ESolidAngle;
            bRec.eta = 1.0f;

            if (_sample.x() > m_ks) {
                // This is the diffuse case
                Point2f sample = Point2f((_sample.x()-m_ks)/(1-m_ks), _sample.y());
                /* Warp a uniformly distributed sample on [0,1]^2
                to a direction on a cosine-weighted hemisphere */
                bRec.wo = Warp::squareToCosineHemisphere(sample);
            }
            else {
                // This is the specular case
                Point2f sample = Point2f(_sample.x()/m_ks, _sample.y());
                Normal3f normal = Warp::squareToBeckmann(sample, m_alpha);

                // Rotate the incidence ray around the normal for 2 pi radiance
                // This is the same as flipping the incidence ray with respect to the normal
                bRec.wo = -bRec.wi + 2 * normal * normal.dot(bRec.wi);
            }

            if (bRec.measure != ESolidAngle || Frame::cosTheta(bRec.wi) <= 0 || Frame::cosTheta(bRec.wo) <= 0) {
                return Color3f(.0f);
            }

            return eval(bRec) * Frame::cosTheta(bRec.wo) / pdf(bRec);
        }

        /**
            * \brief This function implements the Reflection cache. and Optimal importance sampling.
            * from the paper "A Practical Stochastic Algorithm for Rendering Mirror-Like Flakes"
            * This function should be called inside the sample() function (where an outgoing direction
            * has already been sampled).
            *
            * \param bRec The BSDF record that contains the outgoing direction of its (in the paper this should be the light ray) in local coordinates
            * \param footprint The pixel footprint of the current its in the UV space
            *
            * \return
        */
        float computeDiscreteD(const BSDFQueryRecord &bRec, const Point2f footprint[4]) const {
            return quadTree.evalD(footprint, bRec.wi, bRec.wo, m_alpha);
        }

        bool isDiffuse() const {
            /* While microfacet BRDFs are not perfectly diffuse, they can be
               handled by sampling techniques for diffuse/non-specular materials,
               hence we return true here */
            return true;
        }

        bool needRayDifferentials() const {
            return true;
        }

        string toString() const {
            return tfm::format(
                    "Glints[\n"
                    "  alpha = %f,\n"
                    "  intIOR = %f,\n"
                    "  extIOR = %f,\n"
                    "  kd = %s,\n"
                    "  ks = %f\n"
                    "]",
                    m_alpha,
                    m_intIOR,
                    m_extIOR,
                    m_kd.toString(),
                    m_ks
            );
        }
    private:
        float m_alpha;
        float m_intIOR, m_extIOR;
        float m_ks;
        Color3f m_kd;
        QuadTree quadTree = QuadTree();
    };

    NORI_REGISTER_CLASS(Glints, "glints");
NORI_NAMESPACE_END

quadtree.h

#include "pcg32.h"
#include "common.h"
#include "vector.h"
#include "sampler.h"
#include "warp.h"
#include "bsdf.h"
#include "bbox.h"
#include <stack>

using namespace std;

NORI_NAMESPACE_BEGIN
    struct QuadNode {
        Point2f corners[2];
        size_t num_flakes;
        size_t id;

        QuadNode(const Point2f (&input_corners)[2], const size_t num_flakes, const size_t id): num_flakes(num_flakes), id(id){
            copy(begin(input_corners), end(input_corners), begin(corners));
        }
    };

    class QuadTree {
    public:
        QuadTree() = default;

        /**
         * \brief This function implements the "Texture query." and "Reflection cache."
         * from the paper "A Practical Stochastic Algorithm for Rendering Mirror-Like Flakes"
         *
         * \param footprint the projection of a pixel onto the UV space
         * \param wi The incidence direction of the current light path in the local coordinates system
         * \param wo The outgoing direction of the current light path in the local coordinates system
         * \param m_alpha The roughness parameter of the current BSDF
         *
         * \return The discrete microfacet distribution D
        */
        float evalD(const Point2f footprint[4], const Vector3f &wi, const Vector3f &wo, float m_alpha) const {
            float D = .0f;
            stack<QuadNode> stack; // store quadnodes themselves not ptrs (because not allocating somewhere in memory)
            stack.push(root);

            pcg32 sampler;
            const float cosgamma = cos(GAMMA);

            while (!stack.empty()) {
                QuadNode cur_node = stack.top();
                stack.pop();
                sampler.seed(cur_node.id);

                const float min_x = cur_node.corners[0].x();
                const float max_x = cur_node.corners[1].x();
                const float min_y = cur_node.corners[0].y();
                const float edge_length = max_x - min_x;

                if (cur_node.num_flakes < LEAF_FLAKES || stack.size() > MAX_DEPTH) {
                    // generate exact coordinates
                    // check overlapping with footprint
                    // if overlapping, push the point to the result list
                    for (int counter = 0; counter < cur_node.num_flakes; ++counter){
                        const Point2f p_i = Point2f(min_x + edge_length * sampler.nextFloat(), min_y + edge_length * sampler.nextFloat());
                        const Point2f n = Point2f(sampler.nextFloat(), sampler.nextFloat());

                        if (polygonContainsPoint(footprint, p_i)) {
                            const Normal3f m_i = Warp::squareToBeckmann(n, m_alpha);
                            const Vector3f r_i = -wi + 2 * wi.dot(m_i) * m_i;

                            if (r_i.dot(wo) >= cosgamma) {
                                /// Note that r_i.dot(m_i) is an approximation to the local Beckmann PDF
//                                const float cur_cos_angle = r_i.dot(m_i);
                                const float cur_cos_angle = Warp::squareToBeckmannPdf(m_i, m_alpha);
                                if (cur_cos_angle > .0f) {
                                    D += cur_cos_angle * WEIGHT_SCALE;
                                }
                            }
                        }

                    }
                }
                else {
                    Vector4l sub_flakes;
                    if (cur_node.num_flakes < SAMPLING_FLAKES) {
                        sub_flakes = distributeFlakesNaive(cur_node.num_flakes, sampler);
                    }
                    else {
                        sub_flakes = distributeFlakesMultinomial(cur_node.num_flakes, sampler);
                    }

                    QuadNode children[4] = {
                            QuadNode(
                                    {
                                            Point2f(cur_node.corners[0].x(), cur_node.corners[0].y() + edge_length/2.0f),
                                            Point2f(cur_node.corners[0].x() + edge_length/2.0f, cur_node.corners[1].y()),
                                    },
                                    sub_flakes[0],
                                    (cur_node.id << 2) + 0 // shift 2 bits to the left
                            ),
                            QuadNode(
                                    {
                                            Point2f(cur_node.corners[0].x() + edge_length/2.0f, cur_node.corners[0].y() + edge_length/2.0f),
                                            cur_node.corners[1],
                                    },
                                    sub_flakes[1],
                                    (cur_node.id << 2) + 1 // shift 2 bits to the left
                            ),
                            QuadNode(
                                    {
                                            cur_node.corners[0],
                                            Point2f(cur_node.corners[0].x() + edge_length/2.0f, cur_node.corners[0].y() + edge_length/2.0f),
                                    },
                                    sub_flakes[2],
                                    (cur_node.id << 2) + 2 // shift 2 bits to the left
                            ),
                            QuadNode(
                                    {
                                            Point2f(cur_node.corners[0].x() + edge_length/2.0f, cur_node.corners[0].y()),
                                            Point2f(cur_node.corners[1].x(), cur_node.corners[0].y() + edge_length/2.0f),
                                    },
                                    sub_flakes[3],
                                    (cur_node.id << 2) + 3 // shift 2 bits to the left
                            ),
                    };

                    for (QuadNode child : children) {
                        // if the current child overlaps with the footprint, push onto stack
                        BoundingBox2f cur_box = BoundingBox2f(child.corners[0], child.corners[1]);
                        BoundingBox2f footprint_box = getBboxForFootprint(footprint);
                        if (cur_box.overlaps(footprint_box, true)) {
                            stack.push(child);
                        }
                    }
                }
            }

            D /= N;

            return D;
        }

        /**
         * \brief Get the bounding box of the footprint.
         *
         * \param footprint The footprint of the current iteration.
         *
         * \return The bounding box of the footprint.
         * */
        BoundingBox2f getBboxForFootprint(const Point2f footprint[4]) const {
            float diagonal_dist_1 = (footprint[0] - footprint[2]).norm();
            float diagonal_dist_2 = (footprint[1] - footprint[3]).norm();
            BoundingBox2f footprint_box;
            if (diagonal_dist_1 >= diagonal_dist_2) {
                footprint_box = BoundingBox2f(footprint[0]);
                footprint_box.expandBy(footprint[2]);
            } else {
                footprint_box = BoundingBox2f(footprint[1]);
                footprint_box.expandBy(footprint[3]);
            }

            return footprint_box;
        }

        /**
         * \brief Distribute the number of flakes n to 4 bins according to a multi-dimensional normal
         * distribution with Box-Muller estimation. The distirbution has a mean of p and a covariance
         * of n(diag(p)-pp^T).
         *
         * \param num_flakes The number of flakes to distribute.
         * \param sampler The seeded sampler.
         *
         * \return A vector of 4 integers representing the number of flakes in each bin.
         * */
        Vector4l distributeFlakesMultinomial(size_t num_flakes, pcg32 &sampler) const {
//            const Vector4f p = {0.25,0.25,0.25,0.25};
            // get two pairs of samples from a 2D normal distribution with Box-Muller method
            const Point2f s1 = Point2f(sampler.nextFloat(), sampler.nextFloat());
            const Point2f s2 = Point2f(sampler.nextFloat(), sampler.nextFloat());
            const Vector4f Z = {
                    (float) (sqrt(-2*log(s1.x())) * cos(2 * EIGEN_PI * s1.y())),
                    (float) (sqrt(-2*log(s1.x())) * sin(2 * EIGEN_PI * s1.y())),
                    (float) (sqrt(-2*log(s2.x())) * cos(2 * EIGEN_PI * s2.y())),
                    (float) (sqrt(-2*log(s2.x())) * sin(2 * EIGEN_PI * s2.y()))
            };

            // The covariance matrix of the desired distribution
            MatrixXf sigma = MatrixXf::Constant(4,4,-1.0/8);
            sigma.diagonal() << Vector4f::Constant(3.0/8);
            // transform the samples from mean 0 and unit variance to have mean of p and covariance of sigma
            const Vector4f result = Vector4f::Constant(num_flakes * 0.25) + sqrt(num_flakes) * (sigma * Z);

            // cast the results into integers
            Vector4l result_i = result.cast<size_t>();
            size_t sum = result_i.sum();

            // randomly increase or decrease the number of flakes in one of the bins until the total number is equal to n
            for(; sum < num_flakes; ++sum)
                result_i(std::min(3, (int)(sampler.nextFloat() * 4))) += 1;

            while (sum > num_flakes) {
                int idx = std::min(3, (int)(sampler.nextFloat() * 4));

                if (result_i[idx] > 0) {
                    result_i(idx) -= 1;
                    sum -= 1;
                }
            }

            return result_i;
        }

        /**
         * \brief Distribute the number of flakes n to 4 bins according to a multi-dimensional normal distribution with
         * brute force.
         *
         * \param num_flakes The number of flakes to distribute.
         * \param sampler The seeded sampler.
         *
         * \return A vector of 4 integers representing the number of flakes in each bin.
         * */
        Vector4l distributeFlakesNaive(size_t n, pcg32 &sampler) const {
            Vector4l bins = Vector4l::Constant(0);
            for (size_t i = 0; i < n; ++i) {
                bins(std::min(3, (int)(sampler.nextFloat() * 4))) += 1;
            }

            return bins;
        }

    private:
        const size_t N = 1e5; /// the total number of flakes on the entire object. Need to be size_t to avoid overflow
        const int LEAF_FLAKES = 7; /// the maximum number of flakes per leaf node (need to be smaller than SAMPLING_FLAKES) it's 16 in the paper
        const int MAX_DEPTH = 15; /// the maximum depth for the stack
        const int SAMPLING_FLAKES = 8; /// the maximum number of flakes sampled with distributeFlakesNaive
        const float GAMMA = 0.1f; /// The radius of the cone in radians (in paper it's 1 degree (0.0175 radian) and 5 degrees (0.0873 radian))
        const float WEIGHT_SCALE = 10000000; /// Scale up the contribution of each flake to D
        // Note that root has to be initialised after N cuz
        // otherwise N will be zero
        QuadNode root = QuadNode(
                // order goes: bottom-left, top-right
                {Point2f(0.0f, 0.0f),
                 Point2f(1.0f, 1.0f)},
                N,
                1
        ); /// The root of the node, which is the entire texture space
    };

NORI_NAMESPACE_END

Reference

• Atanasov, A. and Koylazov, V. (2016). A practical stochastic algorithm for rendering mirror-like flakes. In ACM SIGGRAPH 2016 Talks, SIGGRAPH ’16, New York, NY, USA. Association for Computing Machinery.
• Nori