intern/cycles/kernel/kernel_random.h

   1 /*
   2  * Copyright 2011-2013 Blender Foundation
   3  *
   4  * Licensed under the Apache License, Version 2.0 (the "License");
   5  * you may not use this file except in compliance with the License.
   6  * You may obtain a copy of the License at
   7  *
   8  * http://www.apache.org/licenses/LICENSE-2.0
   9  *
  10  * Unless required by applicable law or agreed to in writing, software
  11  * distributed under the License is distributed on an "AS IS" BASIS,
  12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  13  * See the License for the specific language governing permissions and
  14  * limitations under the License.
  15  */
  16
  17 #include "kernel_jitter.h"
  18
  19 CCL_NAMESPACE_BEGIN
  20
  21 #ifdef __SOBOL__
  22
  23 /* skip initial numbers that are not as well distributed, especially the
  24  * first sequence is just 0 everywhere, which can be problematic for e.g.
  25  * path termination */
  26 #define SOBOL_SKIP 64
  27
  28 /* High Dimensional Sobol */
  29
  30 /* van der corput radical inverse */
  31 ccl_device uint van_der_corput(uint bits)
  32 {
  33         bits = (bits << 16) | (bits >> 16);
  34         bits = ((bits & 0x00ff00ff) << 8) | ((bits & 0xff00ff00) >> 8);
  35         bits = ((bits & 0x0f0f0f0f) << 4) | ((bits & 0xf0f0f0f0) >> 4);
  36         bits = ((bits & 0x33333333) << 2) | ((bits & 0xcccccccc) >> 2);
  37         bits = ((bits & 0x55555555) << 1) | ((bits & 0xaaaaaaaa) >> 1);
  38         return bits;
  39 }
  40
  41 /* sobol radical inverse */
  42 ccl_device uint sobol(uint i)
  43 {
  44         uint r = 0;
  45
  46         for(uint v = 1U << 31; i; i >>= 1, v ^= v >> 1)
  47                 if(i & 1)
  48                         r ^= v;
  49
  50         return r;
  51 }
  52
  53 /* inverse of sobol radical inverse */
  54 ccl_device uint sobol_inverse(uint i)
  55 {
  56         const uint msb = 1U << 31;
  57         uint r = 0;
  58
  59         for(uint v = 1; i; i <<= 1, v ^= v << 1)
  60                 if(i & msb)
  61                         r ^= v;
  62
  63         return r;
  64 }
  65
  66 /* multidimensional sobol with generator matrices
  67  * dimension 0 and 1 are equal to van_der_corput() and sobol() respectively */
  68 ccl_device uint sobol_dimension(KernelGlobals *kg, int index, int dimension)
  69 {
  70         uint result = 0;
  71         uint i = index;
  72
  73         for(uint j = 0; i; i >>= 1, j++)
  74                 if(i & 1)
  75                         result ^= kernel_tex_fetch(__sobol_directions, 32*dimension + j);
  76
  77         return result;
  78 }
  79
  80 /* lookup index and x/y coordinate, assumes m is a power of two */
  81 ccl_device uint sobol_lookup(const uint m, const uint frame, const uint ex, const uint ey, uint *x, uint *y)
  82 {
  83         /* shift is constant per frame */
  84         const uint shift = frame << (m << 1);
  85         const uint sobol_shift = sobol(shift);
  86         /* van der Corput is its own inverse */
  87         const uint lower = van_der_corput(ex << (32 - m));
  88         /* need to compensate for ey difference and shift */
  89         const uint sobol_lower = sobol(lower);
  90         const uint mask = ~-(1 << m) << (32 - m); /* only m upper bits */
  91         const uint delta = ((ey << (32 - m)) ^ sobol_lower ^ sobol_shift) & mask;
  92         /* only use m upper bits for the index (m is a power of two) */
  93         const uint sobol_result = delta | (delta >> m);
  94         const uint upper = sobol_inverse(sobol_result);
  95         const uint index = shift | upper | lower;
  96         *x = van_der_corput(index);
  97         *y = sobol_shift ^ sobol_result ^ sobol_lower;
  98         return index;
  99 }
 100
 101 ccl_device_inline float path_rng_1D(KernelGlobals *kg, ccl_addr_space RNG *rng, int sample, int num_samples, int dimension)
 102 {
 103 #ifdef __CMJ__
 104         if(kernel_data.integrator.sampling_pattern == SAMPLING_PATTERN_CMJ) {
 105                 /* correlated multi-jittered */
 106                 int p = *rng + dimension;
 107                 return cmj_sample_1D(sample, num_samples, p);
 108         }
 109 #endif
 110
 111 #ifdef __SOBOL_FULL_SCREEN__
 112         uint result = sobol_dimension(kg, *rng, dimension);
 113         float r = (float)result * (1.0f/(float)0xFFFFFFFF);
 114         return r;
 115 #else
 116         /* compute sobol sequence value using direction vectors */
 117         uint result = sobol_dimension(kg, sample + SOBOL_SKIP, dimension);
 118         float r = (float)result * (1.0f/(float)0xFFFFFFFF);
 119
 120         /* Cranly-Patterson rotation using rng seed */
 121         float shift;
 122
 123         /* using the same *rng value to offset seems to give correlation issues,
 124          * we could hash it with the dimension but this has a performance impact,
 125          * we need to find a solution for this */
 126         if(dimension & 1)
 127                 shift = (*rng >> 16) * (1.0f/(float)0xFFFF);
 128         else
 129                 shift = (*rng & 0xFFFF) * (1.0f/(float)0xFFFF);
 130
 131         return r + shift - floorf(r + shift);
 132 #endif
 133 }
 134
 135 ccl_device_inline void path_rng_2D(KernelGlobals *kg, ccl_addr_space RNG *rng, int sample, int num_samples, int dimension, float *fx, float *fy)
 136 {
 137 #ifdef __CMJ__
 138         if(kernel_data.integrator.sampling_pattern == SAMPLING_PATTERN_CMJ) {
 139                 /* correlated multi-jittered */
 140                 int p = *rng + dimension;
 141                 cmj_sample_2D(sample, num_samples, p, fx, fy);
 142         }
 143         else
 144 #endif
 145         {
 146                 /* sobol */
 147                 *fx = path_rng_1D(kg, rng, sample, num_samples, dimension);
 148                 *fy = path_rng_1D(kg, rng, sample, num_samples, dimension + 1);
 149         }
 150 }
 151
 152 ccl_device_inline void path_rng_init(KernelGlobals *kg, ccl_global uint *rng_state, int sample, int num_samples, ccl_addr_space RNG *rng, int x, int y, float *fx, float *fy)
 153 {
 154 #ifdef __SOBOL_FULL_SCREEN__
 155         uint px, py;
 156         uint bits = 16; /* limits us to 65536x65536 and 65536 samples */
 157         uint size = 1 << bits;
 158         uint frame = sample;
 159
 160         *rng = sobol_lookup(bits, frame, x, y, &px, &py);
 161
 162         *rng ^= kernel_data.integrator.seed;
 163
 164         if(sample == 0) {
 165                 *fx = 0.5f;
 166                 *fy = 0.5f;
 167         }
 168         else {
 169                 *fx = size * (float)px * (1.0f/(float)0xFFFFFFFF) - x;
 170                 *fy = size * (float)py * (1.0f/(float)0xFFFFFFFF) - y;
 171         }
 172 #else
 173         *rng = *rng_state;
 174
 175         *rng ^= kernel_data.integrator.seed;
 176
 177         if(sample == 0) {
 178                 *fx = 0.5f;
 179                 *fy = 0.5f;
 180         }
 181         else {
 182                 path_rng_2D(kg, rng, sample, num_samples, PRNG_FILTER_U, fx, fy);
 183         }
 184 #endif
 185 }
 186
 187 ccl_device void path_rng_end(KernelGlobals *kg, ccl_global uint *rng_state, RNG rng)
 188 {
 189         /* nothing to do */
 190 }
 191
 192 #else
 193
 194 /* Linear Congruential Generator */
 195
 196 ccl_device_inline float path_rng_1D(KernelGlobals *kg, RNG& rng, int sample, int num_samples, int dimension)
 197 {
 198         /* implicit mod 2^32 */
 199         rng = (1103515245*(rng) + 12345);
 200         return (float)rng * (1.0f/(float)0xFFFFFFFF);
 201 }
 202
 203 ccl_device_inline void path_rng_2D(KernelGlobals *kg, RNG& rng, int sample, int num_samples, int dimension, float *fx, float *fy)
 204 {
 205         *fx = path_rng_1D(kg, rng, sample, num_samples, dimension);
 206         *fy = path_rng_1D(kg, rng, sample, num_samples, dimension + 1);
 207 }
 208
 209 ccl_device void path_rng_init(KernelGlobals *kg, ccl_global uint *rng_state, int sample, int num_samples, RNG *rng, int x, int y, float *fx, float *fy)
 210 {
 211         /* load state */
 212         *rng = *rng_state;
 213
 214         *rng ^= kernel_data.integrator.seed;
 215
 216         if(sample == 0) {
 217                 *fx = 0.5f;
 218                 *fy = 0.5f;
 219         }
 220         else {
 221                 path_rng_2D(kg, rng, sample, num_samples, PRNG_FILTER_U, fx, fy);
 222         }
 223 }
 224
 225 ccl_device void path_rng_end(KernelGlobals *kg, ccl_global uint *rng_state, RNG rng)
 226 {
 227         /* store state for next sample */
 228         *rng_state = rng;
 229 }
 230
 231 #endif
 232
 233 /* Linear Congruential Generator */
 234
 235 ccl_device uint lcg_step_uint(ccl_addr_space uint *rng)
 236 {
 237         /* implicit mod 2^32 */
 238         *rng = (1103515245*(*rng) + 12345);
 239         return *rng;
 240 }
 241
 242 ccl_device float lcg_step_float(ccl_addr_space uint *rng)
 243 {
 244         /* implicit mod 2^32 */
 245         *rng = (1103515245*(*rng) + 12345);
 246         return (float)*rng * (1.0f/(float)0xFFFFFFFF);
 247 }
 248
 249 ccl_device uint lcg_init(uint seed)
 250 {
 251         uint rng = seed;
 252         lcg_step_uint(&rng);
 253         return rng;
 254 }
 255
 256 /* Path Tracing Utility Functions
 257  *
 258  * For each random number in each step of the path we must have a unique
 259  * dimension to avoid using the same sequence twice.
 260  *
 261  * For branches in the path we must be careful not to reuse the same number
 262  * in a sequence and offset accordingly. */
 263
 264 ccl_device_inline float path_state_rng_1D(KernelGlobals *kg, ccl_addr_space RNG *rng, const ccl_addr_space PathState *state, int dimension)
 265 {
 266         return path_rng_1D(kg, rng, state->sample, state->num_samples, state->rng_offset + dimension);
 267 }
 268
 269 ccl_device_inline float path_state_rng_1D_for_decision(KernelGlobals *kg, ccl_addr_space RNG *rng, const ccl_addr_space PathState *state, int dimension)
 270 {
 271         /* the rng_offset is not increased for transparent bounces. if we do then
 272          * fully transparent objects can become subtly visible by the different
 273          * sampling patterns used where the transparent object is.
 274          *
 275          * however for some random numbers that will determine if we next bounce
 276          * is transparent we do need to increase the offset to avoid always making
 277          * the same decision */
 278         int rng_offset = state->rng_offset + state->transparent_bounce*PRNG_BOUNCE_NUM;
 279         return path_rng_1D(kg, rng, state->sample, state->num_samples, rng_offset + dimension);
 280 }
 281
 282 ccl_device_inline void path_state_rng_2D(KernelGlobals *kg, ccl_addr_space RNG *rng, const ccl_addr_space PathState *state, int dimension, float *fx, float *fy)
 283 {
 284         path_rng_2D(kg, rng, state->sample, state->num_samples, state->rng_offset + dimension, fx, fy);
 285 }
 286
 287 ccl_device_inline float path_branched_rng_1D(KernelGlobals *kg, ccl_addr_space RNG *rng, const PathState *state, int branch, int num_branches, int dimension)
 288 {
 289         return path_rng_1D(kg, rng, state->sample*num_branches + branch, state->num_samples*num_branches, state->rng_offset + dimension);
 290 }
 291
 292 ccl_device_inline float path_branched_rng_1D_for_decision(KernelGlobals *kg, ccl_addr_space RNG *rng, const PathState *state, int branch, int num_branches, int dimension)
 293 {
 294         int rng_offset = state->rng_offset + state->transparent_bounce*PRNG_BOUNCE_NUM;
 295         return path_rng_1D(kg, rng, state->sample*num_branches + branch, state->num_samples*num_branches, rng_offset + dimension);
 296 }
 297
 298 ccl_device_inline void path_branched_rng_2D(KernelGlobals *kg, ccl_addr_space RNG *rng, const PathState *state, int branch, int num_branches, int dimension, float *fx, float *fy)
 299 {
 300         path_rng_2D(kg, rng, state->sample*num_branches + branch, state->num_samples*num_branches, state->rng_offset + dimension, fx, fy);
 301 }
 302
 303 ccl_device_inline void path_state_branch(PathState *state, int branch, int num_branches)
 304 {
 305         /* path is splitting into a branch, adjust so that each branch
 306          * still gets a unique sample from the same sequence */
 307         state->rng_offset += PRNG_BOUNCE_NUM;
 308         state->sample = state->sample*num_branches + branch;
 309         state->num_samples = state->num_samples*num_branches;
 310 }
 311
 312 ccl_device_inline uint lcg_state_init(RNG *rng, const ccl_addr_space PathState *state, uint scramble)
 313 {
 314         return lcg_init(*rng + state->rng_offset + state->sample*scramble);
 315 }
 316
 317 CCL_NAMESPACE_END
 318