b534d9950c552c96749222633fd0a7219b77c9ae
[blender-staging.git] / 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 /* Temporary workaround for Pascal cards, otherwise AA does not work properly. */
136 #if defined(__KERNEL_GPU__) && __CUDA_ARCH__ >= 600
137 __device__ __forceinline__
138 #else
139 ccl_device_inline
140 #endif
141 void path_rng_2D(KernelGlobals *kg, ccl_addr_space RNG *rng, int sample, int num_samples, int dimension, float *fx, float *fy)
142 {
143 #ifdef __CMJ__
144         if(kernel_data.integrator.sampling_pattern == SAMPLING_PATTERN_CMJ) {
145                 /* correlated multi-jittered */
146                 int p = *rng + dimension;
147                 cmj_sample_2D(sample, num_samples, p, fx, fy);
148         }
149         else
150 #endif
151         {
152                 /* sobol */
153                 *fx = path_rng_1D(kg, rng, sample, num_samples, dimension);
154                 *fy = path_rng_1D(kg, rng, sample, num_samples, dimension + 1);
155         }
156 }
157
158 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)
159 {
160 #ifdef __SOBOL_FULL_SCREEN__
161         uint px, py;
162         uint bits = 16; /* limits us to 65536x65536 and 65536 samples */
163         uint size = 1 << bits;
164         uint frame = sample;
165
166         *rng = sobol_lookup(bits, frame, x, y, &px, &py);
167
168         *rng ^= kernel_data.integrator.seed;
169
170         if(sample == 0) {
171                 *fx = 0.5f;
172                 *fy = 0.5f;
173         }
174         else {
175                 *fx = size * (float)px * (1.0f/(float)0xFFFFFFFF) - x;
176                 *fy = size * (float)py * (1.0f/(float)0xFFFFFFFF) - y;
177         }
178 #else
179         *rng = *rng_state;
180
181         *rng ^= kernel_data.integrator.seed;
182
183         if(sample == 0) {
184                 *fx = 0.5f;
185                 *fy = 0.5f;
186         }
187         else {
188                 path_rng_2D(kg, rng, sample, num_samples, PRNG_FILTER_U, fx, fy);
189         }
190 #endif
191 }
192
193 ccl_device void path_rng_end(KernelGlobals *kg, ccl_global uint *rng_state, RNG rng)
194 {
195         /* nothing to do */
196 }
197
198 #else
199
200 /* Linear Congruential Generator */
201
202 ccl_device_inline float path_rng_1D(KernelGlobals *kg, RNG& rng, int sample, int num_samples, int dimension)
203 {
204         /* implicit mod 2^32 */
205         rng = (1103515245*(rng) + 12345);
206         return (float)rng * (1.0f/(float)0xFFFFFFFF);
207 }
208
209 ccl_device_inline void path_rng_2D(KernelGlobals *kg, RNG& rng, int sample, int num_samples, int dimension, float *fx, float *fy)
210 {
211         *fx = path_rng_1D(kg, rng, sample, num_samples, dimension);
212         *fy = path_rng_1D(kg, rng, sample, num_samples, dimension + 1);
213 }
214
215 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)
216 {
217         /* load state */
218         *rng = *rng_state;
219
220         *rng ^= kernel_data.integrator.seed;
221
222         if(sample == 0) {
223                 *fx = 0.5f;
224                 *fy = 0.5f;
225         }
226         else {
227                 path_rng_2D(kg, rng, sample, num_samples, PRNG_FILTER_U, fx, fy);
228         }
229 }
230
231 ccl_device void path_rng_end(KernelGlobals *kg, ccl_global uint *rng_state, RNG rng)
232 {
233         /* store state for next sample */
234         *rng_state = rng;
235 }
236
237 #endif
238
239 /* Linear Congruential Generator */
240
241 ccl_device uint lcg_step_uint(uint *rng)
242 {
243         /* implicit mod 2^32 */
244         *rng = (1103515245*(*rng) + 12345);
245         return *rng;
246 }
247
248 ccl_device float lcg_step_float(uint *rng)
249 {
250         /* implicit mod 2^32 */
251         *rng = (1103515245*(*rng) + 12345);
252         return (float)*rng * (1.0f/(float)0xFFFFFFFF);
253 }
254
255 ccl_device uint lcg_init(uint seed)
256 {
257         uint rng = seed;
258         lcg_step_uint(&rng);
259         return rng;
260 }
261
262 /* Path Tracing Utility Functions
263  *
264  * For each random number in each step of the path we must have a unique
265  * dimension to avoid using the same sequence twice.
266  *
267  * For branches in the path we must be careful not to reuse the same number
268  * in a sequence and offset accordingly. */
269
270 ccl_device_inline float path_state_rng_1D(KernelGlobals *kg, ccl_addr_space RNG *rng, const ccl_addr_space PathState *state, int dimension)
271 {
272         return path_rng_1D(kg, rng, state->sample, state->num_samples, state->rng_offset + dimension);
273 }
274
275 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)
276 {
277         /* the rng_offset is not increased for transparent bounces. if we do then
278          * fully transparent objects can become subtly visible by the different
279          * sampling patterns used where the transparent object is.
280          *
281          * however for some random numbers that will determine if we next bounce
282          * is transparent we do need to increase the offset to avoid always making
283          * the same decision */
284         int rng_offset = state->rng_offset + state->transparent_bounce*PRNG_BOUNCE_NUM;
285         return path_rng_1D(kg, rng, state->sample, state->num_samples, rng_offset + dimension);
286 }
287
288 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)
289 {
290         path_rng_2D(kg, rng, state->sample, state->num_samples, state->rng_offset + dimension, fx, fy);
291 }
292
293 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)
294 {
295         return path_rng_1D(kg, rng, state->sample*num_branches + branch, state->num_samples*num_branches, state->rng_offset + dimension);
296 }
297
298 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)
299 {
300         int rng_offset = state->rng_offset + state->transparent_bounce*PRNG_BOUNCE_NUM;
301         return path_rng_1D(kg, rng, state->sample*num_branches + branch, state->num_samples*num_branches, rng_offset + dimension);
302 }
303
304 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)
305 {
306         path_rng_2D(kg, rng, state->sample*num_branches + branch, state->num_samples*num_branches, state->rng_offset + dimension, fx, fy);
307 }
308
309 ccl_device_inline void path_state_branch(PathState *state, int branch, int num_branches)
310 {
311         /* path is splitting into a branch, adjust so that each branch
312          * still gets a unique sample from the same sequence */
313         state->rng_offset += PRNG_BOUNCE_NUM;
314         state->sample = state->sample*num_branches + branch;
315         state->num_samples = state->num_samples*num_branches;
316 }
317
318 ccl_device_inline uint lcg_state_init(RNG *rng, const PathState *state, uint scramble)
319 {
320         return lcg_init(*rng + state->rng_offset + state->sample*scramble);
321 }
322
323 /* TODO(sergey): For until we can use generic address space from OpenCL 2.0. */
324
325 ccl_device_inline uint lcg_state_init_addrspace(ccl_addr_space RNG *rng,
326                                                 const ccl_addr_space PathState *state,
327                                                 uint scramble)
328 {
329         return lcg_init(*rng + state->rng_offset + state->sample*scramble);
330 }
331
332 ccl_device float lcg_step_float_addrspace(ccl_addr_space uint *rng)
333 {
334         /* implicit mod 2^32 */
335         *rng = (1103515245*(*rng) + 12345);
336         return (float)*rng * (1.0f/(float)0xFFFFFFFF);
337 }
338
339 CCL_NAMESPACE_END
340