fc33e226051e12a842685b1399a72dc096e8b062
[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 /* skip initial numbers that are not as well distributed, especially the
26  * first sequence is just 0 everywhere, which can be problematic for e.g.
27  * path termination */
28 #define SOBOL_SKIP 64
29
30 /* High Dimensional Sobol */
31
32 /* van der corput radical inverse */
33 __device uint van_der_corput(uint bits)
34 {
35         bits = (bits << 16) | (bits >> 16);
36         bits = ((bits & 0x00ff00ff) << 8) | ((bits & 0xff00ff00) >> 8);
37         bits = ((bits & 0x0f0f0f0f) << 4) | ((bits & 0xf0f0f0f0) >> 4);
38         bits = ((bits & 0x33333333) << 2) | ((bits & 0xcccccccc) >> 2);
39         bits = ((bits & 0x55555555) << 1) | ((bits & 0xaaaaaaaa) >> 1);
40         return bits;
41 }
42
43 /* sobol radical inverse */
44 __device uint sobol(uint i)
45 {
46         uint r = 0;
47
48         for(uint v = 1U << 31; i; i >>= 1, v ^= v >> 1)
49                 if(i & 1)
50                         r ^= v;
51
52         return r;
53 }
54
55 /* inverse of sobol radical inverse */
56 __device uint sobol_inverse(uint i)
57 {
58         const uint msb = 1U << 31;
59         uint r = 0;
60
61         for(uint v = 1; i; i <<= 1, v ^= v << 1)
62                 if(i & msb)
63                         r ^= v;
64
65         return r;
66 }
67
68 /* multidimensional sobol with generator matrices
69  * dimension 0 and 1 are equal to van_der_corput() and sobol() respectively */
70 __device uint sobol_dimension(KernelGlobals *kg, int index, int dimension)
71 {
72         uint result = 0;
73         uint i = index;
74
75         for(uint j = 0; i; i >>= 1, j++)
76                 if(i & 1)
77                         result ^= kernel_tex_fetch(__sobol_directions, 32*dimension + j);
78         
79         return result;
80 }
81
82 /* lookup index and x/y coordinate, assumes m is a power of two */
83 __device uint sobol_lookup(const uint m, const uint frame, const uint ex, const uint ey, uint *x, uint *y)
84 {
85         /* shift is constant per frame */
86         const uint shift = frame << (m << 1);
87         const uint sobol_shift = sobol(shift);
88         /* van der Corput is its own inverse */
89         const uint lower = van_der_corput(ex << (32 - m));
90         /* need to compensate for ey difference and shift */
91         const uint sobol_lower = sobol(lower);
92         const uint mask = ~-(1 << m) << (32 - m); /* only m upper bits */
93         const uint delta = ((ey << (32 - m)) ^ sobol_lower ^ sobol_shift) & mask;
94         /* only use m upper bits for the index (m is a power of two) */
95         const uint sobol_result = delta | (delta >> m);
96         const uint upper = sobol_inverse(sobol_result);
97         const uint index = shift | upper | lower;
98         *x = van_der_corput(index);
99         *y = sobol_shift ^ sobol_result ^ sobol_lower;
100         return index;
101 }
102
103 __device_inline float path_rng(KernelGlobals *kg, RNG *rng, int sample, int dimension)
104 {
105 #ifdef __SOBOL_FULL_SCREEN__
106         uint result = sobol_dimension(kg, *rng, dimension);
107         float r = (float)result * (1.0f/(float)0xFFFFFFFF);
108         return r;
109 #else
110         /* compute sobol sequence value using direction vectors */
111         uint result = sobol_dimension(kg, sample + SOBOL_SKIP, dimension);
112         float r = (float)result * (1.0f/(float)0xFFFFFFFF);
113
114         /* Cranly-Patterson rotation using rng seed */
115         float shift;
116
117         if(dimension & 1)
118                 shift = (*rng >> 16)/((float)0xFFFF);
119         else
120                 shift = (*rng & 0xFFFF)/((float)0xFFFF);
121
122         return r + shift - floorf(r + shift);
123 #endif
124 }
125
126 __device_inline void path_rng_init(KernelGlobals *kg, __global uint *rng_state, int sample, RNG *rng, int x, int y, float *fx, float *fy)
127 {
128 #ifdef __SOBOL_FULL_SCREEN__
129         uint px, py;
130         uint bits = 16; /* limits us to 65536x65536 and 65536 samples */
131         uint size = 1 << bits;
132         uint frame = sample;
133
134         *rng = sobol_lookup(bits, frame, x, y, &px, &py);
135
136         *rng ^= kernel_data.integrator.seed;
137
138         if(sample == 0) {
139                 *fx = 0.5f;
140                 *fy = 0.5f;
141         }
142         else {
143                 *fx = size * (float)px * (1.0f/(float)0xFFFFFFFF) - x;
144                 *fy = size * (float)py * (1.0f/(float)0xFFFFFFFF) - y;
145         }
146 #else
147         *rng = *rng_state;
148
149         *rng ^= kernel_data.integrator.seed;
150
151         if(sample == 0) {
152                 *fx = 0.5f;
153                 *fy = 0.5f;
154         }
155         else {
156                 *fx = path_rng(kg, rng, sample, PRNG_FILTER_U);
157                 *fy = path_rng(kg, rng, sample, PRNG_FILTER_V);
158         }
159 #endif
160 }
161
162 __device void path_rng_end(KernelGlobals *kg, __global uint *rng_state, RNG rng)
163 {
164         /* nothing to do */
165 }
166
167 #else
168
169 /* Linear Congruential Generator */
170
171 __device float path_rng(KernelGlobals *kg, RNG *rng, int sample, int dimension)
172 {
173         /* implicit mod 2^32 */
174         *rng = (1103515245*(*rng) + 12345);
175         return (float)*rng * (1.0f/(float)0xFFFFFFFF);
176 }
177
178 __device void path_rng_init(KernelGlobals *kg, __global uint *rng_state, int sample, RNG *rng, int x, int y, float *fx, float *fy)
179 {
180         /* load state */
181         *rng = *rng_state;
182
183         *rng ^= kernel_data.integrator.seed;
184
185         if(sample == 0) {
186                 *fx = 0.5f;
187                 *fy = 0.5f;
188         }
189         else {
190                 *fx = path_rng(kg, rng, sample, PRNG_FILTER_U);
191                 *fy = path_rng(kg, rng, sample, PRNG_FILTER_V);
192         }
193 }
194
195 __device void path_rng_end(KernelGlobals *kg, __global uint *rng_state, RNG rng)
196 {
197         /* store state for next sample */
198         *rng_state = rng;
199 }
200
201 #endif
202
203 __device float lcg_step(uint *rng)
204 {
205         /* implicit mod 2^32 */
206         *rng = (1103515245*(*rng) + 12345);
207         return (float)*rng * (1.0f/(float)0xFFFFFFFF);
208 }
209
210 __device uint lcg_init(uint seed)
211 {
212         uint rng = seed;
213         lcg_step(&rng);
214         return rng;
215 }
216
217 CCL_NAMESPACE_END
218