|
1 /* |
|
2 * shapes.cc: shape classes |
|
3 * |
|
4 * This file is part of Pyrit Ray Tracer. |
|
5 * |
|
6 * Copyright 2006, 2007, 2008 Radek Brich |
|
7 * |
|
8 * Permission is hereby granted, free of charge, to any person obtaining a copy |
|
9 * of this software and associated documentation files (the "Software"), to deal |
|
10 * in the Software without restriction, including without limitation the rights |
|
11 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell |
|
12 * copies of the Software, and to permit persons to whom the Software is |
|
13 * furnished to do so, subject to the following conditions: |
|
14 * |
|
15 * The above copyright notice and this permission notice shall be included in |
|
16 * all copies or substantial portions of the Software. |
|
17 * |
|
18 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR |
|
19 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, |
|
20 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE |
|
21 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER |
|
22 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, |
|
23 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN |
|
24 * THE SOFTWARE. |
|
25 */ |
|
26 |
|
27 #include "shapes.h" |
|
28 #include "serialize.h" |
|
29 |
|
30 bool Sphere::intersect(const Ray &ray, Float &dist) const |
|
31 { |
|
32 Vector3 V = ray.o - center; |
|
33 register Float d = -dot(V, ray.dir); |
|
34 register Float Det = d * d - (dot(V,V) - sqr_radius); |
|
35 register Float t1,t2; |
|
36 if (Det > 0) { |
|
37 Det = sqrtf(Det); |
|
38 t1 = d - Det; |
|
39 t2 = d + Det; |
|
40 if (t1 > 0) |
|
41 { |
|
42 if (t1 < dist) |
|
43 { |
|
44 dist = t1; |
|
45 return true; |
|
46 } |
|
47 } |
|
48 else if (t2 > 0 && t2 < dist) |
|
49 { |
|
50 dist = t2; |
|
51 return true; |
|
52 } |
|
53 } |
|
54 return false; |
|
55 } |
|
56 |
|
57 /* if there should be CSG sometimes, this may be needed... */ |
|
58 bool Sphere::intersect_all(const Ray &ray, Float dist, vector<Float> &allts) const |
|
59 { |
|
60 //allts = new vector<Float>(); |
|
61 |
|
62 Vector3 V = ((Ray)ray).o - center; |
|
63 Float Vd = - dot(V, ray.dir); |
|
64 Float Det = Vd * Vd - (dot(V,V) - sqr_radius); |
|
65 |
|
66 if (Det > 0) { |
|
67 Det = sqrtf(Det); |
|
68 Float t1 = Vd - Det; |
|
69 Float t2 = Vd + Det; |
|
70 if (t1 < 0) |
|
71 { |
|
72 if (t2 > 0) |
|
73 { |
|
74 // ray from inside of the sphere |
|
75 allts.push_back(0.0); |
|
76 allts.push_back(t2); |
|
77 return true; |
|
78 } |
|
79 else |
|
80 return false; |
|
81 } |
|
82 else |
|
83 { |
|
84 allts.push_back(t1); |
|
85 allts.push_back(t2); |
|
86 return true; |
|
87 } |
|
88 } |
|
89 return false; |
|
90 } |
|
91 |
|
92 bool Sphere::intersect_bbox(const BBox &bbox) const |
|
93 { |
|
94 register float dmin = 0; |
|
95 for (int i = 0; i < 3; i++) |
|
96 { |
|
97 if (center[i] < bbox.L[i]) |
|
98 dmin += (center[i] - bbox.L[i])*(center[i] - bbox.L[i]); |
|
99 else |
|
100 if (center[i] > bbox.H[i]) |
|
101 dmin += (center[i] - bbox.H[i])*(center[i] - bbox.H[i]); |
|
102 } |
|
103 if (dmin <= sqr_radius) |
|
104 return true; |
|
105 return false; |
|
106 }; |
|
107 |
|
108 BBox Sphere::get_bbox() const |
|
109 { |
|
110 return BBox(center - radius, center + radius); |
|
111 } |
|
112 |
|
113 bool Box::intersect(const Ray &ray, Float &dist) const |
|
114 { |
|
115 register Float tnear = -Inf; |
|
116 register Float tfar = Inf; |
|
117 register Float t1, t2; |
|
118 |
|
119 for (int i = 0; i < 3; i++) |
|
120 { |
|
121 if (ray.dir[i] == 0) { |
|
122 /* ray is parallel to these planes */ |
|
123 if (ray.o[i] < L[i] || ray.o[i] > H[i]) |
|
124 return false; |
|
125 } |
|
126 else |
|
127 { |
|
128 /* compute the intersection distance of the planes */ |
|
129 t1 = (L[i] - ray.o[i]) / ray.dir[i]; |
|
130 t2 = (H[i] - ray.o[i]) / ray.dir[i]; |
|
131 |
|
132 if (t1 > t2) |
|
133 swap(t1, t2); |
|
134 |
|
135 if (t1 > tnear) |
|
136 tnear = t1; /* want largest Tnear */ |
|
137 if (t2 < tfar) |
|
138 tfar = t2; /* want smallest Tfar */ |
|
139 if (tnear > tfar || tfar < 0) |
|
140 return false; /* box missed; box is behind ray */ |
|
141 } |
|
142 } |
|
143 |
|
144 if (tnear < dist) |
|
145 { |
|
146 dist = tnear; |
|
147 return true; |
|
148 } |
|
149 return false; |
|
150 } |
|
151 |
|
152 bool Box::intersect_bbox(const BBox &bbox) const |
|
153 { |
|
154 return ( |
|
155 H.x > bbox.L.x && L.x < bbox.H.x && |
|
156 H.y > bbox.L.y && L.y < bbox.H.y && |
|
157 H.z > bbox.L.z && L.z < bbox.H.z); |
|
158 } |
|
159 |
|
160 const Vector3 Box::normal(const Vector3 &P) const |
|
161 { |
|
162 register Vector3 l = P - L; |
|
163 register Vector3 h = H - P; |
|
164 |
|
165 if (l.x < h.x) |
|
166 h.x = -1; |
|
167 else |
|
168 { |
|
169 l.x = h.x; |
|
170 h.x = +1; |
|
171 } |
|
172 |
|
173 if (l.y < h.y) |
|
174 h.y = -1; |
|
175 else |
|
176 { |
|
177 l.y = h.y; |
|
178 h.y = +1; |
|
179 } |
|
180 |
|
181 if (l.z < h.z) |
|
182 h.z = -1; |
|
183 else |
|
184 { |
|
185 l.z = h.z; |
|
186 h.z = +1; |
|
187 } |
|
188 |
|
189 if (l.x > l.y) |
|
190 { |
|
191 h.x = 0; |
|
192 if (l.y > l.z) |
|
193 h.y = 0; |
|
194 else |
|
195 h.z = 0; |
|
196 } |
|
197 else |
|
198 { |
|
199 h.y = 0; |
|
200 if (l.x > l.z) |
|
201 h.x = 0; |
|
202 else |
|
203 h.z = 0; |
|
204 } |
|
205 return h; |
|
206 } |
|
207 |
|
208 #ifdef TRI_PLUCKER |
|
209 inline void Plucker(const Vector3 &p, const Vector3 &q, Float* pl) |
|
210 { |
|
211 pl[0] = p.x*q.y - q.x*p.y; |
|
212 pl[1] = p.x*q.z - q.x*p.z; |
|
213 pl[2] = p.x - q.x; |
|
214 pl[3] = p.y*q.z - q.y*p.z; |
|
215 pl[4] = p.z - q.z; |
|
216 pl[5] = q.y - p.y; |
|
217 } |
|
218 |
|
219 inline Float Side(const Float* pla, const Float* plb) |
|
220 { |
|
221 return pla[0]*plb[4] + pla[1]*plb[5] + pla[2]*plb[3] + pla[4]*plb[0] + pla[5]*plb[1] + pla[3]*plb[2]; |
|
222 } |
|
223 #endif |
|
224 |
|
225 Triangle::Triangle(Vertex *aA, Vertex *aB, Vertex *aC, Material *amaterial) |
|
226 : A(aA), B(aB), C(aC) |
|
227 { |
|
228 material = amaterial; |
|
229 |
|
230 const Vector3 c = B->P - A->P; |
|
231 const Vector3 b = C->P - A->P; |
|
232 |
|
233 N = cross(c, b); |
|
234 N.normalize(); |
|
235 |
|
236 #ifdef TRI_PLUCKER |
|
237 Plucker(B->P,C->P,pla); |
|
238 Plucker(C->P,A->P,plb); |
|
239 Plucker(A->P,B->P,plc); |
|
240 #endif |
|
241 |
|
242 #if defined(TRI_BARI) || defined(TRI_BARI_PRE) |
|
243 if (fabsf(N.x) > fabsf(N.y)) |
|
244 { |
|
245 if (fabsf(N.x) > fabsf(N.z)) |
|
246 k = 0; |
|
247 else |
|
248 k = 2; |
|
249 } |
|
250 else |
|
251 { |
|
252 if (fabsf(N.y) > fabsf(N.z)) |
|
253 k = 1; |
|
254 else |
|
255 k = 2; |
|
256 } |
|
257 #endif |
|
258 #ifdef TRI_BARI_PRE |
|
259 int u = (k + 1) % 3; |
|
260 int v = (k + 2) % 3; |
|
261 |
|
262 Float krec = 1.0 / N[k]; |
|
263 nu = N[u] * krec; |
|
264 nv = N[v] * krec; |
|
265 nd = dot(N, A->P) * krec; |
|
266 |
|
267 // first line equation |
|
268 Float reci = 1.0f / (b[u] * c[v] - b[v] * c[u]); |
|
269 bnu = b[u] * reci; |
|
270 bnv = -b[v] * reci; |
|
271 |
|
272 // second line equation |
|
273 cnu = -c[u] * reci; |
|
274 cnv = c[v] * reci; |
|
275 #endif |
|
276 } |
|
277 |
|
278 bool Triangle::intersect(const Ray &ray, Float &dist) const |
|
279 { |
|
280 #ifdef TRI_PLUCKER |
|
281 Float plr[6]; |
|
282 Plucker(ray.o, ray.o+ray.dir, plr); |
|
283 const bool side0 = Side(plr, pla) >= 0.0; |
|
284 const bool side1 = Side(plr, plb) >= 0.0; |
|
285 if (side0 != side1) |
|
286 return false; |
|
287 const bool side2 = Side(plr, plc) >= 0.0; |
|
288 if (side0 != side2) |
|
289 return false; |
|
290 const Float t = - dot( (ray.o - A->P), N) / dot(ray.dir,N); |
|
291 if(t <= Eps || t >= dist) |
|
292 return false; |
|
293 dist = t; |
|
294 return true; |
|
295 #endif |
|
296 |
|
297 #if defined(TRI_BARI) || defined(TRI_BARI_PRE) |
|
298 static const int modulo3[5] = {0,1,2,0,1}; |
|
299 const Vector3 &O = ray.o; |
|
300 const Vector3 &D = ray.dir; |
|
301 register const int u = modulo3[k+1]; |
|
302 register const int v = modulo3[k+2]; |
|
303 #endif |
|
304 #ifdef TRI_BARI_PRE |
|
305 const Float t = (nd - O[k] - nu * O[u] - nv * O[v]) / (D[k] + nu * D[u] + nv * D[v]); |
|
306 |
|
307 if (t >= dist || t < Eps) |
|
308 return false; |
|
309 |
|
310 const Float hu = O[u] + t * D[u] - A->P[u]; |
|
311 const Float hv = O[v] + t * D[v] - A->P[v]; |
|
312 const Float beta = hv * bnu + hu * bnv; |
|
313 |
|
314 if (beta < 0.) |
|
315 return false; |
|
316 |
|
317 const Float gamma = hu * cnv + hv * cnu; |
|
318 if (gamma < 0. || beta + gamma > 1.) |
|
319 return false; |
|
320 |
|
321 dist = t; |
|
322 return true; |
|
323 #endif |
|
324 |
|
325 #ifdef TRI_BARI |
|
326 // original barycentric coordinates based intesection |
|
327 // not optimized, just for reference |
|
328 const Vector3 c = B - A; |
|
329 const Vector3 b = C - A; |
|
330 // distance test |
|
331 const Float t = - dot( (O-A), N) / dot(D,N); |
|
332 if (t < Eps || t > dist) |
|
333 return false; |
|
334 |
|
335 // calc hitpoint |
|
336 const Float Hu = O[u] + t * D[u] - A[u]; |
|
337 const Float Hv = O[v] + t * D[v] - A[v]; |
|
338 const Float beta = (b[u] * Hv - b[v] * Hu) / (b[u] * c[v] - b[v] * c[u]); |
|
339 if (beta < 0) |
|
340 return false; |
|
341 const Float gamma = (c[v] * Hu - c[u] * Hv) / (b[u] * c[v] - b[v] * c[u]); |
|
342 if (gamma < 0) |
|
343 return false; |
|
344 if (beta+gamma > 1) |
|
345 return false; |
|
346 dist = t; |
|
347 return true; |
|
348 #endif |
|
349 } |
|
350 |
|
351 bool Triangle::intersect_bbox(const BBox &bbox) const |
|
352 { |
|
353 const Vector3 boxcenter = (bbox.L+bbox.H)*0.5; |
|
354 const Vector3 boxhalfsize = (bbox.H-bbox.L)*0.5; |
|
355 const Vector3 v0 = A->P - boxcenter; |
|
356 const Vector3 v1 = B->P - boxcenter; |
|
357 const Vector3 v2 = C->P - boxcenter; |
|
358 const Vector3 e0 = v1-v0; |
|
359 const Vector3 e1 = v2-v1; |
|
360 const Vector3 e2 = v0-v2; |
|
361 |
|
362 Float fex = fabsf(e0.x); |
|
363 Float fey = fabsf(e0.y); |
|
364 Float fez = fabsf(e0.z); |
|
365 |
|
366 Float p0,p1,p2,min,max,rad; |
|
367 |
|
368 p0 = e0.z*v0.y - e0.y*v0.z; |
|
369 p2 = e0.z*v2.y - e0.y*v2.z; |
|
370 if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;} |
|
371 rad = fez * boxhalfsize.y + fey * boxhalfsize.z; |
|
372 if(min>rad || max<-rad) return false; |
|
373 |
|
374 p0 = -e0.z*v0.x + e0.x*v0.z; |
|
375 p2 = -e0.z*v2.x + e0.x*v2.z; |
|
376 if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;} |
|
377 rad = fez * boxhalfsize.x + fex * boxhalfsize.z; |
|
378 if(min>rad || max<-rad) return false; |
|
379 |
|
380 p1 = e0.y*v1.x - e0.x*v1.y; |
|
381 p2 = e0.y*v2.x - e0.x*v2.y; |
|
382 if(p2<p1) {min=p2; max=p1;} else {min=p1; max=p2;} |
|
383 rad = fey * boxhalfsize.x + fex * boxhalfsize.y; |
|
384 if(min>rad || max<-rad) return false; |
|
385 |
|
386 fex = fabsf(e1.x); |
|
387 fey = fabsf(e1.y); |
|
388 fez = fabsf(e1.z); |
|
389 |
|
390 p0 = e1.z*v0.y - e1.y*v0.z; |
|
391 p2 = e1.z*v2.y - e1.y*v2.z; |
|
392 if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;} |
|
393 rad = fez * boxhalfsize.y + fey * boxhalfsize.z; |
|
394 if(min>rad || max<-rad) return false; |
|
395 |
|
396 p0 = -e1.z*v0.x + e1.x*v0.z; |
|
397 p2 = -e1.z*v2.x + e1.x*v2.z; |
|
398 if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;} |
|
399 rad = fez * boxhalfsize.x + fex * boxhalfsize.z; |
|
400 if(min>rad || max<-rad) return false; |
|
401 |
|
402 p0 = e1.y*v0.x - e1.x*v0.y; |
|
403 p1 = e1.y*v1.x - e1.x*v1.y; |
|
404 if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} |
|
405 rad = fey * boxhalfsize.x + fex * boxhalfsize.y; |
|
406 if(min>rad || max<-rad) return false; |
|
407 |
|
408 fex = fabsf(e2.x); |
|
409 fey = fabsf(e2.y); |
|
410 fez = fabsf(e2.z); |
|
411 |
|
412 p0 = e2.z*v0.y - e2.y*v0.z; |
|
413 p1 = e2.z*v1.y - e2.y*v1.z; |
|
414 if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} |
|
415 rad = fez * boxhalfsize.y + fey * boxhalfsize.z; |
|
416 if(min>rad || max<-rad) return false; |
|
417 |
|
418 p0 = -e2.z*v0.x + e2.x*v0.z; |
|
419 p1 = -e2.z*v1.x + e2.x*v1.z; |
|
420 if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} |
|
421 rad = fez * boxhalfsize.x + fex * boxhalfsize.z; |
|
422 if(min>rad || max<-rad) return false; |
|
423 |
|
424 p1 = e2.y*v1.x - e2.x*v1.y; |
|
425 p2 = e2.y*v2.x - e2.x*v2.y; |
|
426 if(p2<p1) {min=p2; max=p1;} else {min=p1; max=p2;} |
|
427 rad = fey * boxhalfsize.x + fex * boxhalfsize.y; |
|
428 if(min>rad || max<-rad) return false; |
|
429 |
|
430 /* test overlap in the {x,y,z}-directions */ |
|
431 /* test in X-direction */ |
|
432 min = v0.x; |
|
433 if (v1.x < min) min = v1.x; |
|
434 if (v2.x < min) min = v2.x; |
|
435 max = v0.x; |
|
436 if (v1.x > max) max = v1.x; |
|
437 if (v2.x > max) max = v2.x; |
|
438 if(min>boxhalfsize.x || max<-boxhalfsize.x) return false; |
|
439 |
|
440 /* test in Y-direction */ |
|
441 min = v0.y; |
|
442 if (v1.y < min) min = v1.y; |
|
443 if (v2.y < min) min = v2.y; |
|
444 max = v0.y; |
|
445 if (v1.y > max) max = v1.y; |
|
446 if (v2.y > max) max = v2.y; |
|
447 if(min>boxhalfsize.y || max<-boxhalfsize.y) return false; |
|
448 |
|
449 /* test in Z-direction */ |
|
450 min = v0.z; |
|
451 if (v1.z < min) min = v1.z; |
|
452 if (v2.z < min) min = v2.z; |
|
453 max = v0.z; |
|
454 if (v1.z > max) max = v1.z; |
|
455 if (v2.z > max) max = v2.z; |
|
456 if(min>boxhalfsize.z || max<-boxhalfsize.z) return false; |
|
457 |
|
458 /* test if the box intersects the plane of the triangle */ |
|
459 Vector3 vmin,vmax; |
|
460 Float v; |
|
461 for(int q=0;q<3;q++) |
|
462 { |
|
463 v=v0[q]; |
|
464 if(N[q]>0.0f) |
|
465 { |
|
466 vmin.cell[q]=-boxhalfsize[q] - v; |
|
467 vmax.cell[q]= boxhalfsize[q] - v; |
|
468 } |
|
469 else |
|
470 { |
|
471 vmin.cell[q]= boxhalfsize[q] - v; |
|
472 vmax.cell[q]=-boxhalfsize[q] - v; |
|
473 } |
|
474 } |
|
475 if(dot(N,vmin)>0.0f) return false; |
|
476 if(dot(N,vmax)>=0.0f) return true; |
|
477 |
|
478 return false; |
|
479 } |
|
480 |
|
481 BBox Triangle::get_bbox() const |
|
482 { |
|
483 BBox bbox = BBox(); |
|
484 bbox.L = A->P; |
|
485 if (B->P.x < bbox.L.x) bbox.L.x = B->P.x; |
|
486 if (C->P.x < bbox.L.x) bbox.L.x = C->P.x; |
|
487 if (B->P.y < bbox.L.y) bbox.L.y = B->P.y; |
|
488 if (C->P.y < bbox.L.y) bbox.L.y = C->P.y; |
|
489 if (B->P.z < bbox.L.z) bbox.L.z = B->P.z; |
|
490 if (C->P.z < bbox.L.z) bbox.L.z = C->P.z; |
|
491 bbox.H = A->P; |
|
492 if (B->P.x > bbox.H.x) bbox.H.x = B->P.x; |
|
493 if (C->P.x > bbox.H.x) bbox.H.x = C->P.x; |
|
494 if (B->P.y > bbox.H.y) bbox.H.y = B->P.y; |
|
495 if (C->P.y > bbox.H.y) bbox.H.y = C->P.y; |
|
496 if (B->P.z > bbox.H.z) bbox.H.z = B->P.z; |
|
497 if (C->P.z > bbox.H.z) bbox.H.z = C->P.z; |
|
498 return bbox; |
|
499 } |
|
500 |
|
501 ostream & Sphere::dump(ostream &st) const |
|
502 { |
|
503 return st << "(sphere," << center << "," << radius << ")"; |
|
504 } |
|
505 |
|
506 ostream & Box::dump(ostream &st) const |
|
507 { |
|
508 return st << "(box," << L << "," << H << ")"; |
|
509 } |
|
510 |
|
511 ostream & Vertex::dump(ostream &st) const |
|
512 { |
|
513 return st << "(v," << P << ")"; |
|
514 } |
|
515 |
|
516 ostream & NormalVertex::dump(ostream &st) const |
|
517 { |
|
518 return st << "(vn," << P << "," << N << ")"; |
|
519 } |
|
520 |
|
521 ostream & Triangle::dump(ostream &st) const |
|
522 { |
|
523 int idxA, idxB, idxC; |
|
524 if (!vertex_index.get(A, idxA)) |
|
525 st << *A << "," << endl; |
|
526 if (!vertex_index.get(B, idxB)) |
|
527 st << *B << "," << endl; |
|
528 if (!vertex_index.get(C, idxC)) |
|
529 st << *C << "," << endl; |
|
530 return st << "(t," << idxA << "," << idxB << "," << idxC << ")"; |
|
531 } |