1 //////////////////////////////////////////////////////////////////////
2 //
3 // Pixie
4 //
5 // Copyright � 1999 - 2003, Okan Arikan
6 //
7 // Contact: okan@cs.utexas.edu
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the GNU Lesser General Public
11 // License as published by the Free Software Foundation; either
12 // version 2.1 of the License, or (at your option) any later version.
13 //
14 // This library is distributed in the hope that it will be useful,
15 // but WITHOUT ANY WARRANTY; without even the implied warranty of
16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 // Lesser General Public License for more details.
18 //
19 // You should have received a copy of the GNU Lesser General Public
20 // License along with this library; if not, write to the Free Software
21 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22 //
23 ///////////////////////////////////////////////////////////////////////
24 ///////////////////////////////////////////////////////////////////////
25 //
26 // File : random.cpp
27 // Classes : CSobol, CSphereSampler, CCosineSampler
28 // Description : Several random generators
29 //
30 ////////////////////////////////////////////////////////////////////////
31 #include <math.h>
32
33 #include "random.h"
34
35
36 // primitive polynomials in binary encoding
37 const int primitive_polynomials[SOBOL_MAX_DIMENSION] = {
38 1, 3, 7, 11, 13, 19, 25, 37, 59, 47,
39 61, 55, 41, 67, 97, 91, 109, 103, 115, 131,
40 193, 137, 145, 143, 241, 157, 185, 167, 229, 171,
41 213, 191, 253, 203, 211, 239, 247, 285, 369, 299
42 };
43
44 // degrees of the primitive polynomials
45 const int degree_table[SOBOL_MAX_DIMENSION] = {
46 0, 1, 2, 3, 3, 4, 4, 5, 5, 5,
47 5, 5, 5, 6, 6, 6, 6, 6, 6, 7,
48 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
49 7, 7, 7, 7, 7, 7, 7, 8, 8, 8
50 };
51
52
53 // initial values for direction tables, following
54 // Bratley+Fox, taken from [Sobol+Levitan, preprint 1976]
55 const int v_init[8][SOBOL_MAX_DIMENSION] = {
56 {
57 0, 1, 1, 1, 1, 1, 1, 1, 1, 1,
58 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
59 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
60 1, 1, 1, 1, 1, 1, 1, 1, 1, 1
61 },
62 {
63 0, 0, 1, 3, 1, 3, 1, 3, 3, 1,
64 3, 1, 3, 1, 3, 1, 1, 3, 1, 3,
65 1, 3, 1, 3, 3, 1, 3, 1, 3, 1,
66 3, 1, 1, 3, 1, 3, 1, 3, 1, 3
67 },
68 {
69 0, 0, 0, 7, 5, 1, 3, 3, 7, 5,
70 5, 7, 7, 1, 3, 3, 7, 5, 1, 1,
71 5, 3, 3, 1, 7, 5, 1, 3, 3, 7,
72 5, 1, 1, 5, 7, 7, 5, 1, 3, 3
73 },
74 {
75 0, 0, 0, 0, 0, 1, 7, 9, 13, 11,
76 1, 3, 7, 9, 5, 13, 13, 11, 3, 15,
77 5, 3, 15, 7, 9, 13, 9, 1, 11, 7,
78 5, 15, 1, 15, 11, 5, 3, 1, 7, 9
79 },
80 {
81 0, 0, 0, 0, 0, 0, 0, 9, 3, 27,
82 15, 29, 21, 23, 19, 11, 25, 7, 13, 17,
83 1, 25, 29, 3, 31, 11, 5, 23, 27, 19,
84 21, 5, 1, 17, 13, 7, 15, 9, 31, 9
85 },
86 {
87 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
88 0, 0, 0, 37, 33, 7, 5, 11, 39, 63,
89 27, 17, 15, 23, 29, 3, 21, 13, 31, 25,
90 9, 49, 33, 19, 29, 11, 19, 27, 15, 25
91 },
92 {
93 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
94 0, 0, 0, 0, 0, 0, 0, 0, 0, 13,
95 33, 115, 41, 79, 17, 29, 119, 75, 73, 105,
96 7, 59, 65, 21, 3, 113, 61, 89, 45, 107
97 },
98 {
99 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
100 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
101 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
102 0, 0, 0, 0, 0, 0, 0, 7, 23, 39
103 }
104 };
105
106
107
108
109
110
111
112
113
114
115
116
117
118 ///////////////////////////////////////////////////////////////////////
119 // Function : sampleHemisphere
120 // Description : Sample vectors distributed uniformly in a hemisphere
121 // Return Value :
122 // Comments : Z must be unit
sampleHemisphere(float * R,const float * Z,const float theta,CSobol<4> & generator)123 void sampleHemisphere(float *R,const float *Z,const float theta,CSobol<4> &generator) {
124 float P[4];
125 vector Po;
126 float cosa;
127 float sina;
128
129 while(TRUE) {
130 generator.get(P);
131
132 // Sample a uniformly distributed point on a sphere
133 P[COMP_X] = 2*P[COMP_X]-1;
134 P[COMP_Y] = 2*P[COMP_Y]-1;
135 P[COMP_Z] = 2*P[COMP_Z]-1;
136
137 // did we get something inside the unit sphere and non-zero
138 const float l = dotvv(P,P);
139 if (l < 1 && l > C_EPSILON) break;
140 }
141
142 cosa = 1 - P[3]*(1 - (float) cos(theta));
143 sina = sqrtf(1 - cosa*cosa);
144
145 // Po is orthagonal to N
146 crossvv(Po,P,Z);
147 // Po is unit length
148 normalizev(Po);
149 // Construct the sample vector
150 mulvf(R,Z,cosa);
151 mulvf(Po,sina);
152 addvv(R,Po);
153 }
154
155
156 ///////////////////////////////////////////////////////////////////////
157 // Function : sampleHemisphere
158 // Description : Sample vectors distributed uniformly in a hemisphere
159 // Return Value :
160 // Comments : Z must be unit
sampleCosineHemisphere(float * R,const float * Z,const float theta,CSobol<4> & generator)161 void sampleCosineHemisphere(float *R,const float *Z,const float theta,CSobol<4> &generator) {
162 float P[4];
163 vector Po;
164 float cosa;
165 float sina;
166 const float cosmin = (float) cos(theta);
167
168 while(TRUE) {
169 generator.get(P);
170
171 // Sample a uniformly distributed point on a sphere
172 P[COMP_X] = 2*P[COMP_X]-1;
173 P[COMP_Y] = 2*P[COMP_Y]-1;
174 P[COMP_Z] = 2*P[COMP_Z]-1;
175
176 // did we get something inside the unit sphere and non-zero
177 const float l = dotvv(P,P);
178 if (l < 1 && l > C_EPSILON) break;
179 }
180
181 cosa = sqrtf(P[3])*(1 - cosmin) + cosmin;
182 sina = sqrtf(1 - cosa*cosa);
183
184 // Po is orthagonal to N
185 crossvv(Po,P,Z);
186 // Po is unit length
187 normalizev(Po);
188 // Construct the sample vector
189 mulvf(R,Z,cosa);
190 mulvf(Po,sina);
191 addvv(R,Po);
192 }
193
194
195 ///////////////////////////////////////////////////////////////////////
196 // Function : sampleSphere
197 // Description : Sample a point in a unit sphere
198 // Return Value :
199 // Comments :
sampleSphere(float * P,CSobol<3> & generator)200 void sampleSphere(float *P,CSobol<3> &generator) {
201 float r[3];
202
203 while(TRUE) {
204
205 generator.get(r);
206
207 // Sample a uniformly distributed point on a sphere
208 P[COMP_X] = 2*r[0]-1;
209 P[COMP_Y] = 2*r[1]-1;
210 P[COMP_Z] = 2*r[2]-1;
211
212 if (dotvv(P,P) < 1) break;
213 }
214 }
215
216