RixRNGProgressive.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 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 a pre-generated table of progressive samples.
// We xor-scramble the bits and locally shuffle the order -- this preserves the PMJ
// properties but increases the number of apparently different sample sequences.
// The xor-scrambles and shuffles have been carefully generated to give
// visually pleasing blue noise noise patterns and also improve denoising.
//

#define N_PMJ_SEQS (256*256)   // number of different shuffles+scrambles of pmj sequences
#define N_PMJ_SAMPS 4096   // number of samples in the 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
    }

    // Shuffle the sample order to avoid correlation between samples when sequences are 
    // combined for high-dimensional sampling ("padding").  
    // Local shuffling does not affect the convergence properties of each sequence.
    // This shuffling involves a fine balance: If the samples are not shuffled enough, then
    // the correlation will persist; if the samples are shuffled too much then there is too
    // much noise at samplecounts between powers of two (as in regression tests with e.g.
    // 49 or 900 spp).
    // Here the shuffle tables have been generated so that shuffling is within groups of
    // at most 32, while only allowing swaps within each of the two A/B subclasses.
    // Class order is: AABBAABBBBAABBAA ...
    // (Alternatively, this shuffling could be done using Brent's hashed Owen-shuffling (as
    // long as we keep track of the sample class the sample had before the shuffling); this
    // might better preserve the lowest possible noise for in-between sample counts?)
    PRMAN_INLINE unsigned int
    shuffle(const unsigned int pattern,
            const unsigned int sample) const
    {
        void** const tables = (void**)pmjTables;
        unsigned short* shuffleTable = (unsigned short*)tables[2];  // extra shuffle table
        unsigned int t1, t2, s;
        const unsigned int nSamps = N_PMJ_SAMPS;  // 4096

        // Optimization: first 4 samples are never shuffled, so don't bother looking up in shuffleTable
        if (sample < 4) return sample;

        t1 = pattern & 0xff;
        t2 = (pattern >> 8) & 0xff;
        assert(sample < nSamps);
        s = sample;

        // Do two shuffles by looking up in shuffle table twice -- this gives us a total of
        // 256^2 = 64K different shuffles.
        // (Ideas for future improvement:
        // - The shuffle table could be smaller by storing signed char deltas between -31 and 31
        //   instead of unsigned shorts.  This would require some signed-unsigned conversions here.
        // - The table reads could be more coherent by storing all sample 0s first, then sample 1s,
        //   etc.)
        s = shuffleTable[t1 * nSamps + s];
        s = shuffleTable[t2 * nSamps + s];

        assert(s/32 == sample/32);  // shuffles are local within groups of 32 samples
        assert(s < nSamps);
        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 >> 9); // scramble the 23 fp mantissa bits
        return f2i.f - 1.0f;      // map result to [0,1)
    }

    // Compute a progressive 1D PMJ (actually PJ) sample.
    // Uses shifting of sequences and shuffling and scrambling of samples to generate more
    // sequences out of the original pj sequence in the table, while at the same time
    // preserving the tiled blue noise properties of the samples.
    PRMAN_INLINE float progressive1DSample(
        const unsigned int pattern,
        const unsigned int sample) const
    {
        void** const tables = (void**)pmjTables;
        float* sampleTable = (float*)tables[0];  // pmj02+pj table
        unsigned int* shuffleScrambleTable = (unsigned int*)tables[1];  // shuffle+scramble table
        float result;
        unsigned int t, shift, s, ss;
        unsigned int shuffleBitsZ, scrambleBitsZ;
        const unsigned int nSeqs = N_PMJ_SEQS;  // 256*256 = 65,536
        const unsigned int nSamps = N_PMJ_SAMPS;  // 4096

        if (sample > nSeqs * nSamps)
        {
            // Resort to (repeatable) random when running out of table entries
            // after 256*256*4096 = 67,108,864 samples with same pattern id.
            result = HashToRandom(sample, pattern * 0x51633e2d);
            return result;
        }

        // Select sequence (table) based on the pixel position encoded in the top 16 bits and
        // a "random" shift that depends on the domain (and indirectly on the frame number, to
        // avoid shower-door effect)
        t = pattern >> 16;  // for 256x256 tiles (8+8 bits for pixel pos leaves 16 bits)
        assert(t < nSeqs);
        shift = pattern & 0xffff;  // shift sequence (table) to avoid shower-door effect
        t += shift;  // (% nSeqs is done below)

        // When we have "used up" all 4096 samples in a pixel's sequence we start over
        // in another sequence 10 rows down, 19 columns over.  (2579 is a prime.)
        t += 2579 * (sample / nSamps);
        t = t % nSeqs;
        s = sample % nSamps;

        // Locally pre-shuffle the sample order by carefully swapping samples
        // within groups of up to 32.  Only samples of same class are swapped.
        s = shuffle(pattern, s);

        // Lookup z shuffle and scramble bits from large table for tile
        shuffleBitsZ = shuffleScrambleTable[5*t+3];
        scrambleBitsZ = shuffleScrambleTable[5*t+4];

        // Shuffle the sample number
        assert(shuffleBitsZ < 64);  // we only shuffle z by at most 64 (to preserve 3D stratification)
        ss = s ^ shuffleBitsZ;  // xor the sample index
        assert(ss < nSamps);

        // Lookup sample in pj table
        result = sampleTable[3*ss+2];  // z dim: pj sequence with 1D tiled blue noise
        assert(0.0f <= result && result < 1.0f);

        // Scramble sample value
        result = fpScramble(result, scrambleBitsZ);  // xor the sample value

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

        // Ensure scrambled (and flipped?) value is strictly less than 1.0
        result = std::min(result, 0.999990f);
        assert(0.0f <= result && result < 1.0f);

        return result;
    }

    // Compute a tiled 2D PMJ sample.
    // Uses shifting of sequences and shuffling and scrambling of samples to generate more
    // sequences out of the original pmj02 sequence in the table, while at the same time
    // preserving the tiled blue noise properties of the samples.
    // (This function is not suitable for photon emission, where the 64K sequences -- even
    // with 8 "random" flips -- are not enough unique samples.)
    PRMAN_INLINE RtFloat2 progressive2DSample(
        const unsigned int pattern,
        const unsigned int sample) const
    {
        void** const tables = (void**)pmjTables;
        float* sampleTable = (float*)tables[0];  // pmj02+pj table
        unsigned int* shuffleScrambleTable = (unsigned int*)tables[1];  // shuffle+scramble table
        RtFloat2 result;
        unsigned int t, shift, s, ss;
        unsigned int shuffleBits, scrambleBitsX, scrambleBitsY;
        const unsigned int nSeqs = N_PMJ_SEQS;  // 256*256 = 65,536
        const unsigned int nSamps = N_PMJ_SAMPS;  // 4096

        if (sample > nSeqs * nSamps)
        {
            // Resort to (repeatable) random when running out of table entries
            // after 256*256*4096 = 67,108,864 samples with same pattern id.
            result.x = HashToRandom(sample, pattern * 0x51633e2d);
            result.y = HashToRandom(sample, pattern * 0x68bc21eb);
            return result;
        }

        // Select sequence (table) based on the pixel position encoded in the top 16 bits and
        // a "random" shift that depends on the domain (and indirectly on the frame number, to
        // avoid shower-door effect)
        t = pattern >> 16;  // for 256x256 tiles (8+8 bits for pixel pos leaves 16 bits)
        assert(t < nSeqs);
        shift = pattern & 0xffff;  // shift sequence to avoid shower-door effect
        t += shift;  // (% nSeqs is done below)

        // When we have "used up" all 4096 samples in a pixel's sequence we start over
        // in another sequence 10 rows down, 19 columns over.  (2579 is a prime.)
        t += 2579 * (sample / nSamps);
        t = t % nSeqs;
        s = sample % nSamps;

        // Locally pre-shuffle the sample order by carefully swapping samples
        // within groups of up to 32.  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 sequence will be decorrelated when shuffled
        // differently.
        s = shuffle(pattern, s);

        // Lookup shuffle and scramble bits from large table for tile
        shuffleBits = shuffleScrambleTable[5*t];
        scrambleBitsX = shuffleScrambleTable[5*t+1];
        scrambleBitsY = shuffleScrambleTable[5*t+2];

        // Shuffle the sample number
        assert(shuffleBits < 256);
        ss = s ^ shuffleBits;  // xor the sample index
        assert(ss < nSamps);

        // Lookup sample in pmj02 table
        result.x = sampleTable[3*ss];
        result.y = sampleTable[3*ss+1];
        assert(0.0f <= result.x && result.x < 1.0f);
        assert(0.0f <= result.y && result.y < 1.0f);

        // Scramble sample x and y values
        result.x = fpScramble(result.x, scrambleBitsX);  // xor the sample's x value
        result.y = fpScramble(result.y, scrambleBitsY);  // xor the sample's y value

        // Random (but consistent) flips -- these preserve tiled blue noise (but useless??)
        // Using 0xX000 as scrables makes quality worse; using 0xX0000 does nothing to quality
        //if (pattern & 0x80000) result.x = 1.0f - result.x;
        //if (pattern & 0x40000) result.y = 1.0f - result.y;
        //if (pattern & 0x20000) std::swap(result.x, result.y);

        // Ensure scrambled (and flipped?) values are strictly less than 1.0
        result.x = std::min(result.x, 0.999990f);
        result.y = std::min(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 scrambled progressive 2D PMJ sample.
    // Nearly identical to the function above, but this version does not only use the 16 top
    // bits of pattern to select the table, and it has a last step that does
    // bit-scrambling to generate *many* different sequences from the pmj02 table.
    // This is equivalent to the old progressive2DSample() function from RenderMan 23 and older.
    // This function is used for photon emission, where the original 256*256 sequences -- even
    // with 8 "random" flips = 512K sequences -- do not provide enough unique samples.
    PRMAN_INLINE RtFloat2 progressive2DScrambledSample(
        const unsigned int pattern,
        const unsigned int sample) const
    {
        void** const tables = (void**)pmjTables;
        float* sampleTable = (float*)tables[0];  // pmj02+pj table
        unsigned int* shuffleScrambleTable = (unsigned int*)tables[1];  // shuffle+scramble table
        RtFloat2 result;
        unsigned int t, s, ss;
        unsigned int shuffleBits, pattern1, pattern2;
        const unsigned int nSeqs = N_PMJ_SEQS;  // 256*256 = 65,536
        const unsigned int nSamps = N_PMJ_SAMPS;  // 4096

        if (sample > nSeqs * nSamps)
        {
            // Resort to (repeatable) random when running out of table entries
            // after 256*256*4096 = 67,108,864 samples with same pattern id.
            result.x = HashToRandom(sample, pattern * 0x51633e2d);
            result.y = HashToRandom(sample, pattern * 0x68bc21eb);
            return result;
        }

        // Select sequence (table)
        t = pattern % nSeqs;  // same as in RenderMan 23 and older
        assert(t < nSeqs);

        // When we have "used up" all 4096 samples in a pixel's sequence we start over
        // in another sequence 10 rows down, 19 columns over.  (2579 is a prime.)
        t += 2579 * (sample / nSamps);
        t = t % nSeqs;
        s = sample % nSamps;

        // Locally pre-shuffle the sample order by carefully swapping samples
        // within groups of up to 32.  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 sequence will be decorrelated when shuffled
        // differently.
        s = shuffle(pattern, s);

        // Lookup shuffle bits from large table for tile (ignore scramble bits in table)
        shuffleBits = shuffleScrambleTable[5*t];

        // Shuffle the sample number
        ss = s ^ shuffleBits;  // xor the sample index
        assert(ss < nSamps);

        // Lookup sample in pmj02 table
        result.x = sampleTable[3*ss];
        result.y = sampleTable[3*ss+1];
        assert(0.0f <= result.x && result.x < 1.0f);
        assert(0.0f <= result.y && result.y < 1.0f);

        // Xor-scramble bits of sample with bits of pattern, convert back to float
        pattern1 = HashToRandomUInt(pattern, 0x51633e2d);
        pattern2 = HashToRandomUInt(pattern, 0x68bc21eb);
        result.x = fpScramble(result.x, pattern1);  // xor the sample's x value
        result.y = fpScramble(result.y, pattern2);  // xor the sample's y value

        // Ensure (bit-scrambled) values are strictly less than 1.0
        result.x = std::min(result.x, 0.999990f);
        result.y = std::min(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.
    // Uses shifting of sequences and shuffling and scrambling of samples to generate more
    // sequences out of the original pmj02+pj 3D sequence in the table, while at the same time
    // preserving the tiled blue noise properties of the samples.
    PRMAN_INLINE RtFloat3 progressive3DSample(
        const unsigned int pattern,
        const unsigned int sample) const
    {
        void** const tables = (void**)pmjTables;
        float* sampleTable = (float*)tables[0];  // pmj02+pj table
        unsigned int* shuffleScrambleTable = (unsigned int*)tables[1];  // shuffle+scramble table
        RtFloat3 result;
        unsigned int t, shift, s, ss, sz;
        unsigned int shuffleBits, scrambleBitsX, scrambleBitsY, shuffleBitsZ, scrambleBitsZ;
        const unsigned int nSeqs = N_PMJ_SEQS;  // 256*256 = 65,536
        const unsigned int nSamps = N_PMJ_SAMPS;  // 4096

        if (sample > nSeqs * nSamps)
        {
            // Resort to (repeatable) random when running out of table entries
            // after 256*256*4096 = 67,108,864 samples with same pattern id.
            result.x = HashToRandom(sample, pattern * 0x51633e2d);
            result.y = HashToRandom(sample, pattern * 0x68bc21eb);
            result.z = HashToRandom(sample, pattern * 0x02e5be93);
            return result;
        }

        // Select sequence (table) based on the pixel position encoded in the top 16 bits and
        // a "random" shift that depends on the domain (and indirectly on the frame number, to
        // avoid shower-door effect)
        t = pattern >> 16;  // for 256x256 tiles (8+8 bits for pixel pos)
        assert(t < nSeqs);
        shift = pattern & 0xffff;  // shift sequence (table) to avoid shower-door effect
        t += shift;  // (% nSeqs is done below)

        // When we have "used up" all 4096 samples in a pixel's sequence we start over
        // in another sequence 10 rows down, 19 columns over.  (2579 is a prime.)
        t += 2579 * (sample / nSamps);
        t = t % nSeqs;
        s = sample % nSamps;

        // Locally pre-shuffle the sample order by carefully swapping samples
        // within groups of up to 32.  Only samples of same class are swapped.
        s = shuffle(pattern, s);

        // Lookup shuffle and scramble bits from large table for tile
        shuffleBits = shuffleScrambleTable[5*t];
        scrambleBitsX = shuffleScrambleTable[5*t+1];
        scrambleBitsY = shuffleScrambleTable[5*t+2];
        shuffleBitsZ = shuffleScrambleTable[5*t+3];
        scrambleBitsZ = shuffleScrambleTable[5*t+4];

        // Shuffle the sample number (z has different shuffle than x and y)
        ss = s ^ shuffleBits;  // xor the xy sample index
        assert(ss < nSamps);
        sz = s ^ shuffleBitsZ;  // xor the z sample index
        assert(sz < nSamps);

        // Lookup sample in pmj02+pj 3D table
        result.x = sampleTable[3*ss];  // progressive multijittered (0,2) samples + z strafied in 3D
        result.y = sampleTable[3*ss+1];
        result.z = sampleTable[3*sz+2];
        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);

        // Scramble sample x, y, and z values
        result.x = fpScramble(result.x, scrambleBitsX);  // xor the sample's x value
        result.y = fpScramble(result.y, scrambleBitsY);  // xor the sample's y value
        result.z = fpScramble(result.z, scrambleBitsZ);  // xor the sample's z value

        // Random (but consistent) flips.  (Could also swap xy xz yz.)
        //if (pattern & 0x80000) result.x = 1.0f - result.x;
        //if (pattern & 0x40000) result.y = 1.0f - result.y;
        //if (pattern & 0x20000) result.z = 1.0f - result.z;

        // Ensure scramble (and flipped?) values are strictly less than 1.0
        result.x = std::min(result.x, 0.999990f);
        result.y = std::min(result.y, 0.999990f);
        result.z = std::min(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 RtFloat2 ScrambledSample2D(const SampleCtx& rctx, unsigned i) const
    {
        PIXAR_ARGUSED(i);
        return progressive2DScrambledSample(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.
    // (These functions are candidates 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 MultiScrambledSample2D(
        unsigned int n,
        const SampleCtx* rctx,
        RtFloat2* xis) const
    {
        for (unsigned int i = 0; i < n; ++i)
        {
            xis[i] = progressive2DScrambledSample(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