RixRNGProgressive.h

Go to the documentation of this file.

/*
# ------------------------------------------------------------------------------
#
# Copyright (c) 1986-2019 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 RixRNGProgressive_h
#define RixRNGProgressive_h

// IWYU pragma: no_include <assert.h>          // for assert
// IWYU pragma: no_include "RiTypesHelper.h"   // for RtFloat3, RtFloat2
// IWYU pragma: no_include "prmanapi.h"        // for PRMAN_INLINE

// This file is intended solely for inclusion by RixInterfaces.h and should
// not be included directly.  Here we define inline implementations of private
// RixRNG methods.  Conceptually the contents of this file are "opaque" but
// we provide them to deliver the highest performance.  Each release may
// differ from the last, so in order to obtain the new benefits, clients
// will need to recompile.
//
// These routines collectively implement progressive multi-jittered (PMJ) sample
// sequences based on lookups in precomputed tables.  The progressive multi-jittered
// sample sequences are described in detail in:
//   "Progressive multi-jittered sample sequences" by
//   Per Christensen, Andrew Kensler, and Charlie Kilpatrick,
//   Proc. Eurographics Symposium on Rendering, 2018.
// (http://graphics.pixar.com/library/ProgressiveMultiJitteredSampling/).
//
// See also the related work on correlated multi-jittered sample sets described
// in Pixar Technical Memo #13-01, "Correlated Multi-Jittered Sampling"
// (http://graphics.pixar.com/library/MultiJitteredSampling/).
//
// Alternatively, you may use the fields in Ctx as inputs to your own
// sampling sequence instead if you prefer.
//

//
// Progressive multi-jittered samples (PMJ).
// Lookups in pre-generated tables of progressive samples.
// We scramble the bits (xor) -- this preserves the PMJ properties
// but further increases the number of apparently different sample sequences.
// x and y could be flipped in order to add further apparent randomness to
// the samples, thus reducing aliasing from using a limited number of tables
// -- but this alone is not sufficient; bit scrambling gives many more
// different sequences.
// (In contrast, regular Cranley-Patterson rotation is a very bad idea
// for PMJ samples: it increases variance by a lot.  Probably because the
// strata on edges get wrapped around introducing large jumps for the
// samples in those strata.  If we really needed rotation, we could
// restrict the rotations to e.g. multiples of 0.25 or 0.125.)
//

#define N_PMJ_DIMS 24   // number of pmj dimensions (8*3D)
#define N_PMJ_SEQS 48   // number of different pmj sequences
#define N_PMJ_SAMPS 4096   // number of samples in each pmj sequence

class ProgressiveSampler : public Generator
{
public:
    void* pmjTables;

    ProgressiveSampler(void* tables) : pmjTables(tables) {}

private:

    // Compute a repeatable random float in [0,1) based on value and scramble
    PRMAN_INLINE static float
    HashToRandom(const unsigned value, const unsigned scramble)
    {
        // clang-format off
        unsigned result = value;
        result ^= scramble;
        result ^= result >> 17;
        result ^= result >> 10;  result *= 0xb36534e5;
        result ^= result >> 12;
        result ^= result >> 21;  result *= 0x93fc4795;
        result ^= 0xdf6e307f;
        result ^= result >> 17;  result *= 1 | scramble >> 18;
        return static_cast<float>(result) / 4298115584.0f;
        // clang-format on
    }

    // Compute a repeatable random unsigned int based on value and scramble.
    // Same as HashToRandom() but returns unsigned int instead of float.
    // Has avalanche property: if a single input bit changes, an average
    // of half the output bits change.
    PRMAN_INLINE static unsigned int 
    HashToRandomUInt(unsigned int const value, unsigned int const scramble)
    {
        // clang-format off
        unsigned int result = value;
        result ^= scramble;
        result ^= result >> 17;
        result ^= result >> 10;  result *= 0xb36534e5;
        result ^= result >> 12;
        result ^= result >> 21;  result *= 0x93fc4795;
        result ^= 0xdf6e307f;
        result ^= result >> 17;  result *= 1 | scramble >> 18;
        return result;
        // clang-format on
    }

    // Locally shuffle the sample order by carefully swapping neighbor samples
    // (and groups of two samples and reversal of order): only samples of same
    // class are swapped, and only within local groups of 16 samples.
    // Class order is: AABBAABBBBAABBAA ...
    // Local shuffling does not affect the convergence properties of each
    // sequence, but decorrelates 2D sequences so they can be safely combined 
    // to a higher-dimensional sample (aka. "padding").  Even two patterns
    // that happen to map to the same table and dimension will be decorrelated
    // when shuffled differently.
    // (In RenderMan 21 we locally shuffled the sample order by swapping neighbor
    // samples, groups of two samples, groups of four samples, and groups of 
    // eight samples.  Simpler, but breaks multiclass classification.)
    PRMAN_INLINE static unsigned int
    shuffle(unsigned int pattern, unsigned int sample)
    {
        unsigned int sample16, s;
        s = sample;
        sample16 = sample % 16;

        if (sample < 16) return sample;

        // Randomly swap neighbor samples (with same class)
        if (sample16 ==  0 && (pattern & 0x0100)) s++;
        if (sample16 ==  1 && (pattern & 0x0100)) s--;
        if (sample16 ==  2 && (pattern & 0x0200)) s++;
        if (sample16 ==  3 && (pattern & 0x0200)) s--;
        if (sample16 ==  4 && (pattern & 0x0400)) s++;
        if (sample16 ==  5 && (pattern & 0x0400)) s--;
        if (sample16 ==  6 && (pattern & 0x0800)) s++;
        if (sample16 ==  7 && (pattern & 0x0800)) s--;
        if (sample16 ==  8 && (pattern & 0x1000)) s++;
        if (sample16 ==  9 && (pattern & 0x1000)) s--;
        if (sample16 == 10 && (pattern & 0x2000)) s++;
        if (sample16 == 11 && (pattern & 0x2000)) s--;
        if (sample16 == 12 && (pattern & 0x4000)) s++;
        if (sample16 == 13 && (pattern & 0x4000)) s--;
        if (sample16 == 14 && (pattern & 0x8000)) s++;
        if (sample16 == 15 && (pattern & 0x8000)) s--;

        // Randomly swap sample pairs offset by 4 (such pairs have same class)
        if ((sample16 ==  0 || sample16 ==  1) && (pattern & 0x10000)) s += 4;
        if ((sample16 ==  4 || sample16 ==  5) && (pattern & 0x10000)) s -= 4;
        if ((sample16 ==  2 || sample16 ==  3) && (pattern & 0x20000)) s += 4;
        if ((sample16 ==  6 || sample16 ==  7) && (pattern & 0x20000)) s -= 4;
        if ((sample16 ==  8 || sample16 ==  9) && (pattern & 0x40000)) s += 4;
        if ((sample16 == 12 || sample16 == 13) && (pattern & 0x40000)) s -= 4;
        if ((sample16 == 10 || sample16 == 11) && (pattern & 0x80000)) s += 4;
        if ((sample16 == 14 || sample16 == 15) && (pattern & 0x80000)) s -= 4;

        // Randomly reverse sample order within groups of 16 samples
        if (pattern & 0x100000) s = s + 15 - 2*(s%16);

        // Swap adjacent groups of 16 samples (probably not necessary?)
        /*
        if (pattern & 0x200000 && sample >= 32) {
            if ((sample/16) % 2 == 0)
                sample += 16;
            else
                sample -= 16;
        }
        */

        assert(sample/16 == s/16); // original and shuffled sample should be in same group of 16

        return s;
    }

    // Scramble bits of sample f with bits of 'scramble' and convert back to
    // float.  Both the input f and the result are floats between 0.0 and 1.0.
    PRMAN_INLINE static float 
    fpScramble(float f, unsigned int scramble)
    {
        union
        {
            float f;
            unsigned int i;
        } f2i;
        f2i.f = f + 1.0f;         // map f to [1,2) so we know the exponent is 0
        f2i.i ^= (scramble & 0x007fffff); // scramble the fp mantissa bits
        return f2i.f - 1.0f;      // map result to [0,1)
    }

    // Compute a progressive 1D PMJ (actually PJ) sample
    PRMAN_INLINE float progressive1DSample(
        const unsigned int pattern,
        const unsigned int sample) const
    {
        float result;
        float* sampleTable;
        unsigned int d, t, s;
        const unsigned int nDims  = N_PMJ_DIMS;  // 24
        const unsigned int nSeqs  = N_PMJ_SEQS;  // 48
        const unsigned int nSamps = N_PMJ_SAMPS; // 4096
        //const unsigned int totalTableSize = nDims * nSeqs * nSamps;

        if (sample > nSeqs * nSamps)
        {
            // Resort to (repeatable) random when running out of table entries
            // after 48*4096 = 196,608 samples with same pattern id.
            // (Could be 24 times more if wrapping to next dim.)
            result = HashToRandom(sample, pattern * 0x51633e2d);
            return result;
        }

        // Select sequence (table).
        // (Note that we cannot use a 'table' specified by pixel number encoded
        // e.g. in first 6 bits of pattern -- Woodcock tracking needs less 
        // correlation in volume tracing within in each pixel.)
        t = pattern % nSeqs;
        assert(t < nSeqs);

        // Select dimension
        d = (pattern >> 7) % nDims;
        if (d < 3) d += 9; // avoid dimension 0--2: correlated with pixel pos
        assert(d < nDims);

        // Locally shuffle the sample order by carefully swapping neighbor samples
        // (and groups of four samples) -- only samples of same class are swapped.
        s = shuffle(pattern, sample);

        // When we have "used up" all 4096 samples in a table we start over
        // in the next table.
        t += s / nSamps;
        t = t % nSeqs;
        s = s % nSamps;

        // Select table entry -- note that these are ordered with sample number
        // having the largest stride in the table, this is for improved memory
        // access coherency during incremental (progressive) rendering.
        s = s * nSeqs * nDims + t * nDims + d;
        assert(s < nDims * nSeqs * nSamps);

        // Lookup sample in table
        float** const tables = (float**)pmjTables;
        sampleTable = tables[0];   // a single 24D table
        result = sampleTable[s];   // progressive jittered 1D
        assert(0.0f <= result && result <= 1.0f); // 1.0 is okay here

        // Random (but consistent) flips
        // if (pattern & 0x1) result = 1.0f-result;

        // Scramble bits of sample with bits of pattern, convert back to float
        unsigned int pattern1 = HashToRandomUInt(pattern, 0x51633e2d);
        result = fpScramble(result, pattern1);

        // Ensure (bit-scrambled) values are strictly less than 1.0
        if (result > 0.999990f) result = 0.999990f;
        assert(0.0f <= result && result < 1.0f);

        return result;
    }

    // Compute a progressive 2D PMJ sample.
    // Uses bit-scrambling to generate many sequences from only 48 tables.
    PRMAN_INLINE RtFloat2 progressive2DSample(
        const unsigned int pattern,
        const unsigned int sample) const
    {
        RtFloat2 result;
        float* sampleTable;
        unsigned int d, t, s;
        const unsigned int nDims  = N_PMJ_DIMS;  // 24
        const unsigned int nSeqs  = N_PMJ_SEQS;  // 48
        const unsigned int nSamps = N_PMJ_SAMPS; // 4096
        //const unsigned int totalTableSize = nDims * nSeqs * nSamps;

        if (sample > nSeqs * nSamps)
        {
            // Resort to (repeatable) random when running out of table entries
            // after 48*4096 = 196,608 samples with same pattern id.
            // (Could be 8 times more if wrapping to next dim.)
            result.x = HashToRandom(sample, pattern * 0x51633e2d);
            result.y = HashToRandom(sample, pattern * 0x68bc21eb);
            return result;
        }

        // Select sequence (table).
        // (Note that we cannot use a 'table' specified by pixel number encoded
        // e.g. in first 6 bits of pattern -- subsurface scattering domains etc 
        // need less correlation within in each pixel.)
        t = pattern % nSeqs;
        assert(t < nSeqs);

        // Select dimension based on 3 bits of pattern --
        // each pair of dimensions 0+1, 3+4, 6+7, etc. are a pmj (0,2) sequence
        d = 3 * ((pattern >> 23) & 0x7);
        if (d == 0) d = 12; // avoid dimensions 0+1: correlated with pixel pos
        assert(d < nDims);

        // Locally shuffle the sample order by carefully swapping neighbor samples
        // (and groups of four samples) -- only samples of same class are swapped.
        // Local shuffling does not affect the convergence properties of each
        // sequence, but decorrelates 2D sequences so they can be safely combined 
        // to a higher-dimensional sample (aka. "padding").  Even two patterns
        // that happen to map to the same table and dimension will be decorrelated
        // when shuffled differently.
        s = shuffle(pattern, sample);

        // When we have "used up" all 4096 samples in a table we start over
        // in the next table.
        t += s / nSamps;
        t = t % nSeqs;
        s = s % nSamps;

        // Select table entry -- note that these are ordered with sample number
        // having the largest stride in the table, this is for improved memory
        // access coherency during incremental (progressive) rendering.
        s = s * nSeqs * nDims + t * nDims + d;
        assert(s < nDims * nSeqs * nSamps);

        // Lookup sample in table
        float** const tables = (float**)pmjTables;
        sampleTable = tables[0];     // a single 24D table
        result.x = sampleTable[s];   // 2D progressive multijittered (0,2) samples
        result.y = sampleTable[s+1];
        assert(0.0f <= result.x && result.x <= 1.0f); // 1.0 is okay here
        assert(0.0f <= result.y && result.y <= 1.0f);

        // Random (but consistent) flips
        // if (pattern & 0x1) result.x = 1.0f-result.x;
        // if (pattern & 0x2) result.y = 1.0f-result.y;
        // if (pattern & 0x4) std::swap(result.x, result.y);

        // Scramble bits of sample with bits of pattern, convert back to float
        unsigned int pattern1 = HashToRandomUInt(pattern, 0x51633e2d);
        unsigned int pattern2 = HashToRandomUInt(pattern, 0x68bc21eb);
        result.x = fpScramble(result.x, pattern1);
        result.y = fpScramble(result.y, pattern2);

        // Ensure (bit-scrambled) values are strictly less than 1.0
        if (result.x > 0.999990f) result.x = 0.999990f;
        if (result.y > 0.999990f) result.y = 0.999990f;
        assert(0.0f <= result.x && result.x < 1.0f);
        assert(0.0f <= result.y && result.y < 1.0f);

        return result;
    }

    // Compute a progressive 3D PMJ sample
    PRMAN_INLINE RtFloat3 progressive3DSample(
        const unsigned int pattern,
        const unsigned int sample) const
    {
        RtFloat3 result;
        float* sampleTable;
        unsigned int d, t, s;
        const unsigned int nDims  = N_PMJ_DIMS;  // 24
        const unsigned int nSeqs  = N_PMJ_SEQS;  // 48
        const unsigned int nSamps = N_PMJ_SAMPS; // 4096
        //const unsigned int totalTableSize = nDims * nSeqs * nSamps;

        if (sample > nSeqs * nSamps)
        {
            // Resort to (repeatable) random when running out of table entries
            // after 48*4096 = 196,608 samples with same pattern id.
            // (Could be 8 times more if wrapping to next dim.)
            result.x = HashToRandom(sample, pattern * 0x51633e2d);
            result.y = HashToRandom(sample, pattern * 0x68bc21eb);
            result.z = HashToRandom(sample, pattern * 0x02e5be93);
            return result;
        }

        // Select sequence (table).
        // (Note that we cannot use a 'table' specified by pixel number encoded
        // e.g. in first 6 bits of pattern -- subsurface scattering domains etc
        // need less correlation within in each pixel.)
        t = pattern % nSeqs;
        assert(t < nSeqs);

        // Select dimension based on 3 bits of pattern --
        // each triple of dimensions 0--2, 3--5, 6--8, etc. are stratified in 3D
        d = 3 * ((pattern >> 23) & 0x7);
        if (d == 0) d = 12; // avoid dimensions 0+1: correlated with pixel pos
        assert(d < nDims);

        // Locally shuffle the sample order by carefully swapping neighbor samples
        // (and groups of two samples) -- only samples of same class are swapped.
        s = shuffle(pattern, sample);

        // When we have "used up" all 4096 samples in a table we start over
        // in the next table.
        t += s / nSamps;
        t = t % nSeqs;
        s = s % nSamps;

        // Select table entry -- note that these are ordered with sample number
        // having the largest stride in the table, this is for improved memory
        // access coherency during incremental (progressive) rendering.
        s = s * nSeqs * nDims + t * nDims + d;
        assert(s < nDims * nSeqs * nSamps);

        // Lookup sample in table
        float** const tables = (float**)pmjTables;
        sampleTable = tables[0];     // a single 24D table
        result.x = sampleTable[s];   // stratified in 3D
        result.y = sampleTable[s+1];
        result.z = sampleTable[s+2];
        assert(0.0f <= result.x && result.x <= 1.0f); // 1.0 is okay here
        assert(0.0f <= result.y && result.y <= 1.0f);
        assert(0.0f <= result.z && result.z <= 1.0f);

        // Random (but consistent) flips (could also swap xy xz yz)
        // if (pattern & 0x1) result.x = 1.0f-result.x;
        // if (pattern & 0x2) result.y = 1.0f-result.y;
        // if (pattern & 0x4) result.z = 1.0f-result.z;

        // Scramble bits of sample with bits of patterns, convert back to float
        unsigned int pattern1 = HashToRandomUInt(pattern, 0x51633e2d);
        unsigned int pattern2 = HashToRandomUInt(pattern, 0x68bc21eb);
        unsigned int pattern3 = HashToRandomUInt(pattern, 0x02e5be93);
        result.x = fpScramble(result.x, pattern1);
        result.y = fpScramble(result.y, pattern2);
        result.z = fpScramble(result.z, pattern3);

        // Ensure (bit-scrambled) values are strictly less than 1.0
        if (result.x > 0.999990f) result.x = 0.999990f;
        if (result.y > 0.999990f) result.y = 0.999990f;
        if (result.z > 0.999990f) result.z = 0.999990f;
        assert(0.0f <= result.x && result.x < 1.0f);
        assert(0.0f <= result.y && result.y < 1.0f);
        assert(0.0f <= result.z && result.z < 1.0f);

        return result;
    }

public:

    // Functions to draw a single PMJ sample from a sequence.
    // The 'i' parameter (shading context index) is part of the Generator class 
    // interface, but ignored for PMJ samples.
    virtual float Sample1D(const SampleCtx& rctx, unsigned i) const
    {
        PIXAR_ARGUSED(i);
        return progressive1DSample(rctx.patternid, rctx.sampleid);
    }

    virtual RtFloat2 Sample2D(const SampleCtx& rctx, unsigned i) const
    {
        PIXAR_ARGUSED(i);
        return progressive2DSample(rctx.patternid, rctx.sampleid);
    }

    virtual RtFloat3 Sample3D(const SampleCtx& rctx, unsigned i) const
    {
        PIXAR_ARGUSED(i);
        return progressive3DSample(rctx.patternid, rctx.sampleid);
    }

    // Generator functions to draw multiple samples, each from a different sequence.
    // (Candidates for for future simd implementation)
    virtual void MultiSample1D(
        unsigned int n,
        const SampleCtx* rctx,
        float* xis) const
    {
        for (unsigned int i = 0; i < n; ++i)
        {
            xis[i] = progressive1DSample(rctx[i].patternid, rctx[i].sampleid);
        }
    }
    virtual void MultiSample2D(
        unsigned int n,
        const SampleCtx* rctx,
        RtFloat2* xis) const
    {
        for (unsigned int i = 0; i < n; ++i)
        {
            xis[i] = progressive2DSample(rctx[i].patternid, rctx[i].sampleid);
        }
    }
    virtual void MultiSample3D(
        unsigned int n,
        const SampleCtx* rctx,
        RtFloat3* xis) const
    {
        for (unsigned int i = 0; i < n; ++i)
        {
            xis[i] = progressive3DSample(rctx[i].patternid, rctx[i].sampleid);
        }
    }
};

#endif