mmp / pbrt-v2

Source code for the version of pbrt described in the second edition of "Physically Based Rendering"
http://pbrt.org
989 stars 343 forks source link

PermutedRadicalInverse() is horribly broken #52

Closed mmp closed 8 years ago

mmp commented 8 years ago

Per Christensen noted the following. (Text slightly edited.)

I came across an issue with the PermutedRadicalInverse() function from the 2nd edition. To explain the issue as simply as possible, assume we’re computing in base 2 and the permutation array p is {1,0} so ones get converted to zeroes and vice versa. (I know that for the particular case of base 2 there are much more efficient ways of doing digit scrambling using e.g. bit-wise xor, but let’s ignore that for now.)

Here’s my version of your function — pretty much verbatim straight out of the book:

float 
PermutedRadicalInverse(unsigned n, const unsigned base, unsigned *perm)
{
    float val = 0.0f;
    float invBase = 1.0f / base, invBi = invBase;
    while (n > 0) {
        unsigned d_i = perm[n%base];
        val += d_i * invBi;
        n *= invBase;
        invBi *= invBase;
    }
    return val;
}

The results I get for n=0..15 are not what I’d expect:

0 0.000000
1 0.000000
2 0.500000
3 0.000000
4 0.750000
5 0.250000
6 0.500000
7 0.000000
8 0.875000
9 0.375000
10 0.625000
11 0.125000
12 0.750000
13 0.250000
14 0.500000
15 0.000000

There are several repeats of 0 and 0.5. It seems to me that the issue is that the loop exits too early — just because n is 0 doesn’t mean we can skip the rest of the digits. Here is a slightly modified version that loops over all digits to give them a chance to be permuted into non-zero values:

float 
PermutedRadicalInverseBetter(unsigned n, const unsigned base, unsigned *perm)
{
    float val = 0.0f;
    float invBase = 1.0f / base, invBi = invBase;
    for (unsigned d = 0; d < 32; d++) { // loop over all digits
        unsigned d_i = perm[n%base];
        val += d_i * invBi;
        n *= invBase;
        invBi *= invBase;
    }
    return val;
}

With this change I get the following values (which seem more reasonable):

0 1.000000
1 0.500000
2 0.750000
3 0.250000
4 0.875000
5 0.375000
6 0.625000
7 0.125000
8 0.937500
9 0.437500
10 0.687500
11 0.187500
12 0.812500
13 0.312500
14 0.562500
15 0.062500