1+ #include < stdlib.h>
2+ #include < fstream>
3+ #include < vector>
4+
5+ #include " glm\glm.hpp"
6+ using glm::vec3;
7+ using namespace glm ;
8+
9+ using namespace std ;
10+
11+ #include " CIE.h"
12+
13+ const int IMAGE_WIDTH = 400 ;
14+ const int IMAGE_HEIGHT = 400 ;
15+
16+ const int NUM_SAMPLES = 50 ;
17+
18+ const float MAX_DIST = 8000 ;
19+ const float PI = 3.14159 ;
20+
21+ vec3 buffer[IMAGE_WIDTH][IMAGE_HEIGHT];
22+
23+ float nrand ()
24+ {
25+ return (float )rand () / RAND_MAX;
26+ }
27+
28+ struct material
29+ {
30+ float mean = 500 ;
31+ float std = 100 ;
32+
33+ material (float m, float s)
34+ : mean(m), std(s)
35+ {}
36+
37+ material ()
38+ : material(0 , 0 )
39+ {}
40+ };
41+
42+ struct light_path_node
43+ {
44+ vec3 normal ;
45+ vec3 ray;
46+ material m;
47+ float weight;
48+
49+ light_path_node () {}
50+ };
51+
52+ // non-normalized bell curve
53+ float bell (float x, float m, float s)
54+ {
55+ return exp (-(x - m)*(x - m) / (2 * s*s));
56+ }
57+
58+ // -----------------------------------------------------------------------------
59+ // Intersection stuff
60+ // -----------------------------------------------------------------------------
61+ float rPlane (vec3 p, vec3 n, vec3 ray)
62+ {
63+ if (abs (dot (ray, n)) < 0.001 )
64+ return MAX_DIST;
65+ float d = dot (p, n) / dot (ray, n);
66+ if (d < 0 )
67+ return MAX_DIST;
68+ return d;
69+ }
70+ float rSphere (vec3 p, float rad, vec3 ray)
71+ {
72+ if (dot (p, ray) < 0 )
73+ return MAX_DIST;
74+ vec3 q = dot (p, ray) * ray;
75+ float l = length (p - q);
76+ if (l > rad)
77+ return MAX_DIST;
78+
79+ float m = length (q);
80+ float ret = m - sqrt (rad*rad - l*l);
81+ if (ret < 0 )
82+ return MAX_DIST;
83+ return ret;
84+ }
85+
86+ float BRDF (float lambda, vec3 p, vec3 inDir, vec3 outDir)
87+ {
88+ float ret;
89+ if (abs (p.x + 5 ) < 0.1 )
90+ ret = 0.8 * bell (lambda, 600 , 20 );
91+ else if (abs (p.x - 5 ) < 0.1 )
92+ ret = 0.8 * bell (lambda, 550 , 20 );
93+ else
94+ ret = 0.8 ;
95+
96+ return ret / PI;
97+ }
98+ float emmision (vec3 p)
99+ {
100+ if (abs (p.y - 5 ) < 0.1 && length (vec2 (p.x , p.z + 15 )) < 1.414 )
101+ {
102+ return 50 .f ;
103+ }
104+ return 0 .f ;
105+ }
106+
107+ float intersect (vec3 o, vec3 r, vec3& normal )
108+ {
109+ float ret = MAX_DIST;
110+ float q;
111+
112+ q = rPlane (vec3 (0 , -5 , 0 ) - o, vec3 (0 , 1 , 0 ), r);
113+ if (q < ret)
114+ {
115+ ret = q;
116+ normal = vec3 (0 , 1 , 0 );
117+ }
118+ q = rPlane (vec3 (0 , 5 , 0 ) - o, vec3 (0 , -1 , 0 ), r);
119+ if (q < ret)
120+ {
121+ ret = q;
122+ normal = vec3 (0 , -1 , 0 );
123+ }
124+
125+
126+ q = rPlane (vec3 (5 , 0 , 0 ) - o, vec3 (-1 , 0 , 0 ), r);
127+ if (q < ret)
128+ {
129+ ret = q;
130+ normal = vec3 (-1 , 0 , 0 );
131+ }
132+ q = rPlane (vec3 (-5 , 0 , 0 ) - o, vec3 (1 , 0 , 0 ), r);
133+ if (q < ret)
134+ {
135+ ret = q;
136+ normal = vec3 (1 , 0 , 0 );
137+ }
138+
139+
140+ q = rPlane (vec3 (0 , 0 , -20 ) - o, vec3 (0 , 0 , 1 ), r);
141+ if (q < ret)
142+ {
143+ ret = q;
144+ normal = vec3 (0 , 0 , 1 );
145+ }
146+ q = rPlane (vec3 (0 , 0 , 1 ) - o, vec3 (0 , 0 , -1 ), r);
147+ if (q < ret)
148+ {
149+ ret = q;
150+ normal = vec3 (0 , 0 , -1 );
151+ }
152+
153+
154+ q = rSphere (vec3 (-2 , -3 , -17 ) - o, 2 , r);
155+ if (q < ret)
156+ {
157+ ret = q;
158+ normal = normalize (r * q + o - vec3 (-2 , -3 , -17 ));
159+ }
160+
161+ return ret;
162+ }
163+
164+ // -----------------------------------------------------------------------------
165+ // Coloring and ray tracing stuff
166+ // -----------------------------------------------------------------------------
167+ vec3 getTangent (vec3 norm)
168+ {
169+ vec3 tangent;
170+ vec3 c1 = cross (norm, vec3 (0 , 0 , 1 ));
171+ vec3 c2 = cross (norm, vec3 (0 , 1 , 0 ));
172+ if (dot (c1, c1) > dot (c2, c2))
173+ tangent = c1;
174+ else
175+ tangent = c2;
176+ return tangent;
177+ }
178+ vec3 randCosineWeightedRay (vec3 norm)
179+ {
180+ float rx = 1 , rz = 1 ;
181+ while (rx*rx + rz*rz >= 1 )
182+ {
183+ rx = 2 * nrand () - 1 .0f ;
184+ rz = 2 * nrand () - 1 .0f ;
185+ }
186+ float ry = sqrt (1 - rx*rx - rz*rz);
187+
188+ vec3 tangent = getTangent (norm);
189+ vec3 bitangent = cross (norm, tangent);
190+
191+ vec3 castRay = normalize (tangent*rx + bitangent*rz + norm*ry);
192+ return castRay;
193+ }
194+
195+ vec3 radiance (vec3 o, vec3 ray)
196+ {
197+ // Decide a-priori what the wavelength of this sample will be
198+ // Note: this only works if all lights output the exact same spectrums
199+ float lambda = nrand () * 225 + 425 ;
200+
201+ float accum = 1 ;
202+ while (true )
203+ {
204+ // Russian Roulette
205+ float r = nrand ();
206+ float russian = min (1 .f , accum);
207+ if (russian < r)
208+ break ;
209+
210+ vec3 normal ;
211+ float dist = intersect (o, ray, normal ) - 0.001 ;
212+ vec3 p = o + ray * dist;
213+
214+ if (dist > 100 )
215+ return vec3 (0 , 0 , 0 );
216+
217+ vec3 emm = emmision (p) * xyz_to_rgb (wavelength_to_xyz (1 .f , lambda));
218+ if (dot (emm, emm) > 0.001 )
219+ return accum * emm / russian;
220+
221+ vec3 nextRay = randCosineWeightedRay (normal );
222+ float pdf = dot (nextRay, normal ) / PI;
223+
224+ accum *= BRDF (lambda, p, nextRay, -ray) * max (0 .f , dot (normal , nextRay)) / (pdf * russian);
225+
226+ o = p;
227+ ray = nextRay;
228+ }
229+
230+ return vec3 ();
231+ }
232+
233+ int main ()
234+ {
235+ for (int i = 0 ; i < NUM_SAMPLES; ++i)
236+ {
237+ printf (" Iteration %d\n "
238+ for (int x = 0 ; x < IMAGE_WIDTH; ++x)
239+ {
240+ for (int y = 0 ; y < IMAGE_HEIGHT; ++y)
241+ {
242+ vec3 ray = vec3 ((float )(x - IMAGE_WIDTH / 2 ) / IMAGE_WIDTH, (float )(y - IMAGE_HEIGHT / 2 ) / IMAGE_HEIGHT, -1.0 );
243+ ray = normalize (ray);
244+
245+ // buffer[x][y] = vec3((float)x / IMAGE_WIDTH, (float)y / IMAGE_HEIGHT, 0);
246+ buffer[x][y] += radiance (vec3 (0 , -1 , 0 ), ray) / (float )NUM_SAMPLES;
247+ }
248+ }
249+ }
250+
251+ // World's worst tonemapping
252+ for (int x = 0 ; x < IMAGE_WIDTH; ++x)
253+ {
254+ for (int y = 0 ; y < IMAGE_HEIGHT; ++y)
255+ {
256+ buffer[x][y] /= buffer[x][y] + vec3 (1 , 1 , 1 );
257+ }
258+ }
259+
260+ ofstream file (" image.ppm"
261+ file << " P3 " " " " 255"
262+
263+ for (int y = IMAGE_HEIGHT - 1 ; y >= 0 ; --y)
264+ {
265+ for (int x = 0 ; x < IMAGE_WIDTH; ++x)
266+ {
267+ if ((int )(buffer[x][y].x * 255 ) == 0x80000000 ||
268+ (int )(buffer[x][y].y * 255 ) == 0x80000000 ||
269+ (int )(buffer[x][y].z * 255 ) == 0x80000000 )
270+ {
271+ // cout << "a divide by zero happened" << endl;
272+ /* for (int i = 0; i < 50; ++i)
273+ file << "\n";*/
274+ file << 0 << " " 0 << " " 0 << " "
275+ }
276+ else
277+ file << (int )(buffer[x][y].x * 255 ) << " " int )(buffer[x][y].y * 255 ) << " " int )(buffer[x][y].z * 255 ) << " "
278+ }
279+ }
280+
281+ file.close ();
282+
283+ printf (" Finished\n "
284+ // getchar();
285+ }
0 commit comments