using aabb planes as mirror planes
[laserbrain_demo] / src / geom.cc
1 #include <stdlib.h>
2 #include <float.h>
3 #include <algorithm>
4 #include "geom.h"
5 #include "app.h"
6
7 #define SPHERE(ptr)     ((Sphere*)ptr)
8 #define AABOX(ptr)      ((AABox*)ptr)
9 #define BOX(ptr)        ((Box*)ptr)
10 #define PLANE(ptr)      ((Plane*)ptr)
11
12 GeomObject::GeomObject()
13 {
14         type = GOBJ_UNKNOWN;
15 }
16
17 GeomObject::~GeomObject()
18 {
19 }
20
21 bool GeomObject::valid() const
22 {
23         return true;
24 }
25
26 void GeomObject::invalidate()
27 {
28 }
29
30 float GeomObject::distance(const Vec3 &v) const
31 {
32         return sqrt(distance_sq(v));
33 }
34
35 float GeomObject::signed_distance(const Vec3 &v) const
36 {
37         return sqrt(signed_distance_sq(v));
38 }
39
40 Sphere::Sphere()
41 {
42         type = GOBJ_SPHERE;
43         radius = 1.0;
44 }
45
46 Sphere::Sphere(const Vec3 &cent, float radius)
47         : center(cent)
48 {
49         type = GOBJ_SPHERE;
50         this->radius = radius;
51 }
52
53 bool Sphere::valid() const
54 {
55         return radius >= 0.0f;
56 }
57
58 void Sphere::invalidate()
59 {
60         center = Vec3(0, 0, 0);
61         radius = -1;
62 }
63
64 bool Sphere::intersect(const Ray &ray, HitPoint *hit) const
65 {
66         float a = dot(ray.dir, ray.dir);
67         float b = 2.0 * ray.dir.x * (ray.origin.x - center.x) +
68                 2.0 * ray.dir.y * (ray.origin.y - center.y) +
69                 2.0 * ray.dir.z * (ray.origin.z - center.z);
70         float c = dot(ray.origin, ray.origin) + dot(center, center) -
71                 2.0 * dot(ray.origin, center) - radius * radius;
72
73         float discr = b * b - 4.0 * a * c;
74         if(discr < 1e-4) {
75                 return false;
76         }
77
78         float sqrt_discr = sqrt(discr);
79         float t0 = (-b + sqrt_discr) / (2.0 * a);
80         float t1 = (-b - sqrt_discr) / (2.0 * a);
81
82         if(t0 < 1e-4)
83                 t0 = t1;
84         if(t1 < 1e-4)
85                 t1 = t0;
86
87         float t = t0 < t1 ? t0 : t1;
88         if(t < 1e-4) {
89                 return false;
90         }
91
92         // fill the HitPoint structure
93         if(hit) {
94                 hit->obj = this;
95                 hit->dist = t;
96                 hit->pos = ray.origin + ray.dir * t;
97                 hit->normal = (hit->pos - center) / radius;
98         }
99         return true;
100 }
101
102 bool Sphere::contains(const Vec3 &pt) const
103 {
104         return length_sq(pt - center) <= radius * radius;
105 }
106
107 float Sphere::distance_sq(const Vec3 &v) const
108 {
109         return std::max(length_sq(v - center) - radius * radius, 0.0f);
110 }
111
112 float Sphere::signed_distance_sq(const Vec3 &v) const
113 {
114         return length_sq(v - center) - radius * radius;
115 }
116
117 AABox::AABox()
118 {
119         type = GOBJ_AABOX;
120 }
121
122 AABox::AABox(const Vec3 &vmin, const Vec3 &vmax)
123         : min(vmin), max(vmax)
124 {
125         type = GOBJ_AABOX;
126 }
127
128 bool AABox::valid() const
129 {
130         return min.x <= max.x && min.y <= max.y && min.z <= max.z;
131 }
132
133 void AABox::invalidate()
134 {
135         min = Vec3(FLT_MAX, FLT_MAX, FLT_MAX);
136         max = Vec3(-FLT_MAX, -FLT_MAX, -FLT_MAX);
137 }
138
139 Vec3 AABox::get_corner(int idx) const
140 {
141         Vec3 v[] = {min, max};
142         static const int xidx[] = {0, 1, 1, 0, 0, 1, 1, 0};
143         static const int yidx[] = {0, 0, 0, 0, 1, 1, 1, 1};
144         static const int zidx[] = {0, 0, 1, 1, 0, 0, 1, 1};
145         return Vec3(v[xidx[idx]].x, v[yidx[idx]].y, v[zidx[idx]].z);
146 }
147
148 Plane AABox::get_plane(int pidx) const
149 {
150         switch(pidx) {
151         case AABOX_PLANE_PX:
152                 return Plane(Vec3(max.x, 0, 0), Vec3(1, 0, 0));
153         case AABOX_PLANE_NX:
154                 return Plane(Vec3(min.x, 0, 0), Vec3(-1, 0, 0));
155         case AABOX_PLANE_PY:
156                 return Plane(Vec3(0, max.x, 0), Vec3(0, 1, 0));
157         case AABOX_PLANE_NY:
158                 return Plane(Vec3(0, min.x, 0), Vec3(0, -1, 0));
159         case AABOX_PLANE_PZ:
160                 return Plane(Vec3(0, 0, max.z), Vec3(0, 0, 1));
161         case AABOX_PLANE_NZ:
162                 return Plane(Vec3(0, 0, min.z), Vec3(0, 0, -1));
163         default:
164                 break;
165         }
166         abort();
167         return Plane();
168 }
169
170 bool AABox::intersect(const Ray &ray, HitPoint *hit) const
171 {
172         Vec3 param[2] = {min, max};
173 #ifndef NDEBUG
174         Vec3 inv_dir;
175         if(fpexcept_enabled) {
176                 inv_dir.x = ray.dir.x == 0.0f ? 1.0f : 1.0f / ray.dir.x;
177                 inv_dir.y = ray.dir.y == 0.0f ? 1.0f : 1.0f / ray.dir.y;
178                 inv_dir.z = ray.dir.z == 0.0f ? 1.0f : 1.0f / ray.dir.z;
179         } else {
180                 inv_dir = Vec3(1.0 / ray.dir.x, 1.0 / ray.dir.y, 1.0 / ray.dir.z);
181         }
182 #else
183         Vec3 inv_dir(1.0 / ray.dir.x, 1.0 / ray.dir.y, 1.0 / ray.dir.z);
184 #endif
185         int sign[3] = {inv_dir.x < 0, inv_dir.y < 0, inv_dir.z < 0};
186
187         float tmin = (param[sign[0]].x - ray.origin.x) * inv_dir.x;
188         float tmax = (param[1 - sign[0]].x - ray.origin.x) * inv_dir.x;
189         float tymin = (param[sign[1]].y - ray.origin.y) * inv_dir.y;
190         float tymax = (param[1 - sign[1]].y - ray.origin.y) * inv_dir.y;
191
192         if(tmin > tymax || tymin > tmax) {
193                 return false;
194         }
195         if(tymin > tmin) {
196                 tmin = tymin;
197         }
198         if(tymax < tmax) {
199                 tmax = tymax;
200         }
201
202         float tzmin = (param[sign[2]].z - ray.origin.z) * inv_dir.z;
203         float tzmax = (param[1 - sign[2]].z - ray.origin.z) * inv_dir.z;
204
205         if(tmin > tzmax || tzmin > tmax) {
206                 return false;
207         }
208         if(tzmin > tmin) {
209                 tmin = tzmin;
210         }
211         if(tzmax < tmax) {
212                 tmax = tzmax;
213         }
214
215         float t = tmin < 1e-4 ? tmax : tmin;
216         if(t >= 1e-4) {
217
218                 if(hit) {
219                         hit->obj = this;
220                         hit->dist = t;
221                         hit->pos = ray.origin + ray.dir * t;
222
223                         float min_dist = FLT_MAX;
224                         Vec3 offs = min + (max - min) / 2.0;
225                         Vec3 local_hit = hit->pos - offs;
226
227                         static const Vec3 axis[] = {
228                                 Vec3(1, 0, 0), Vec3(0, 1, 0), Vec3(0, 0, 1)
229                         };
230                         //int tcidx[][2] = {{2, 1}, {0, 2}, {0, 1}};
231
232                         for(int i=0; i<3; i++) {
233                                 float dist = fabs((max[i] - offs[i]) - fabs(local_hit[i]));
234                                 if(dist < min_dist) {
235                                         min_dist = dist;
236                                         hit->normal = axis[i] * (local_hit[i] < 0.0 ? 1.0 : -1.0);
237                                         //hit->texcoord = Vec2(hit->pos[tcidx[i][0]], hit->pos[tcidx[i][1]]);
238                                 }
239                         }
240                 }
241                 return true;
242         }
243         return false;
244 }
245
246 bool AABox::contains(const Vec3 &v) const
247 {
248         return v.x >= min.x && v.y >= min.y && v.z >= min.z &&
249                 v.x <= max.x && v.y <= max.y && v.z <= max.z;
250 }
251
252 #define SQ(x)   ((x) * (x))
253 float AABox::distance_sq(const Vec3 &v) const
254 {
255         float dsq = 0.0f;
256
257         for(int i=0; i<3; i++) {
258                 if(v[i] < min[i]) dsq += SQ(min[i] - v[i]);
259                 if(v[i] > max[i]) dsq += SQ(v[i] - max[i]);
260         }
261         return dsq;
262 }
263
264 float AABox::signed_distance_sq(const Vec3 &v) const
265 {
266         if(!contains(v)) {
267                 return distance_sq(v);
268         }
269
270         float dsq = 0.0f;
271         for(int i=0; i<3; i++) {
272                 float dmin = v[i] - min[i];
273                 float dmax = max[i] - v[i];
274                 dsq -= dmin < dmax ? SQ(dmin) : SQ(dmax);
275         }
276         return dsq;
277 }
278
279 Box::Box()
280 {
281         type = GOBJ_BOX;
282 }
283
284 Box::Box(const AABox &aabox, const Mat4 &xform)
285         : xform(xform)
286 {
287         type = GOBJ_BOX;
288         min = aabox.min;
289         max = aabox.max;
290 }
291
292 void Box::invalidate()
293 {
294         AABox::invalidate();
295         xform = Mat4::identity;
296 }
297
298 Box::Box(const Vec3 &min, const Vec3 &max)
299         : AABox(min, max)
300 {
301         type = GOBJ_BOX;
302 }
303
304 Box::Box(const Vec3 &min, const Vec3 &max, const Mat4 &xform)
305         : AABox(min, max), xform(xform)
306 {
307         type = GOBJ_BOX;
308 }
309
310 // XXX all this shit is completely untested
311 Box::Box(const Vec3 &pos, const Vec3 &vi, const Vec3 &vj, const Vec3 &vk)
312 {
313         type = GOBJ_BOX;
314         float ilen = length(vi);
315         float jlen = length(vj);
316         float klen = length(vk);
317
318         min = Vec3(-ilen, -jlen, -klen);
319         max = Vec3(ilen, jlen, klen);
320
321         float si = ilen == 0.0 ? 1.0 : 1.0 / ilen;
322         float sj = jlen == 0.0 ? 1.0 : 1.0 / jlen;
323         float sk = klen == 0.0 ? 1.0 : 1.0 / klen;
324
325         xform = Mat4(vi * si, vj * sj, vk * sk);
326         xform.translate(pos);
327 }
328
329 Box::Box(const Vec3 *varr, int vcount)
330 {
331         type = GOBJ_BOX;
332         calc_bounding_aabox(this, varr, vcount);
333 }
334
335 Vec3 Box::get_corner(int idx) const
336 {
337         return xform * AABox::get_corner(idx);
338 }
339
340 bool Box::intersect(const Ray &ray, HitPoint *hit) const
341 {
342         Mat4 inv_xform = inverse(xform);
343         Mat4 dir_inv_xform = inv_xform.upper3x3();
344         Mat4 dir_xform = transpose(dir_inv_xform);
345         Ray local_ray = Ray(inv_xform * ray.origin, dir_inv_xform * ray.dir);
346
347         bool res = AABox::intersect(local_ray, hit);
348         if(!res || !hit) return res;
349
350         hit->pos = xform * hit->pos;
351         hit->normal = dir_xform * hit->normal;
352         hit->local_ray = local_ray;
353         hit->ray = ray;
354         return true;
355 }
356
357 bool Box::contains(const Vec3 &pt) const
358 {
359         // XXX is it faster to extract 6 planes and do dot products? sounds marginal
360         return AABox::contains(inverse(xform) * pt);
361 }
362
363 float Box::distance_sq(const Vec3 &v) const
364 {
365         return 0.0f;    // TODO
366 }
367
368 float Box::signed_distance_sq(const Vec3 &v) const
369 {
370         return 0.0f;    // TODO
371 }
372
373 Plane::Plane()
374         : normal(0.0, 1.0, 0.0)
375 {
376         type = GOBJ_PLANE;
377 }
378
379 Plane::Plane(const Vec3 &p, const Vec3 &norm)
380         : pt(p)
381 {
382         type = GOBJ_PLANE;
383         normal = normalize(norm);
384 }
385
386 Plane::Plane(const Vec3 &p1, const Vec3 &p2, const Vec3 &p3)
387         : pt(p1)
388 {
389         type = GOBJ_PLANE;
390         normal = normalize(cross(p2 - p1, p3 - p1));
391 }
392
393 Plane::Plane(const Vec3 &normal, float dist)
394 {
395         type = GOBJ_PLANE;
396         this->normal = normalize(normal);
397         pt = this->normal * dist;
398 }
399
400 bool Plane::intersect(const Ray &ray, HitPoint *hit) const
401 {
402         float ndotdir = dot(normal, ray.dir);
403         if(fabs(ndotdir) < 1e-4) {
404                 return false;
405         }
406
407         if(hit) {
408                 Vec3 ptdir = pt - ray.origin;
409                 float t = dot(normal, ptdir) / ndotdir;
410
411                 hit->dist = t;
412                 hit->pos = ray.origin + ray.dir * t;
413                 hit->normal = normal;
414                 hit->obj = this;
415         }
416         return true;
417 }
418
419 bool Plane::contains(const Vec3 &v) const
420 {
421         return dot(v, normal) <= 0.0;
422 }
423
424 float Plane::distance(const Vec3 &v) const
425 {
426         return std::max(dot(v - pt, normal), 0.0f);
427 }
428
429 float Plane::signed_distance(const Vec3 &v) const
430 {
431         return dot(v - pt, normal);
432 }
433
434 float Plane::distance_sq(const Vec3 &v) const
435 {
436         float d = distance(v);
437         return d * d;
438 }
439
440 float Plane::signed_distance_sq(const Vec3 &v) const
441 {
442         float d = distance(v);
443         return dot(v, normal) >= 0.0f ? d * d : -(d * d);
444 }
445
446
447 Disc::Disc()
448 {
449         type = GOBJ_DISC;
450         radius = 1.0;
451 }
452
453 Disc::Disc(const Vec3 &pt, const Vec3 &normal, float rad)
454         : Plane(pt, normal)
455 {
456         type = GOBJ_DISC;
457         radius = rad;
458 }
459
460 Disc::Disc(const Vec3 &normal, float dist, float rad)
461         : Plane(normal, dist)
462 {
463         type = GOBJ_DISC;
464         radius = rad;
465 }
466
467 bool Disc::valid() const
468 {
469         return radius >= 0.0f;
470 }
471
472 void Disc::invalidate()
473 {
474         radius = -1;
475 }
476
477 bool Disc::intersect(const Ray &ray, HitPoint *hit) const
478 {
479         HitPoint phit;
480         if(Plane::intersect(ray, &phit)) {
481                 if(length_sq(phit.pos - pt) <= radius * radius) {
482                         *hit = phit;
483                         return true;
484                 }
485         }
486         return false;
487 }
488
489 bool Disc::contains(const Vec3 &pt) const
490 {
491         Vec3 pj = proj_point_plane(pt, *this);
492         return length_sq(pj - this->pt) <= radius * radius;
493 }
494
495 float Disc::distance_sq(const Vec3 &v) const
496 {
497         return 0.0;     // TODO
498 }
499
500 float Disc::signed_distance_sq(const Vec3 &v) const
501 {
502         return 0.0;     // TODO
503 }
504
505
506 Vec3 proj_point_plane(const Vec3 &pt, const Plane &plane)
507 {
508         float dist = plane.signed_distance(pt);
509         return pt - plane.normal * dist;
510 }
511
512 // ---- bounding sphere calculations ----
513
514 bool calc_bounding_sphere(Sphere *sph, const GeomObject *obj)
515 {
516         if(!obj->valid()) {
517                 sph->invalidate();
518                 return true;
519         }
520
521         switch(obj->type) {
522         case GOBJ_SPHERE:
523                 *sph = *(Sphere*)obj;
524                 break;
525
526         case GOBJ_AABOX:
527                 sph->center = (AABOX(obj)->min + AABOX(obj)->max) * 0.5;
528                 sph->radius = length(AABOX(obj)->max - AABOX(obj)->min) * 0.5;
529                 break;
530
531         case GOBJ_BOX:
532                 sph->center = (BOX(obj)->min + BOX(obj)->max) * 0.5 + BOX(obj)->xform.get_translation();
533                 sph->radius = length(BOX(obj)->max - BOX(obj)->min) * 0.5;
534                 break;
535
536         case GOBJ_PLANE:
537         default:
538                 return false;
539         }
540         return true;
541 }
542
543 bool calc_bounding_sphere(Sphere *sph, const GeomObject *a, const GeomObject *b)
544 {
545         Sphere bsa, bsb;
546
547         if(!calc_bounding_sphere(&bsa, a) || !calc_bounding_sphere(&bsb, b)) {
548                 return false;
549         }
550
551         float dist = length(bsa.center - bsb.center);
552         float surf_dist = dist - (bsa.radius + bsb.radius);
553         float d1 = bsa.radius + surf_dist / 2.0;
554         float d2 = bsb.radius + surf_dist / 2.0;
555         float t = d1 / (d1 + d2);
556
557         if(t < 0.0) t = 0.0;
558         if(t > 1.0) t = 1.0;
559
560         sph->center = bsa.center * t + bsb.center * (1.0 - t);
561         sph->radius = std::max(dist * t + bsb.radius, dist * (1.0f - t) + bsa.radius);
562         return true;
563 }
564
565 bool calc_bounding_sphere(Sphere *sph, const GeomObject **objv, int num)
566 {
567         if(num <= 0) return false;
568
569         if(!calc_bounding_sphere(sph, objv[0])) {
570                 return false;
571         }
572
573         for(int i=1; i<num; i++) {
574                 if(!calc_bounding_sphere(sph, sph, objv[i])) {
575                         return false;
576                 }
577         }
578         return true;
579 }
580
581 bool calc_bounding_sphere(Sphere *sph, const Vec3 *v, int num, const Mat4 &xform)
582 {
583         if(num <= 0) return false;
584
585         sph->center = Vec3(0.0, 0.0, 0.0);
586         for(int i=0; i<num; i++) {
587                 sph->center += xform * v[i];
588         }
589         sph->center /= (float)num;
590
591         float rad_sq = 0.0f;
592         for(int i=0; i<num; i++) {
593                 Vec3 dir = xform * v[i] - sph->center;
594                 rad_sq = std::max(rad_sq, dot(dir, dir));
595         }
596         sph->radius = sqrt(rad_sq);
597         return true;
598 }
599
600 bool calc_bounding_aabox(AABox *box, const GeomObject *obj)
601 {
602         if(!obj->valid()) {
603                 box->invalidate();
604                 return true;
605         }
606
607         switch(obj->type) {
608         case GOBJ_AABOX:
609                 *box = *(AABox*)obj;
610                 break;
611
612         case GOBJ_BOX:
613                 {
614                         Vec3 v[8];
615                         for(int i=0; i<8; i++) {
616                                 v[i] = BOX(obj)->get_corner(i);
617                         }
618                         calc_bounding_aabox(box, v, 8);
619                 }
620                 break;
621
622         case GOBJ_SPHERE:
623                 {
624                         float r = SPHERE(obj)->radius;
625                         box->min = SPHERE(obj)->center - Vec3(r, r, r);
626                         box->max = SPHERE(obj)->center + Vec3(r, r, r);
627                 }
628                 break;
629
630         case GOBJ_PLANE:
631         default:
632                 return false;
633         }
634         return true;
635 }
636
637 bool calc_bounding_aabox(AABox *box, const GeomObject *a, const GeomObject *b)
638 {
639         AABox bba, bbb;
640
641         if(!calc_bounding_aabox(&bba, a) || !calc_bounding_aabox(&bbb, b)) {
642                 return false;
643         }
644
645         for(int i=0; i<3; i++) {
646                 box->min[i] = std::min(bba.min[i], bbb.min[i]);
647                 box->max[i] = std::max(bba.max[i], bbb.max[i]);
648         }
649         return true;
650 }
651
652 bool calc_bounding_aabox(AABox *box, const GeomObject **objv, int num)
653 {
654         if(num <= 0) return false;
655
656         if(!calc_bounding_aabox(box, objv[0])) {
657                 return false;
658         }
659
660         for(int i=1; i<num; i++) {
661                 if(!calc_bounding_aabox(box, box, objv[i])) {
662                         return false;
663                 }
664         }
665         return true;
666 }
667
668 bool calc_bounding_aabox(AABox *box, const Vec3 *v, int num, const Mat4 &xform)
669 {
670         if(num <= 0) return false;
671
672         box->min = box->max = xform * v[0];
673         for(int i=1; i<num; i++) {
674                 Vec3 p = xform * v[i];
675
676                 for(int j=0; j<3; j++) {
677                         box->min[j] = std::min(box->min[j], p[j]);
678                         box->max[j] = std::max(box->max[j], p[j]);
679                 }
680         }
681         return true;
682 }
683
684 bool calc_bounding_box(Box *box, const GeomObject *obj)
685 {
686         if(!obj->valid()) {
687                 box->invalidate();
688                 return true;
689         }
690
691         switch(obj->type) {
692         case GOBJ_BOX:
693                 *box = *(Box*)obj;
694                 break;
695
696         case GOBJ_AABOX:
697                 box->min = BOX(obj)->min;
698                 box->max = BOX(obj)->max;
699                 box->xform = Mat4::identity;
700                 break;
701
702         case GOBJ_SPHERE:
703                 {
704                         float r = SPHERE(obj)->radius;
705                         box->min = SPHERE(obj)->center - Vec3(r, r, r);
706                         box->max = SPHERE(obj)->center + Vec3(r, r, r);
707                         box->xform = Mat4::identity;
708                 }
709                 break;
710
711         case GOBJ_PLANE:
712         default:
713                 return false;
714         }
715         return true;
716 }
717
718 bool intersect_sphere_sphere(Disc *result, const Sphere &a, const Sphere &b)
719 {
720         Vec3 dir = b.center - a.center;
721
722         float dist_sq = length_sq(dir);
723         if(dist_sq <= 1e-8) return false;
724
725         float rsum = a.radius + b.radius;
726         float rdif = fabs(a.radius - b.radius);
727         if(dist_sq > rsum * rsum || dist_sq < rdif * rdif) {
728                 return false;
729         }
730
731         float dist = sqrt(dist_sq);
732         float t = (dist_sq + a.radius * a.radius - b.radius * b.radius) / (2.0 * sqrt(dist_sq));
733
734         result->pt = a.center + dir * t;
735         result->normal = dir / dist;
736         result->radius = sin(acos(t)) * a.radius;
737         return true;
738 }
739
740 bool intersect_plane_plane(Ray *result, const Plane &a, const Plane &b)
741 {
742         return false;   // TODO
743 }
744
745 bool intersect_sphere_plane(Sphere *result, const Sphere &s, const Plane &p)
746 {
747         return false;   // TODO
748 }
749
750 bool intersect_plane_sphere(Sphere *result, const Plane &p, const Sphere &s)
751 {
752         return false;   // TODO
753 }
754
755 bool intersect_aabox_aabox(AABox *res, const AABox &a, const AABox &b)
756 {
757         for(int i=0; i<3; i++) {
758                 res->min[i] = std::max(a.min[i], b.min[i]);
759                 res->max[i] = std::min(a.max[i], b.max[i]);
760
761                 if(res->max[i] < res->min[i]) {
762                         res->max[i] = res->min[i];
763                 }
764         }
765         return res->min.x != res->max.x && res->min.y != res->max.y && res->min.z != res->max.z;
766 }
767
768 bool collision_sphere_aabox(const Sphere &s, const AABox &b)
769 {
770         return b.distance_sq(s.center) <= s.radius * s.radius;
771 }