Cycles: svn merge -r41225:41232 ^/trunk/blender
[blender.git] / intern / cycles / kernel / kernel_random.h
1 /*
2  * Copyright 2011, Blender Foundation.
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
18
19 CCL_NAMESPACE_BEGIN
20
21 typedef uint RNG;
22
23 #ifdef __SOBOL__
24
25 /* High Dimensional Sobol */
26
27 /* van der corput radical inverse */
28 __device uint van_der_corput(uint bits)
29 {
30         bits = (bits << 16) | (bits >> 16);
31         bits = ((bits & 0x00ff00ff) << 8) | ((bits & 0xff00ff00) >> 8);
32         bits = ((bits & 0x0f0f0f0f) << 4) | ((bits & 0xf0f0f0f0) >> 4);
33         bits = ((bits & 0x33333333) << 2) | ((bits & 0xcccccccc) >> 2);
34         bits = ((bits & 0x55555555) << 1) | ((bits & 0xaaaaaaaa) >> 1);
35         return bits;
36 }
37
38 /* sobol radical inverse */
39 __device uint sobol(uint i)
40 {
41         uint r = 0;
42
43         for(uint v = 1U << 31; i; i >>= 1, v ^= v >> 1)
44                 if(i & 1)
45                         r ^= v;
46
47         return r;
48 }
49
50 /* inverse of sobol radical inverse */
51 __device uint sobol_inverse(uint i)
52 {
53         const uint msb = 1U << 31;
54         uint r = 0;
55
56         for(uint v = 1; i; i <<= 1, v ^= v << 1)
57                 if(i & msb)
58                         r ^= v;
59
60         return r;
61 }
62
63 /* multidimensional sobol with generator matrices
64    dimension 0 and 1 are equal to van_der_corput() and sobol() respectively */
65 __device uint sobol_dimension(KernelGlobals *kg, int index, int dimension)
66 {
67         uint result = 0;
68         uint i = index;
69
70         for(uint j = 0; i; i >>= 1, j++)
71                 if(i & 1)
72                         result ^= kernel_tex_fetch(__sobol_directions, 32*dimension + j);
73         
74         return result;
75 }
76
77 /* lookup index and x/y coordinate, assumes m is a power of two */
78 __device uint sobol_lookup(const uint m, const uint frame, const uint ex, const uint ey, uint *x, uint *y)
79 {
80         /* shift is constant per frame */
81         const uint shift = frame << (m << 1);
82         const uint sobol_shift = sobol(shift);
83         /* van der Corput is its own inverse */
84         const uint lower = van_der_corput(ex << (32 - m));
85         /* need to compensate for ey difference and shift */
86         const uint sobol_lower = sobol(lower);
87         const uint mask = ~-(1 << m) << (32 - m); /* only m upper bits */
88         const uint delta = ((ey << (32 - m)) ^ sobol_lower ^ sobol_shift) & mask;
89         /* only use m upper bits for the index (m is a power of two) */
90         const uint sobol_result = delta | (delta >> m);
91         const uint upper = sobol_inverse(sobol_result);
92         const uint index = shift | upper | lower;
93         *x = van_der_corput(index);
94         *y = sobol_shift ^ sobol_result ^ sobol_lower;
95         return index;
96 }
97
98 __device_inline float path_rng(KernelGlobals *kg, RNG *rng, int sample, int dimension)
99 {
100 #ifdef __SOBOL_FULL_SCREEN__
101         uint result = sobol_dimension(kg, *rng, dimension);
102         float r = (float)result * (1.0f/(float)0xFFFFFFFF);
103         return r;
104 #else
105         /* compute sobol sequence value using direction vectors */
106         uint result = sobol_dimension(kg, sample, dimension);
107         float r = (float)result * (1.0f/(float)0xFFFFFFFF);
108
109         /* Cranly-Patterson rotation using rng seed */
110         float shift;
111
112         if(dimension & 1)
113                 shift = (*rng >> 16)/((float)0xFFFF);
114         else
115                 shift = (*rng & 0xFFFF)/((float)0xFFFF);
116
117         return r + shift - floor(r + shift);
118 #endif
119 }
120
121 __device_inline void path_rng_init(KernelGlobals *kg, __global uint *rng_state, int sample, RNG *rng, int x, int y, float *fx, float *fy)
122 {
123 #ifdef __SOBOL_FULL_SCREEN__
124         uint px, py;
125         uint bits = 16; /* limits us to 65536x65536 and 65536 samples */
126         uint size = 1 << bits;
127         uint frame = sample;
128
129         *rng = sobol_lookup(bits, frame, x, y, &px, &py);
130
131         *fx = size * (float)px * (1.0f/(float)0xFFFFFFFF) - x;
132         *fy = size * (float)py * (1.0f/(float)0xFFFFFFFF) - y;
133 #else
134         *rng = rng_state[x + y*kernel_data.cam.width];
135
136         *fx = path_rng(kg, rng, sample, PRNG_FILTER_U);
137         *fy = path_rng(kg, rng, sample, PRNG_FILTER_V);
138 #endif
139 }
140
141 __device void path_rng_end(KernelGlobals *kg, __global uint *rng_state, RNG rng, int x, int y)
142 {
143         /* nothing to do */
144 }
145
146 #else
147
148 /* Linear Congruential Generator */
149
150 __device float path_rng(KernelGlobals *kg, RNG *rng, int sample, int dimension)
151 {
152         /* implicit mod 2^32 */
153         *rng = (1103515245*(*rng) + 12345);
154         return (float)*rng * (1.0f/(float)0xFFFFFFFF);
155 }
156
157 __device void path_rng_init(KernelGlobals *kg, __global uint *rng_state, int sample, RNG *rng, int x, int y, float *fx, float *fy)
158 {
159         /* load state */
160         *rng = rng_state[x + y*kernel_data.cam.width];
161
162         *fx = path_rng(kg, rng, sample, PRNG_FILTER_U);
163         *fy = path_rng(kg, rng, sample, PRNG_FILTER_V);
164 }
165
166 __device void path_rng_end(KernelGlobals *kg, __global uint *rng_state, RNG rng, int x, int y)
167 {
168         /* store state for next sample */
169         rng_state[x + y*kernel_data.cam.width] = rng;
170 }
171
172 #endif
173
174 CCL_NAMESPACE_END
175