source: GTP/trunk/Lib/Vis/Preprocessing/src/AxisAlignedBox3.cpp @ 704

Revision 704, 57.4 KB checked in by mattausch, 18 years ago (diff)

implemented first version of the visibility filter
Line 
1
2// GOLEM library
3#include <assert.h>
4#include <iostream>
5using namespace std;
6#include "AxisAlignedBox3.h"
7#include "Ray.h"
8#include "Polygon3.h"
9#include "Mesh.h"
10
11#define FATAL Debug
12#define FATAL_ABORT exit(1)
13
14
15// AxisAlignedBox3 implementations
16
17// Overload << operator for C++-style output
18ostream&
19operator<< (ostream &s, const AxisAlignedBox3 &A)
20{
21  return s << '[' << A.mMin.x << ", " << A.mMin.y << ", " << A.mMin.z << "]["
22           << A.mMax.x << ", " << A.mMax.y << ", " << A.mMax.z << ']';
23}
24
25// Overload >> operator for C++-style input
26istream&
27operator>> (istream &s, AxisAlignedBox3 &A)
28{
29  char a;
30  // read "[min.x, min.y, min.z][mMax.x, mMax.y, mMax.z]"
31  return s >> a >> A.mMin.x >> a >> A.mMin.y >> a >> A.mMin.z >> a >> a
32    >> A.mMax.x >> a >> A.mMax.y >> a >> A.mMax.z >> a;
33}
34
35bool
36AxisAlignedBox3::Unbounded() const
37{
38  return (mMin == Vector3(-MAXFLOAT)) ||
39         (mMax == Vector3(-MAXFLOAT));
40}
41
42void
43AxisAlignedBox3::Include(const Vector3 &newpt)
44{
45  Minimize(mMin, newpt);
46  Maximize(mMax, newpt);
47}
48
49void
50AxisAlignedBox3::Include(const Polygon3 &newpoly)
51{
52        VertexContainer::const_iterator it, it_end = newpoly.mVertices.end();
53
54        for (it = newpoly.mVertices.begin(); it != it_end; ++ it)
55                Include(*it);
56}
57
58void
59AxisAlignedBox3::Include(Mesh *mesh)
60{
61        VertexContainer::const_iterator it, it_end = mesh->mVertices.end();
62
63        for (it = mesh->mVertices.begin(); it != it_end; ++ it)
64                Include(*it);
65}
66
67
68void
69AxisAlignedBox3::Include(const AxisAlignedBox3 &bbox)
70{
71  Minimize(mMin, bbox.mMin);
72  Maximize(mMax, bbox.mMax);
73}
74
75bool
76AxisAlignedBox3::IsCorrect()
77{
78  if ( (mMin.x > mMax.x) ||
79       (mMin.y > mMax.y) ||
80       (mMin.z > mMax.z) )
81    return false; // box is not formed
82  return true;
83}
84
85void
86AxisAlignedBox3::GetEdge(const int edge, Vector3 *a, Vector3 *b) const
87{
88  switch(edge) {
89  case 0:
90    a->SetValue(mMin.x, mMin.y, mMin.z);
91    b->SetValue(mMin.x, mMin.y, mMax.z);
92    break;
93  case 1:
94    a->SetValue(mMin.x, mMin.y, mMin.z);
95    b->SetValue(mMin.x, mMax.y, mMin.z);
96    break;
97  case 2:
98    a->SetValue(mMin.x, mMin.y, mMin.z);
99    b->SetValue(mMax.x, mMin.y, mMin.z);
100    break;
101  case 3:
102    a->SetValue(mMax.x, mMax.y, mMax.z);
103    b->SetValue(mMax.x, mMax.y, mMin.z);
104    break;
105  case 4:
106    a->SetValue(mMax.x, mMax.y, mMax.z);
107    b->SetValue(mMax.x, mMin.y, mMax.z);
108    break;
109  case 5:
110    a->SetValue(mMax.x, mMax.y, mMax.z);
111    b->SetValue(mMin.x, mMax.y, mMax.z);
112    break;
113   
114  case 6:
115    a->SetValue(mMin.x, mMin.y, mMax.z);
116    b->SetValue(mMin.x, mMax.y, mMax.z);
117    break;
118  case 7:
119    a->SetValue(mMin.x, mMin.y, mMax.z);
120    b->SetValue(mMax.x, mMin.y, mMax.z);
121    break;
122  case 8:
123    a->SetValue(mMin.x, mMax.y, mMin.z);
124    b->SetValue(mMin.x, mMax.y, mMax.z);
125    break;
126  case 9:
127    a->SetValue(mMin.x, mMax.y, mMin.z);
128    b->SetValue(mMax.x, mMax.y, mMin.z);
129    break;
130  case 10:
131    a->SetValue(mMax.x, mMin.y, mMin.z);
132    b->SetValue(mMax.x, mMax.y, mMin.z);
133    break;
134  case 11:
135    a->SetValue(mMax.x, mMin.y, mMin.z);
136    b->SetValue(mMax.x, mMin.y, mMax.z);
137    break;
138  }
139}
140
141// returns the vertex indices in the range <0..7>, v = 4.x + 2.y + z, where
142// x,y,z are either 0 or 1; (0 .. min coordinate, 1 .. max coordinate)
143void
144AxisAlignedBox3::GetEdge(const int edge, int  &aIdx, int &bIdx) const
145{
146  switch(edge) {
147  case 0: aIdx = 0; bIdx = 1; break;
148  case 1: aIdx = 0; bIdx = 2; break;
149  case 2: aIdx = 0; bIdx = 4; break;
150  case 3: aIdx = 7; bIdx = 6; break;
151  case 4: aIdx = 7; bIdx = 5; break;
152  case 5: aIdx = 7; bIdx = 3; break;
153  case 6: aIdx = 1; bIdx = 3; break;
154  case 7: aIdx = 1; bIdx = 5; break;
155  case 8: aIdx = 2; bIdx = 3; break;
156  case 9: aIdx = 2; bIdx = 6; break;
157  case 10: aIdx = 4; bIdx = 6; break;
158  case 11: aIdx = 4; bIdx = 5; break;
159  }
160}
161
162void
163AxisAlignedBox3::Include(const int &axis, const float &newBound)
164{
165  switch (axis) {
166    case 0: { // x-axis
167      if (mMin.x > newBound)
168        mMin.x = newBound;
169      if (mMax.x < newBound)
170        mMax.x = newBound;
171      break;
172    }
173    case 1: { // y-axis
174      if (mMin.y > newBound)
175        mMin.y = newBound;
176      if (mMax.y < newBound)
177        mMax.y = newBound;
178      break;
179    }
180    case 2: { // z-axis
181      if (mMin.z > newBound)
182        mMin.z = newBound;
183      if (mMax.z < newBound)
184        mMax.z = newBound;
185      break;
186    }
187  }
188}
189
190#if 0
191// ComputeMinMaxT computes the minimum and maximum signed distances
192// of intersection with the ray; it returns 1 if the ray hits
193// the bounding box and 0 if it does not.
194int
195AxisAlignedBox3::ComputeMinMaxT(const Ray &ray,
196                                float *tmin,
197                                float *tmax) const
198{
199  float minx, maxx, miny, maxy, minz, maxz;
200
201  if (fabs(ray.dir.x) < 0.001) {
202    if (mMin.x < ray.loc.x && mMax.x > ray.loc.x) {
203      minx = -MAXFLOAT;
204      maxx = MAXFLOAT;
205    }
206    else
207      return 0;
208  }
209  else {
210    float t1 = (mMin.x - ray.loc.x) / ray.dir.x;
211    float t2 = (mMax.x - ray.loc.x) / ray.dir.x;
212    if (t1 < t2) {
213      minx = t1;
214      maxx = t2;
215    }
216    else {
217      minx = t2;
218      maxx = t1;
219    }
220    if (maxx < 0)
221      return 0;
222  }
223
224  if (fabs(ray.dir.y) < 0.001) {
225    if (mMin.y < ray.loc.y && mMax.y > ray.loc.y) {
226      miny = -MAXFLOAT;
227      maxy = MAXFLOAT;
228    }
229    else
230      return 0;
231  }
232  else {
233    float t1 = (mMin.y - ray.loc.y) / ray.dir.y;
234    float t2 = (mMax.y - ray.loc.y) / ray.dir.y;
235    if (t1 < t2) {
236      miny = t1;
237      maxy = t2;
238    }
239    else {
240      miny = t2;
241      maxy = t1;
242    }
243    if (maxy < 0.0)
244      return 0;
245  }
246
247  if (fabs(ray.dir.z) < 0.001) {
248    if (mMin.z < ray.loc.z && mMax.z > ray.loc.z) {
249      minz = -MAXFLOAT;
250      maxz = MAXFLOAT;
251    }
252    else
253      return 0;
254  }
255  else {
256    float t1 = (mMin.z - ray.loc.z) / ray.dir.z;
257    float t2 = (mMax.z - ray.loc.z) / ray.dir.z;
258    if (t1 < t2) {
259      minz = t1;
260      maxz = t2;
261    }
262    else {
263      minz = t2;
264      maxz = t1;
265    }
266    if (maxz < 0.0)
267      return 0;
268  }
269
270  *tmin = minx;
271  if (miny > *tmin)
272    *tmin = miny;
273  if (minz > *tmin)
274    *tmin = minz;
275
276  *tmax = maxx;
277  if (maxy < *tmax)
278    *tmax = maxy;
279  if (maxz < *tmax)
280    *tmax = maxz;
281
282  return 1; // yes, intersection was found
283}
284#else
285// another variant of the same, with less variables
286int
287AxisAlignedBox3::ComputeMinMaxT(const Ray &ray,
288                                float *tmin,
289                                float *tmax) const
290{
291  const float dirEps = 1e-8f;
292  register float minx, maxx;
293  ray.ComputeInvertedDir();
294 
295  if (fabs(ray.dir.x) < dirEps) {
296    if (mMin.x < ray.loc.x && mMax.x > ray.loc.x) {
297      minx = -MAXFLOAT;
298      maxx = MAXFLOAT;
299    }
300    else
301      return 0;
302  }
303  else {
304    float t1 = (mMin.x - ray.loc.x) * ray.invDir.x;
305    float t2 = (mMax.x - ray.loc.x) * ray.invDir.x;
306    if (t1 < t2) {
307      minx = t1;
308      maxx = t2;
309    }
310    else {
311      minx = t2;
312      maxx = t1;
313    }
314        //    if (maxx < 0.0)
315        //      return 0;
316  }
317
318  *tmin = minx;
319  *tmax = maxx;
320 
321  if (fabs(ray.dir.y) < dirEps) {
322    if (mMin.y < ray.loc.y && mMax.y > ray.loc.y) {
323      minx = -MAXFLOAT;
324      maxx = MAXFLOAT;
325    }
326    else
327      return 0;
328  }
329  else {
330    float t1 = (mMin.y - ray.loc.y) * ray.invDir.y;
331    float t2 = (mMax.y - ray.loc.y) * ray.invDir.y;
332    if (t1 < t2) {
333      minx = t1;
334      maxx = t2;
335    }
336    else {
337      minx = t2;
338      maxx = t1;
339    }
340        //    if (maxx < 0.0)
341        //      return 0;
342  }
343
344  if (minx > *tmin)
345    *tmin = minx;
346  if (maxx < *tmax)
347    *tmax = maxx;
348 
349  if (fabs(ray.dir.z) < dirEps) {
350    if (mMin.z < ray.loc.z && mMax.z > ray.loc.z) {
351      minx = -MAXFLOAT;
352      maxx = MAXFLOAT;
353    }
354    else
355      return 0;
356  }
357  else {
358    float t1 = (mMin.z - ray.loc.z) * ray.invDir.z;
359    float t2 = (mMax.z - ray.loc.z) * ray.invDir.z;
360    if (t1 < t2) {
361      minx = t1;
362      maxx = t2;
363    }
364    else {
365      minx = t2;
366      maxx = t1;
367    }
368        //    if (maxx < 0.0)
369        //      return 0;
370  }
371
372  if (minx > *tmin)
373    *tmin = minx;
374  if (maxx < *tmax)
375    *tmax = maxx;
376
377  return 1; // yes, intersection was found
378}
379#endif
380
381// ComputeMinMaxT computes the minimum and maximum parameters
382// of intersection with the ray; it returns 1 if the ray hits
383// the bounding box and 0 if it does not.
384int
385AxisAlignedBox3::ComputeMinMaxT(const Ray &ray, float *tmin, float *tmax,
386                                EFaces &entryFace, EFaces &exitFace) const
387{
388  float minx, maxx, miny, maxy, minz, maxz;
389  int swapped[3];
390 
391  if (fabs(ray.dir.x) < 0.001) {
392    if (mMin.x < ray.loc.x && mMax.x > ray.loc.x) {
393      minx = -MAXFLOAT;
394      maxx = MAXFLOAT;
395    }
396    else
397      return 0;
398  }
399  else {
400    float t1 = (mMin.x - ray.loc.x) / ray.dir.x;
401    float t2 = (mMax.x - ray.loc.x) / ray.dir.x;
402    if (t1 < t2) {
403      minx = t1;
404      maxx = t2;
405      swapped[0] = 0;
406    }
407    else {
408      minx = t2;
409      maxx = t1;
410      swapped[0] = 1;
411    }
412    if (maxx < 0)
413      return 0;
414  }
415
416  if (fabs(ray.dir.y) < 0.001) {
417    if (mMin.y < ray.loc.y && mMax.y > ray.loc.y) {
418      miny = -MAXFLOAT;
419      maxy = MAXFLOAT;
420    }
421    else
422      return 0;
423  }
424  else {
425    float t1 = (mMin.y - ray.loc.y) / ray.dir.y;
426    float t2 = (mMax.y - ray.loc.y) / ray.dir.y;
427    if (t1 < t2) {
428      miny = t1;
429      maxy = t2;
430      swapped[1] = 0;
431    }
432    else {
433      miny = t2;
434      maxy = t1;
435      swapped[1] = 1;
436    }
437    if (maxy < 0.0)
438      return 0;
439  }
440
441  if (fabs(ray.dir.z) < 0.001) {
442    if (mMin.z < ray.loc.z && mMax.z > ray.loc.z) {
443      minz = -MAXFLOAT;
444      maxz = MAXFLOAT;
445    }
446    else
447      return 0;
448  }
449  else {
450    float t1 = (mMin.z - ray.loc.z) / ray.dir.z;
451    float t2 = (mMax.z - ray.loc.z) / ray.dir.z;
452    if (t1 < t2) {
453      minz = t1;
454      maxz = t2;
455      swapped[2] = 0;
456    }
457    else {
458      minz = t2;
459      maxz = t1;
460      swapped[2] = 1;
461    }
462    if (maxz < 0.0)
463      return 0;
464  }
465
466  *tmin = minx;
467  entryFace = ID_Back;
468  if (miny > *tmin) {
469    *tmin = miny;
470    entryFace = ID_Left;
471  }
472  if (minz > *tmin) {
473    *tmin = minz;
474    entryFace = ID_Bottom;
475  }
476 
477  *tmax = maxx;
478  exitFace = ID_Back;
479  if (maxy < *tmax) {
480    *tmax = maxy;
481    exitFace = ID_Left;
482  }
483  if (maxz < *tmax) {
484    *tmax = maxz;
485    exitFace = ID_Bottom;
486  }
487
488  if (swapped[entryFace])
489    entryFace = (EFaces)(entryFace + 3);
490 
491  if (!swapped[exitFace])
492    exitFace = (EFaces)(exitFace + 3);
493 
494  return 1; // yes, intersection was found
495}
496
497void
498AxisAlignedBox3::Describe(ostream &app, int ind) const
499{
500  indent(app, ind);
501  app << "AxisAlignedBox3: min at(" << mMin << "), max at(" << mMax << ")\n";
502}
503
504// computes the passing through parameters for case tmin<tmax and tmax>0
505int
506AxisAlignedBox3::GetMinMaxT(const Ray &ray, float *tmin, float *tmax) const
507{
508  if (!ComputeMinMaxT(ray, tmin, tmax))
509    return 0;
510  if ( *tmax < *tmin)
511    return 0; // the ray passes outside the box
512
513  if ( *tmax < 0.0)
514    return 0; // the intersection is not on the positive halfline
515
516  return 1; // ray hits the box .. origin can be outside or inside the box
517}
518
519// computes the signed distances for case tmin<tmax and tmax>0
520int
521AxisAlignedBox3::GetMinMaxT(const Ray &ray, float *tmin, float *tmax,
522                            EFaces &entryFace, EFaces &exitFace) const
523{
524  if (!ComputeMinMaxT(ray, tmin, tmax, entryFace, exitFace))
525    return 0;
526  if ( *tmax < *tmin)
527    return 0; // the ray passes outside the box
528 
529  if ( *tmax < 0.0)
530    return 0; // the intersection is not on the positive halfline
531 
532  return 1; // ray hits the box .. origin can be outside or inside the box
533}
534
535#if 0
536int
537AxisAlignedBox3::IsInside(const Vector3 &point) const
538{
539  return (point.x >= mMin.x) && (point.x <= mMax.x) &&
540         (point.y >= mMin.y) && (point.y <= mMax.y) &&
541         (point.z >= mMin.z) && (point.z <= mMax.z);
542}
543#else
544int
545AxisAlignedBox3::IsInside(const Vector3 &v) const
546{
547  return ! (v.x < mMin.x ||
548            v.x > mMax.x ||
549            v.y < mMin.y ||
550            v.y > mMax.y ||
551            v.z < mMin.z ||
552            v.z > mMax.z);
553}
554#endif
555
556bool
557AxisAlignedBox3::Includes(const AxisAlignedBox3 &b) const
558{
559  return (b.mMin.x >= mMin.x &&
560      b.mMin.y >= mMin.y &&
561      b.mMin.z >= mMin.z &&
562      b.mMax.x <= mMax.x &&
563      b.mMax.y <= mMax.y &&
564      b.mMax.z <= mMax.z);
565
566}
567
568
569// compute the coordinates of one vertex of the box
570Vector3
571AxisAlignedBox3::GetVertex(int xAxis, int yAxis, int zAxis) const
572{
573  Vector3 p;
574  if (xAxis)
575    p.x = mMax.x;
576  else
577    p.x = mMin.x;
578
579  if (yAxis)
580    p.y = mMax.y;
581  else
582    p.y = mMin.y;
583
584  if (zAxis)
585    p.z = mMax.z;
586  else
587    p.z = mMin.z;
588  return p;
589}
590
591// compute the vertex for number N = <0..7>, N = 4.x + 2.y + z, where
592// x,y,z are either 0 or 1; (0 .. min coordinate, 1 .. max coordinate)
593void
594AxisAlignedBox3::GetVertex(const int N, Vector3 &vertex) const
595{
596  switch (N) {
597    case 0: vertex = mMin; break;
598    case 1: vertex.SetValue(mMin.x, mMin.y, mMax.z); break;
599    case 2: vertex.SetValue(mMin.x, mMax.y, mMin.z); break;
600    case 3: vertex.SetValue(mMin.x, mMax.y, mMax.z); break;
601    case 4: vertex.SetValue(mMax.x, mMin.y, mMin.z); break;
602    case 5: vertex.SetValue(mMax.x, mMin.y, mMax.z); break;
603    case 6: vertex.SetValue(mMax.x, mMax.y, mMin.z); break;
604    case 7: vertex = mMax; break;
605    default: {
606      FATAL << "ERROR in AxisAlignedBox3::GetVertex N=" << N <<  "\n";
607      FATAL_ABORT;
608    }
609  }
610}
611
612// Returns region 0 .. 26 ; R = 9*x + 3*y + z ; (x,y,z) \in {0,1,2}
613int
614AxisAlignedBox3::GetRegionID(const Vector3 &point) const
615{
616  int ID = 0;
617
618  if (point.z >= mMin.z) {
619    if (point.z <= mMax.z)
620      ID += 1; // inside the two boundary planes
621    else
622      ID += 2; // outside
623  }
624
625  if (point.y >= mMin.y) {
626    if (point.y <= mMax.y)
627      ID += 3; // inside the two boundary planes
628    else
629      ID += 6; // outside
630  }
631
632  if (point.x >= mMin.x) {
633    if (point.x <= mMax.x)
634      ID += 9; // inside the two boundary planes
635    else
636      ID += 18; // outside
637  }
638  return ID;
639}
640
641
642
643// computes if a given box (smaller) lies at least in one
644// projection whole in the box (larger) = this encompasses given
645// ;-----;
646// | ;-; |
647// | '-' |
648// '-----'
649int
650AxisAlignedBox3::IsPiercedByBox(const AxisAlignedBox3 &box, int &axis) const
651{
652  // test on x-axis
653  if ( (mMax.x < box.mMin.x) ||
654       (mMin.x > box.mMax.x) )
655    return 0; // the boxes do not overlap at all at x-axis
656  if ( (box.mMin.y > mMin.y) &&
657       (box.mMax.y < mMax.y) &&
658       (box.mMin.z > mMin.z) &&
659       (box.mMax.z < mMax.z) ) {
660    axis = 0;
661    return 1; // the boxes overlap in x-axis
662  }
663  // test on y-axis
664  if ( (mMax.y < box.mMin.y) ||
665       (mMin.y > box.mMax.y) )
666    return 0; // the boxes do not overlap at all at y-axis
667  if ( (box.mMin.x > mMin.x) &&
668       (box.mMax.x < mMax.x) &&
669       (box.mMin.z > mMin.z) &&
670       (box.mMax.z < mMax.z) ) {
671    axis = 1;
672    return 1; // the boxes overlap in y-axis
673  }
674  // test on z-axis
675  if ( (mMax.z < box.mMin.z) ||
676       (mMin.z > box.mMax.z) )
677    return 0; // the boxes do not overlap at all at y-axis
678  if ( (box.mMin.x > mMin.x) &&
679       (box.mMax.x < mMax.x) &&
680       (box.mMin.y > mMin.y) &&
681       (box.mMax.y < mMax.y) ) {
682    axis = 2;
683    return 1; // the boxes overlap in z-axis
684  }
685  return 0;
686}
687
688float
689AxisAlignedBox3::SurfaceArea() const
690{
691  Vector3 ext = mMax - mMin;
692 
693  return 2.0f * (ext.x * ext.y +
694                ext.x * ext.z +
695                ext.y * ext.z);
696}
697
698const int AxisAlignedBox3::bvertices[27][9] =
699{                   // region number.. position
700  {5,1,3,2,6,4,-1,-1,-1},      // 0 .. x=0 y=0 z=0
701  {4,5,1,3,2,0,-1,-1,-1},      // 1 .. x=0 y=0 z=1
702  {4,5,7,3,2,0,-1,-1,-1},      // 2 .. x=0 y=0 z=2
703
704  {0,1,3,2,6,4,-1,-1,-1},      // 3 .. x=0 y=1 z=0
705  {0,1,3,2,-1,-1,-1,-1,-1},    // 4 .. x=0 y=1 z=1
706  {1,5,7,3,2,0,-1,-1,-1},      // 5 .. x=0 y=1 z=2
707
708  {0,1,3,7,6,4,-1,-1,-1},      // 6 .. x=0 y=2 z=0
709  {0,1,3,7,6,2,-1,-1,-1},      // 7 .. x=0 y=2 z=1
710  {1,5,7,6,2,0,-1,-1,-1},      // 8 .. x=0 y=2 z=2
711
712  // the regions number <9,17>
713  {5,1,0,2,6,4,-1,-1,-1},      // 9  .. x=1 y=0 z=0
714  {5,1,0,4,-1,-1,-1,-1,-1},    // 10 .. x=1 y=0 z=1
715  {7,3,1,0,4,5,-1,-1,-1},      // 11 .. x=1 y=0 z=2
716
717  {4,0,2,6,-1,-1,-1,-1,-1},    // 12 .. x=1 y=1 z=0
718  {0,2,3,1,5,4,6,7,-1},        // 13 .. x=1 y=1 z=1 .. inside the box
719  {1,5,7,3,-1,-1,-1,-1,-1},    // 14 .. x=1 y=1 z=2
720
721  {4,0,2,3,7,6,-1,-1,-1},      // 15 .. x=1 y=2 z=0
722  {6,2,3,7,-1,-1,-1,-1,-1},    // 16 .. x=1 y=2 z=1
723  {1,5,7,6,2,3,-1,-1,-1},      // 17 .. x=1 y=2 z=2
724
725  // the regions number <18,26>
726  {1,0,2,6,7,5,-1,-1,-1},      // 18 .. x=2 y=0 z=0
727  {1,0,4,6,7,5,-1,-1,-1},      // 19 .. x=2 y=0 z=1
728  {0,4,6,7,3,1,-1,-1,-1},      // 20 .. x=2 y=0 z=2
729
730  {4,0,2,6,7,5,-1,-1,-1},      // 21 .. x=2 y=1 z=0
731  {5,4,6,7,-1,-1,-1,-1,-1},    // 22 .. x=2 y=1 z=1
732  {5,4,6,7,3,1,-1,-1,-1},      // 23 .. x=2 y=1 z=2
733
734  {4,0,2,3,7,5,-1,-1,-1},      // 24 .. x=2 y=2 z=0
735  {5,4,6,2,3,7,-1,-1,-1},      // 25 .. x=2 y=2 z=1
736  {5,4,6,2,3,1,-1,-1,-1},      // 26 .. x=2 y=2 z=2
737};
738
739// the visibility of boundary faces from a given region
740// one to three triples: (axis, min_vertex, max_vertex), axis==-1(terminator)
741const int AxisAlignedBox3::bfaces[27][10] =
742{                              // region number .. position
743  {0,0,3,1,0,5,2,0,6,-1},      // 0 .. x=0 y=0 z=0
744  {0,0,3,1,0,5,-1,-1,-1,-1},   // 1 .. x=0 y=0 z=1
745  {0,0,3,1,0,5,2,1,7,-1},      // 2 .. x=0 y=0 z=2
746
747  {0,0,3,2,0,6,-1,-1,-1,-1},   // 3 .. x=0 y=1 z=0
748  {0,0,3,-1,-1,-1,-1,-1,-1,-1},// 4 .. x=0 y=1 z=1
749  {0,0,3,2,1,7,-1,-1,-1,-1},   // 5 .. x=0 y=1 z=2
750
751  {0,0,3,1,2,7,2,0,6,-1},      // 6 .. x=0 y=2 z=0
752  {0,0,3,1,2,7,-1,-1,-1,-1},   // 7 .. x=0 y=2 z=1
753  {0,0,3,1,2,7,2,1,7,-1},      // 8 .. x=0 y=2 z=2
754
755  // the regions number <9,17>
756  {1,0,5,2,0,6,-1,-1,-1,-1},   // 9  .. x=1 y=0 z=0
757  {1,0,5,-1,-1,-1,-1,-1,-1,-1},// 10 .. x=1 y=0 z=1
758  {1,0,5,2,1,7,-1,-1,-1,-1},   // 11 .. x=1 y=0 z=2
759
760  {2,0,6,-1,-1,-1,-1,-1,-1,-1},// 12 .. x=1 y=1 z=0
761  {-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},// 13 .. x=1 y=1 z=1 .. inside the box
762  {2,1,7,-1,-1,-1,-1,-1,-1,-1},// 14 .. x=1 y=1 z=2
763
764  {1,2,7,2,0,6,-1,-1,-1,-1},   // 15 .. x=1 y=2 z=0
765  {1,2,7,-1,-1,-1,-1,-1,-1,-1},// 16 .. x=1 y=2 z=1
766  {1,2,7,2,1,7,-1,-1,-1,-1},   // 17 .. x=1 y=2 z=2
767
768  // the region number <18,26>
769  {0,4,7,1,0,5,2,0,6,-1},      // 18 .. x=2 y=0 z=0
770  {0,4,7,1,0,5,-1,-1,-1,-1},   // 19 .. x=2 y=0 z=1
771  {0,4,7,1,0,5,2,1,7,-1},      // 20 .. x=2 y=0 z=2
772
773  {0,4,7,2,0,6,-1,-1,-1,-1},   // 21 .. x=2 y=1 z=0
774  {0,4,7,-1,-1,-1,-1,-1,-1,-1},// 22 .. x=2 y=1 z=1
775  {0,4,7,2,1,7,-1,-1,-1,-1},   // 23 .. x=2 y=1 z=2
776
777  {0,4,7,1,2,7,2,0,6,-1},      // 24 .. x=2 y=2 z=0
778  {0,4,7,1,2,7,-1,-1,-1,-1},   // 25 .. x=2 y=2 z=1
779  {0,4,7,1,2,7,2,1,7,-1},      // 26 .. x=2 y=2 z=2
780};
781
782// the correct corners indexing from entry face to exit face
783// first index determines entry face, second index exit face, and
784// the two numbers (indx, inc) determines: ind = the index on the exit
785// face, when starting from the vertex 0 on entry face, 'inc' is
786// the increment when we go on entry face in order 0,1,2,3 to create
787// convex shaft with the rectangle on exit face. That is, inc = -1 or 1.
788const int AxisAlignedBox3::pairFaceRects[6][6][2] = {
789  { // entry face = 0
790    {-1,0}, // exit face 0 .. no meaning
791    {0,-1}, // 1
792    {0,-1}, // 2
793    {0,1}, // 3 .. opposite face
794    {3,1}, // 4
795    {1,1} // 5
796  },
797  { // entry face = 1
798    {0,-1}, // exit face 0
799    {-1,0}, // 1 .. no meaning
800    {0,-1}, // 2
801    {1,1}, // 3
802    {0,1}, // 4 .. opposite face
803    {3,1} // 5
804  },
805  { // entry face = 2
806    {0,-1}, // 0
807    {0,-1}, // 1
808    {-1,0}, // 2 .. no meaning
809    {3,1}, // 3
810    {1,1}, // 4
811    {0,1} // 5 .. opposite face
812  },
813  { // entry face = 3
814    {0,1}, // 0 .. opposite face
815    {3,-1}, // 1
816    {1,1}, // 2
817    {-1,0}, // 3  .. no meaning
818    {0,-1}, // 4
819    {0,-1} // 5
820  },
821  { // entry face = 4
822    {1,1}, // 0
823    {0,1}, // 1 .. opposite face
824    {3,1}, // 2
825    {0,-1}, // 3
826    {-1,0}, // 4 .. no meaning
827    {0,-1} // 5
828  },
829  { // entry face = 5
830    {3,-1}, // 0
831    {1,1}, // 1
832    {0,1}, // 2  .. opposite face
833    {0,-1}, // 3
834    {0,-1}, // 4
835    {-1,0} // 5 .. no meaning
836  }
837};
838
839
840// ------------------------------------------------------------
841// The vertices that form CLOSEST points with respect to the region
842// for all the regions possible, number of regions is 3^3 = 27,
843// since two parallel sides of bbox forms three disjoint spaces.
844// The vertices are given in anti-clockwise order, stopped by value 15,
845// at most 8 points, at least 1 point.
846// The table includes the closest 1/2/4/8 points, followed possibly
847// by the set of coordinates that should be used for testing for
848// the proximity queries. The coordinates to be tested are described by
849// the pair (a,b), when a=0, we want to test min vector of the box,
850//                 when a=1, we want to test max vector of the box
851//                 b=0,1,2 corresponds to the axis (0=x,1=y,2=z)
852// The sequence is ended by 15, number -1 is used as the separator
853// between the vertices and coordinates.
854const int
855AxisAlignedBox3::cvertices[27][9] =
856{                   // region number.. position
857  {0,15,15,15,15,15,15,15,15}, // 0 .. x=0 y=0 z=0 D one vertex
858  {0,1,-1,0,0,0,1,15,15},      // 1 .. x=0 y=0 z=1 D two vertices foll. by 2
859  {1,15,15,15,15,15,15,15,15}, // 2 .. x=0 y=0 z=2 D one vertex
860
861  {0,2,-1,0,0,0,2,15,15},      // 3 .. x=0 y=1 z=0 D
862  {0,1,3,2,-1,0,0,15,15},      // 4 .. x=0 y=1 z=1 D
863  {1,3,-1,0,0,1,2,15,15},      // 5 .. x=0 y=1 z=2 D
864
865  {2,15,15,15,15,15,15,15,15}, // 6 .. x=0 y=2 z=0 D
866  {2,3,-1,0,0,1,1,15,15},      // 7 .. x=0 y=2 z=1 D
867  {3,15,15,15,15,15,15,15,15}, // 8 .. x=0 y=2 z=2 D
868
869  // the regions number <9,17>
870  {0,4,-1,0,1,0,2,15,15},      // 9  .. x=1 y=0 z=0 D
871  {5,1,0,4,-1,0,1,15,15},      // 10 .. x=1 y=0 z=1 D
872  {1,5,-1,0,1,1,2,15,15},      // 11 .. x=1 y=0 z=2 D
873
874  {4,0,2,6,-1,0,2,15,15},      // 12 .. x=1 y=1 z=0 D
875  {0,2,3,1,5,4,6,7,15},        // 13 .. x=1 y=1 z=1 .. inside the box
876  {1,5,7,3,-1,1,2,15,15},      // 14 .. x=1 y=1 z=2 D
877
878  {6,2,-1,0,2,1,1,15,15},      // 15 .. x=1 y=2 z=0 D
879  {6,2,3,7,-1,1,1,15,15},      // 16 .. x=1 y=2 z=1 D
880  {3,7,-1,1,1,1,2,15,15},      // 17 .. x=1 y=2 z=2 D
881
882  // the regions number <18,26>
883  {4,15,15,15,15,15,15,15,15}, // 18 .. x=2 y=0 z=0 D
884  {4,5,-1,0,1,1,0,15,15},      // 19 .. x=2 y=0 z=1 D
885  {5,15,15,15,15,15,15,15,15}, // 20 .. x=2 y=0 z=2 D
886
887  {4,6,-1,0,2,1,0,15,15},      // 21 .. x=2 y=1 z=0 D
888  {5,4,6,7,-1,1,0,15,15},      // 22 .. x=2 y=1 z=1 D
889  {7,5,-1,1,0,1,2,15,15},      // 23 .. x=2 y=1 z=2 D
890
891  {6,15,15,15,15,15,15,15,15}, // 24 .. x=2 y=2 z=0 D
892  {6,7,-1,1,0,1,1,15,15},      // 25 .. x=2 y=2 z=1 D
893  {7,15,15,15,15,15,15,15,15}, // 26 .. x=2 y=2 z=2 D
894};
895
896// Table for Sphere-AABB intersection based on the region knowledge
897// Similar array to previous cvertices, but we omit the surfaces
898// which are not necessary for testing. First are vertices,
899// they are finished with -1. Second, there are indexes in
900// the pair (a,b), when a=0, we want to test min vector of the box,
901//                 when a=1, we want to test max vector of the box
902//                 b=0,1,2 corresponds to the axis (0=x,1=y,2=z)
903//
904// So either we check the vertices or only the distance in specified
905// dimensions. There are at all four possible cases:
906//
907//   1) we check one vertex - then sequence start with non-negative index
908//      and is finished with 15
909//   2) we check two coordinates of min/max vector describe by the pair
910//         (a,b) .. a=min/max(0/1) b=x/y/z (0/1/2), sequence starts with 8
911//      and finishes with 15
912//   3) we check only one coordinate of min/max, as for 2), sequence start
913//      with 9 and ends with 15
914//   4) Position 13 - sphere is inside the box, intersection always exist
915//        the sequence start with 15 .. no further testing is necessary
916//        in this case
917const int
918AxisAlignedBox3::csvertices[27][6] =
919{                   // region number.. position
920  {0,15,15,15,15,15},  // 0 .. x=0 y=0 z=0 D vertex only
921  {8,0,0,0,1,15},      // 1 .. x=0 y=0 z=1 D two coords.
922  {1,15,15,15,15,15},  // 2 .. x=0 y=0 z=2 D vertex only
923
924  {8,0,0,0,2,15},      // 3 .. x=0 y=1 z=0 D two coords
925  {9,0,0,15,15,15},    // 4 .. x=0 y=1 z=1 D one coord
926  {8,0,0,1,2,15},      // 5 .. x=0 y=1 z=2 D two coords.
927
928  {2,15,15,15,15,15},  // 6 .. x=0 y=2 z=0 D one vertex
929  {8,0,0,1,1,15},      // 7 .. x=0 y=2 z=1 D two coords
930  {3,15,15,15,15,15},  // 8 .. x=0 y=2 z=2 D one vertex
931
932  // the regions number <9,17>
933  {8,0,1,0,2,15},      // 9  .. x=1 y=0 z=0 D two coords
934  {9,0,1,15,15,15},    // 10 .. x=1 y=0 z=1 D one coord
935  {8,0,1,1,2,15},      // 11 .. x=1 y=0 z=2 D two coords
936
937  {9,0,2,15,15,15},    // 12 .. x=1 y=1 z=0 D one coord
938  {15,15,15,15,15,15}, // 13 .. x=1 y=1 z=1 inside the box, special case/value
939  {9,1,2,15,15,15},    // 14 .. x=1 y=1 z=2 D one corrd
940
941  {8,0,2,1,1,15},      // 15 .. x=1 y=2 z=0 D two coords
942  {9,1,1,15,15},       // 16 .. x=1 y=2 z=1 D one coord
943  {8,1,1,1,2,15},      // 17 .. x=1 y=2 z=2 D two coords
944
945  // the regions number <18,26>
946  {4,15,15,15,15,15},  // 18 .. x=2 y=0 z=0 D one vertex
947  {8,0,1,1,0,15},      // 19 .. x=2 y=0 z=1 D two coords
948  {5,15,15,15,15,15},  // 20 .. x=2 y=0 z=2 D one vertex
949
950  {8,0,2,1,0,15},      // 21 .. x=2 y=1 z=0 D two coords
951  {9,1,0,15,15,15},    // 22 .. x=2 y=1 z=1 D one coord
952  {8,1,0,1,2,15},      // 23 .. x=2 y=1 z=2 D two coords
953
954  {6,15,15,15,15,15},  // 24 .. x=2 y=2 z=0 D one vertex
955  {8,1,0,1,1,15},      // 25 .. x=2 y=2 z=1 D two coords
956  {7,15,15,15,15,15},  // 26 .. x=2 y=2 z=2 D one vertex
957};
958
959
960// The vertices that form all FARTHEST points with respect to the region
961// for all the regions possible, number of regions is 3^3 = 27,
962// since two parallel sides of bbox forms three disjoint spaces.
963// The vertices are given in anti-clockwise order, stopped by value 15,
964// at most 8 points, at least 1 point.
965// For testing, if the AABB is whole in the sphere, it is enough
966// to test only vertices, either 1,2,4, or 8.
967const int
968AxisAlignedBox3::fvertices[27][9] =
969{                   // region number.. position
970  {7,15,15,15,15,15,15,15,15}, // 0 .. x=0 y=0 z=0 D
971  {6,7,15,15,15,15,15,15,15},  // 1 .. x=0 y=0 z=1 D
972  {6,15,15,15,15,15,15,15,15}, // 2 .. x=0 y=0 z=2 D
973
974  {5,7,15,15,15,15,15,15,15},  // 3 .. x=0 y=1 z=0 D
975  {4,5,7,6,15,15,15,15,15},    // 4 .. x=0 y=1 z=1 D
976  {4,6,15,15,15,15,15,15,15},  // 5 .. x=0 y=1 z=2 D
977
978  {5,15,15,15,15,15,15,15,15}, // 6 .. x=0 y=2 z=0 D
979  {4,5,15,15,15,15,15,15,15},  // 7 .. x=0 y=2 z=1 D
980  {4,15,15,15,15,15,15,15,15}, // 8 .. x=0 y=2 z=2 D
981
982  // the regions number <9,17>
983  {3,7,15,15,15,15,15,15,15},  // 9  .. x=1 y=0 z=0 D
984  {7,3,2,6,15,15,15,15,15},    // 10 .. x=1 y=0 z=1 D
985  {2,6,15,15,15,15,15,15,15},  // 11 .. x=1 y=0 z=2 D
986
987  {5,1,3,7,15,15,15,15,15},    // 12 .. x=1 y=1 z=0 D
988  {0,7,1,6,3,4,5,2,15},        // 13 .. x=1 y=1 z=1 .. inside the box
989  {0,4,6,2,15,15,15,15,15},    // 14 .. x=1 y=1 z=2 D
990
991  {5,1,15,15,15,15,15,15,15},  // 15 .. x=1 y=2 z=0 D
992  {4,0,1,5,15,15,15,15,15},    // 16 .. x=1 y=2 z=1 D
993  {4,0,15,15,15,15,15,15,15},  // 17 .. x=1 y=2 z=2 D
994
995  // the regions number <18,26>
996  {3,15,15,15,15,15,15,15,15}, // 18 .. x=2 y=0 z=0 D
997  {2,3,15,15,15,15,15,15,15},  // 19 .. x=2 y=0 z=1 D
998  {2,15,15,15,15,15,15,15,15}, // 20 .. x=2 y=0 z=2 D
999
1000  {1,3,15,15,15,15,15,15,15},  // 21 .. x=2 y=1 z=0 D
1001  {1,0,2,3,15,15,15,15,15},    // 22 .. x=2 y=1 z=1 D
1002  {2,0,15,15,15,15,15,15,15},  // 23 .. x=2 y=1 z=2 D
1003
1004  {1,15,15,15,15,15,15,15,15}, // 24 .. x=2 y=2 z=0 D
1005  {0,1,15,15,15,15,15,15,15},  // 25 .. x=2 y=2 z=1 D
1006  {0,15,15,15,15,15,15,15,15}, // 26 .. x=2 y=2 z=2 D
1007};
1008
1009// Similar table as above, farthest points, but only the ones
1010// necessary for testing the intersection problem. If we do
1011// not consider the case 13, center of the sphere is inside the
1012// box, then we can always test at most 2 box vertices to say whether
1013// the whole box is inside the sphere.
1014// The number of vertices is minimized using some assumptions
1015// about the ortogonality of vertices and sphere properties.
1016const int
1017AxisAlignedBox3::fsvertices[27][9] =
1018{                   // region number.. position
1019  {7,15,15,15,15,15,15,15,15}, // 0 .. x=0 y=0 z=0 D 1 vertex
1020  {6,7,15,15,15,15,15,15,15},  // 1 .. x=0 y=0 z=1 D 2 vertices
1021  {6,15,15,15,15,15,15,15,15}, // 2 .. x=0 y=0 z=2 D 1 vertex
1022
1023  {5,7,15,15,15,15,15,15,15},  // 3 .. x=0 y=1 z=0 D 2 vertices
1024  {4,7,15,5,6,15,15,15,15},    // 4 .. x=0 y=1 z=1 D 4/2 vertices
1025  {4,6,15,15,15,15,15,15,15},  // 5 .. x=0 y=1 z=2 D 2 vertices
1026
1027  {5,15,15,15,15,15,15,15,15}, // 6 .. x=0 y=2 z=0 D 1 vertex
1028  {4,5,15,15,15,15,15,15,15},  // 7 .. x=0 y=2 z=1 D 2 vertices
1029  {4,15,15,15,15,15,15,15,15}, // 8 .. x=0 y=2 z=2 D 1 vertex
1030
1031  // the regions number <9,17>
1032  {3,7,15,15,15,15,15,15,15},  // 9  .. x=1 y=0 z=0 D 2 vertices
1033  {7,2,15,3,6,15,15,15,15},    // 10 .. x=1 y=0 z=1 D 4/2 vertices
1034  {2,6,15,15,15,15,15,15,15},  // 11 .. x=1 y=0 z=2 D 2 vertices
1035
1036  {5,3,15,1,7,15,15,15,15},    // 12 .. x=1 y=1 z=0 D 4/2 vertices
1037  {0,7,1,6,3,4,5,2,15},        // 13 .. x=1 y=1 z=1 .. inside the box
1038  {0,6,15,4,2,15,15,15,15},    // 14 .. x=1 y=1 z=2 D 4/2 vertices
1039
1040  {5,1,15,15,15,15,15,15,15},  // 15 .. x=1 y=2 z=0 D 2 vertices
1041  {4,1,15,0,5,15,15,15,15},    // 16 .. x=1 y=2 z=1 D 4/2 vertices
1042  {4,0,15,15,15,15,15,15,15},  // 17 .. x=1 y=2 z=2 D 2 vertices
1043
1044  // the regions number <18,26>
1045  {3,15,15,15,15,15,15,15,15}, // 18 .. x=2 y=0 z=0 D 1 vertex
1046  {2,3,15,15,15,15,15,15,15},  // 19 .. x=2 y=0 z=1 D 2 vertices
1047  {2,15,15,15,15,15,15,15,15}, // 20 .. x=2 y=0 z=2 D 1 vertex
1048
1049  {1,3,15,15,15,15,15,15,15},  // 21 .. x=2 y=1 z=0 D 2 vertices
1050  {1,2,15,0,3,15,15,15,15},    // 22 .. x=2 y=1 z=1 D 4/2 vertices
1051  {2,0,15,15,15,15,15,15,15},  // 23 .. x=2 y=1 z=2 D 2 vertices
1052
1053  {1,15,15,15,15,15,15,15,15}, // 24 .. x=2 y=2 z=0 D 1 vertex
1054  {0,1,15,15,15,15,15,15,15},  // 25 .. x=2 y=2 z=1 D 2 vertices
1055  {0,15,15,15,15,15,15,15,15}, // 26 .. x=2 y=2 z=2 D 1 vertex
1056};
1057
1058
1059// The fast computation of arctangent .. the maximal error is less
1060// than 4.1 degrees, according to Graphics GEMSII, 1991, pages 389--391
1061// Ron Capelli: "Fast approximation to the arctangent"
1062float
1063atan22(const float& y)
1064{
1065  const float x = 1.0;
1066  const float c = (float)(M_PI * 0.25);
1067 
1068  if (y < 0.0) {
1069    if (y < -1.0)
1070      return c * (-2.0f + x / y); // for angle in <-PI/2, -PI/4)
1071    else
1072      return c * (y / x); // for angle in <-PI/4 , 0>
1073  }
1074  else {
1075    if (y > 1.0)
1076      return c * (2.0f - x / y); // for angle in <PI/4, PI/2>
1077    else
1078      return c * (y / x); // for angle in <0, PI/2>
1079  }
1080}
1081
1082
1083float
1084AxisAlignedBox3::ProjectToSphereSA(const Vector3 &viewpoint, int *tcase) const
1085{
1086  int id = GetRegionID(viewpoint);
1087  *tcase = id;
1088
1089  // spherical projection .. SA represents solid angle
1090  if (id == 13) // .. inside the box
1091    return (float)(4.0*M_PI); // the whole sphere
1092  float SA = 0.0; // inital value
1093   
1094  int i = 0; // the pointer in the array of vertices
1095  while (bfaces[id][i] >= 0) {
1096    int axisO = bfaces[id][i++];
1097    int minvIdx = bfaces[id][i++];
1098    int maxvIdx = bfaces[id][i++];
1099    Vector3 vmin, vmax;
1100    GetVertex(minvIdx, vmin);
1101    GetVertex(maxvIdx, vmax);
1102    float h = fabs(vmin[axisO] - viewpoint[axisO]);
1103    int axis = (axisO + 1) % 3; // next axis
1104    float a = (vmin[axis] - viewpoint[axis]) / h; // minimum for v-range
1105    float b = (vmax[axis] - viewpoint[axis]) / h; // maximum for v-range
1106    //if (a > b) {
1107    //  FATAL << "ProjectToSphereSA::Error a > b\n";
1108    //  FATAL_ABORT;
1109    //}
1110    //if (vmin[axisO] != vmax[axisO]) {
1111    //  FATAL << "ProjectToSphereSA::Error a-axis != b-axis\n";
1112    //  FATAL_ABORT;     
1113    //}
1114    axis = (axisO + 2) % 3; // next second axis       
1115    float c = (vmin[axis] - viewpoint[axis]) / h; // minimum for u-range
1116    float d = (vmax[axis] - viewpoint[axis]) / h; // maximum for u-range
1117    //if (c > d) {
1118    //  FATAL << "ProjectToSphereSA::Error c > d\n";
1119    //  FATAL_ABORT;
1120    //}
1121    SA +=atan22(d*b/sqrt(b*b + d*d + 1.0f)) - atan22(b*c/sqrt(b*b + c*c + 1.0f))
1122      - atan22(d*a/sqrt(a*a + d*d + 1.0f)) + atan22(a*c/sqrt(a*a + c*c + 1.0f));
1123  }
1124
1125#if 0
1126  if ((SA > 2.0*M_PI) ||
1127      (SA < 0.0)) {
1128    FATAL << "The solid angle has strange value: ";
1129    FATAL << "SA = "<< SA << endl;
1130    FATAL_ABORT;
1131  }
1132#endif
1133 
1134  return SA;
1135}
1136
1137// Projects the box to a plane given a normal vector only and
1138// computes the surface area of the projected silhouette
1139// no clipping of the box is performed.
1140float
1141AxisAlignedBox3::ProjectToPlaneSA(const Vector3 &normal) const
1142{
1143  Vector3 size = Size();
1144
1145  // the surface area of the box to a yz-plane - perpendicular to x-axis
1146  float sax = size.y * size.z;
1147
1148  // the surface area of the box to a zx-plane - perpendicular to y-axis
1149  float say = size.z * size.x;
1150
1151  // the surface area of the box to a xy-plane - perpendicular to z-axis
1152  float saz = size.x * size.y;
1153 
1154  return sax * fabs(normal.x) + say * fabs(normal.y) + saz * fabs(normal.z);
1155}
1156
1157
1158
1159// This definition allows to be a point when answering true
1160bool
1161AxisAlignedBox3::IsCorrectAndNotPoint() const
1162{
1163  if ( (mMin.x > mMax.x) ||
1164       (mMin.y > mMax.y) ||
1165       (mMin.z > mMax.z) )
1166    return false; // box is not formed
1167
1168  if ( (mMin.x == mMax.x) &&
1169       (mMin.y == mMax.y) &&
1170       (mMin.z == mMax.z) )
1171    return false; // degenerates to a point
1172 
1173  return true;
1174}
1175
1176// This definition allows to be a point when answering true
1177bool
1178AxisAlignedBox3::IsPoint() const
1179{
1180  if ( (mMin.x == mMax.x) &&
1181       (mMin.y == mMax.y) &&
1182       (mMin.z == mMax.z) )
1183    return true; // degenerates to a point
1184 
1185  return false;
1186}
1187
1188// This definition requires shape of non-zero volume
1189bool
1190AxisAlignedBox3::IsSingularOrIncorrect() const
1191{
1192  if ( (mMin.x >= mMax.x) ||
1193       (mMin.y >= mMax.y) ||
1194       (mMin.z >= mMax.z) )
1195    return true; // box is not formed
1196
1197  return false; // has non-zero volume
1198}
1199
1200// returns true, when the sphere specified by the origin and radius
1201// fully contains the box
1202bool
1203AxisAlignedBox3::IsFullyContainedInSphere(const Vector3 &center, float rad) const
1204{
1205  int region = GetRegionID(center);
1206  float rad2 = rad*rad;
1207
1208  // vertex of the box
1209  Vector3 vertex;
1210
1211  int i = 0;
1212  for (i = 0 ; ; i++) {
1213    int a = fsvertices[region][i];
1214    if (a == 15)
1215      return true; // if was not false untill now, it must be contained
1216
1217    assert( (a>=0) && (a<8) );
1218   
1219    // normal vertex
1220    GetVertex(a, vertex);
1221
1222    if (SqrMagnitude(vertex - center) > rad2)
1223      return false;   
1224  } // for
1225 
1226}
1227
1228// returns true, when the volume of the sphere and volume of the
1229// axis aligned box has no intersection
1230bool
1231AxisAlignedBox3::HasNoIntersectionWithSphere(const Vector3 &center, float rad) const
1232{
1233  int region = GetRegionID(center);
1234  float rad2 = rad*rad;
1235
1236  // vertex of the box
1237  Vector3 vertex;
1238
1239  switch (csvertices[region][0]) {
1240  case 8: {
1241    // test two coordinates described within the field
1242    int face = 3*csvertices[region][1] + csvertices[region][2];
1243    float dist = GetExtent(face) - center[csvertices[region][2]];
1244    dist *= dist;
1245    face = 3 * (csvertices[region][3]) + csvertices[region][4];
1246    float dist2 = GetExtent(face) - center[csvertices[region][4]];
1247    dist += (dist2 * dist2);
1248    if (dist > rad2)
1249      return true; // no intersection is possible
1250  }
1251  case 9: {
1252    // test one coordinate described within the field
1253    int face = 3*csvertices[region][1] + csvertices[region][2];
1254    float dist = fabs(GetExtent(face) - center[csvertices[region][2]]);
1255    if (dist > rad)
1256      return true; // no intersection is possible
1257  }
1258  case 15:
1259    return false; // box and sphere surely has intersection
1260  default: {
1261    // test using normal vertices
1262    assert( (csvertices[region][0]>=0) && (csvertices[region][0]<8) );
1263
1264    // normal vertex
1265    GetVertex(csvertices[region][0], vertex);
1266
1267    if (SqrMagnitude(vertex - center) > rad2)
1268      return true; // no intersectino is possible   
1269    }
1270  } // switch
1271
1272  return false; // partial or full containtment
1273}
1274
1275#if 0
1276// Given the sphere, determine the mutual position between the
1277// sphere and box
1278
1279// SOME BUG IS INSIDE !!!! V.H. 25/4/2001
1280int
1281AxisAlignedBox3::MutualPositionWithSphere(const Vector3 &center, float rad) const
1282{
1283  int region = GetRegionID(center);
1284  float rad2 = rad*rad;
1285
1286  // vertex of the box
1287  Vector3 vertex;
1288
1289  // first testing for  full containtment - whether sphere fully
1290  // contains the box
1291  int countInside = 0; // how many points were found inside
1292 
1293  int i = 0;
1294  for (i = 0 ; ; i++) {
1295    int a = fsvertices[region][i];
1296    if (a == 15)
1297      return 1; // the sphere fully contain the box
1298
1299    assert( (a>=0) && (a<8) );
1300   
1301    // normal vertex
1302    GetVertex(a, vertex);
1303
1304    if (SqrMagnitude(vertex - center) <= rad2)
1305      countInside++; // the number of vertices inside the sphere
1306    else {
1307      if (countInside)
1308        return 0; // partiall overlap has been found
1309      // the sphere does not fully contain the box .. only way to break
1310      // this loop and go for other testing
1311      break;
1312    }
1313  } // for
1314
1315  // now only box and sphere can partially overlap or no intersection
1316  switch (csvertices[region][0]) {
1317  case 8: {
1318    // test two coordinates described within the field
1319    int face = 3*csvertices[region][1] + csvertices[region][2];
1320    float dist = GetExtent(face) - center[csvertices[region][2]];
1321    dist *= dist;
1322    face = 3 * (csvertices[region][3]) + csvertices[region][4];
1323    float dist2 = GetExtent(face) - center[csvertices[region][4]];
1324    dist += (dist2 * dist2);
1325    if (dist > rad2 )
1326      return -1; // no intersection is possible
1327  }
1328  case 9: {
1329    // test one coordinate described within the field
1330    int face = 3*csvertices[region][1] + csvertices[region][2];
1331    float dist = fabs(GetExtent(face) - center[csvertices[region][2]]);
1332    if (dist > rad)
1333      return -1; // no intersection is possible
1334  }
1335  case 15:
1336    return 0 ; // partial overlap is now guaranteed
1337  default: {
1338    // test using normal vertices
1339    assert( (csvertices[region][0]>=0) && (csvertices[region][0]<8) );
1340   
1341    // normal vertex
1342    GetVertex(csvertices[region][0], vertex);
1343
1344    if (SqrMagnitude(vertex - center) > rad2)
1345      return -1; // no intersection is possible   
1346    }
1347  } // switch
1348
1349  return 0; // partial intersection is guaranteed
1350}
1351#else
1352
1353// Some maybe smarter version, extendible easily to d-dimensional
1354// space!
1355// Given a sphere described by the center and radius,
1356// the fullowing function returns:
1357//   -1 ... the sphere and the box are completely separate
1358//    0 ... the sphere and the box only partially overlap
1359//    1 ... the sphere contains fully the box
1360//  Note: the case when box fully contains the sphere is not reported
1361//        since it was not required.
1362int
1363AxisAlignedBox3::MutualPositionWithSphere(const Vector3 &center, float rad) const
1364{
1365//#define SPEED_UP
1366
1367#ifndef SPEED_UP
1368  // slow version, instructively written 
1369#if 0
1370  // does it make sense to test
1371  // checking the sides of the box for possible non-intersection
1372  if ( ((center.x + rad) < mMin.x) ||
1373       ((center.x - rad) > mMax.x) ||
1374       ((center.y + rad) < mMin.y) ||
1375       ((center.y - rad) > mMax.y) ||
1376       ((center.z + rad) < mMin.z) ||
1377       ((center.z - rad) > mMax.z) ) {
1378    // cout << "r ";
1379    return -1;  // no overlap is possible
1380  }
1381#endif
1382 
1383  // someoverlap is possible, check the distance of vertices
1384  rad = rad*rad;
1385  float sumMin = 0;
1386  // Try to minimize the function of a distance
1387  // from the sphere center
1388
1389  // for x-axis
1390  float minSqrX = sqr(mMin.x - center.x);
1391  float maxSqrX = sqr(mMax.x - center.x);
1392  if (center.x < mMin.x)
1393    sumMin = minSqrX;
1394  else
1395    if (center.x > mMax.x)
1396      sumMin = maxSqrX;
1397
1398  // for y-axis
1399  float minSqrY = sqr(mMin.y - center.y);
1400  float maxSqrY = sqr(mMax.y - center.y);
1401  if (center.y < mMin.y)
1402    sumMin += minSqrY;
1403  else
1404    if (center.y > mMax.y)
1405      sumMin += maxSqrY;
1406
1407  // for z-axis
1408  float minSqrZ = sqr(mMin.z - center.z);
1409  float maxSqrZ = sqr(mMax.z - center.z);
1410  if (center.z < mMin.z)
1411    sumMin += minSqrZ;
1412  else
1413    if (center.z > mMax.z)
1414      sumMin += maxSqrZ;
1415
1416  if (sumMin > rad)
1417    return -1; // no intersection between sphere and box
1418
1419  // try to find out the maximum distance between the
1420  // sphere center and vertices
1421  float sumMax = 0;
1422 
1423  if (minSqrX > maxSqrX)
1424    sumMax = minSqrX;
1425  else
1426    sumMax = maxSqrX;
1427 
1428  if (minSqrY > maxSqrY)
1429    sumMax += minSqrY;
1430  else
1431    sumMax += maxSqrY;
1432 
1433  if (minSqrZ > maxSqrZ)
1434    sumMax += minSqrZ;
1435  else
1436    sumMax += maxSqrZ;
1437 
1438  // sumMin < rad
1439  if (sumMax < rad)
1440    return 1; // the sphere contains the box completely
1441
1442  // partial intersection, part of the box is outside the sphere
1443  return 0;
1444#else
1445
1446  // Optimized version of the test
1447 
1448#ifndef __VECTOR_HACK
1449#error "__VECTOR_HACK for Vector3 was not defined"
1450#endif
1451
1452  // some overlap is possible, check the distance of vertices
1453  rad = rad*rad;
1454  float sumMin = 0;
1455  float sumMax = 0;
1456  // Try to minimize the function of a distance
1457  // from the sphere center
1458
1459  const float *minp = &(min[0]);
1460  const float *maxp = &(max[0]);
1461  const float *pcenter = &(center[0]);
1462 
1463  // for x-axis
1464  for (int i = 0; i < 3; i++, minp++, maxp++, pcenter++) {
1465    float minsqr = sqr(*minp - *pcenter);
1466    float maxsqr = sqr(*maxp - *pcenter);
1467    if (*pcenter < *minp)
1468      sumMin += minsqr;
1469    else
1470      if (*pcenter > *maxp)
1471        sumMin += maxsqr;
1472    sumMax += (minsqr > maxsqr) ? minsqr : maxsqr;
1473  }
1474
1475  if (sumMin > rad)
1476    return -1; // no intersection between sphere and box
1477
1478  // sumMin < rad
1479  if (sumMax < rad)
1480    return 1; // the sphere contains the box completely
1481
1482  // partial intersection, part of the box is outside the sphere
1483  return 0;
1484#endif
1485}
1486#endif
1487
1488// Given the cube describe by the center and the half-sie,
1489// determine the mutual position between the cube and the box
1490int
1491AxisAlignedBox3::MutualPositionWithCube(const Vector3 &center, float radius) const
1492{
1493  // the cube is described by the center and the distance to the any face
1494  // along the axes
1495
1496  // Note on efficiency!
1497  // Can be quite optimized using tables, but I do not have time
1498  // V.H. 18/11/2001
1499
1500  AxisAlignedBox3 a =
1501    AxisAlignedBox3(Vector3(center.x - radius, center.y - radius, center.z - radius),
1502           Vector3(center.x + radius, center.y + radius, center.z + radius));
1503
1504  if (a.Includes(*this))
1505    return 1; // cube contains the box
1506
1507  if (OverlapS(a,*this))
1508    return 0; // cube partially overlap the box
1509
1510  return -1; // completely separate 
1511}
1512
1513void
1514AxisAlignedBox3::GetSqrDistances(const Vector3 &point,
1515                                 float &minDistance,
1516                                 float &maxDistance
1517                                 ) const
1518{
1519
1520 
1521#ifndef __VECTOR_HACK
1522#error "__VECTOR_HACK for Vector3 was not defined"
1523#endif
1524
1525  // some overlap is possible, check the distance of vertices
1526  float sumMin = 0;
1527  float sumMax = 0;
1528
1529  // Try to minimize the function of a distance
1530  // from the sphere center
1531
1532  const float *minp = &(mMin[0]);
1533  const float *maxp = &(mMax[0]);
1534  const float *pcenter = &(point[0]);
1535 
1536  // for x-axis
1537  for (int i = 0; i < 3; i++, minp++, maxp++, pcenter++) {
1538    float minsqr = sqr(*minp - *pcenter);
1539    float maxsqr = sqr(*maxp - *pcenter);
1540    if (*pcenter < *minp)
1541      sumMin += minsqr;
1542    else
1543      if (*pcenter > *maxp)
1544        sumMin += maxsqr;
1545    sumMax += (minsqr > maxsqr) ? minsqr : maxsqr;
1546  }
1547
1548  minDistance = sumMin;
1549  maxDistance = sumMax;
1550}
1551
1552
1553int
1554AxisAlignedBox3::Side(const Plane3 &plane) const
1555{
1556  Vector3 v;
1557  int i, m=3, M=-3, s;
1558 
1559  for (i=0;i<8;i++) {
1560    GetVertex(i, v);
1561    if((s = plane.Side(v)) < m)
1562      m=s;
1563    if(s > M)
1564      M=s;
1565    if (m && m==-M)
1566      return 0;
1567  }
1568 
1569  return (m == M) ? m : m + M;
1570}
1571
1572
1573Plane3 AxisAlignedBox3::GetPlane(const int face) const
1574{
1575        switch (face)
1576        {
1577       
1578        case 0:
1579                return Plane3(Vector3(-1, 0, 0), mMin);
1580        case 1:
1581                return Plane3(Vector3(1, 0, 0), mMax);
1582        case 2:
1583                return Plane3(Vector3(0, -1, 0), mMin);
1584        case 3:
1585                return Plane3(Vector3(0, 1, 0), mMax);
1586        case 4:
1587                return Plane3(Vector3(0, 0, -1), mMin);
1588        case 5:
1589                return Plane3(Vector3(0, 0, 1), mMax);
1590        }
1591
1592        // should not come here
1593        return Plane3();
1594}
1595
1596
1597int
1598AxisAlignedBox3::GetFaceVisibilityMask(const Rectangle3 &rectangle) const
1599{
1600  int mask = 0;
1601  for (int i=0; i < 4; i++)
1602    mask |= GetFaceVisibilityMask(rectangle.mVertices[i]);
1603  return mask;
1604}
1605
1606
1607int
1608AxisAlignedBox3::GetFaceVisibilityMask(const Vector3 &position) const {
1609 
1610  // assume that we are not inside the box
1611  int c=0;
1612 
1613  if (position.x<(mMin.x-Limits::Small))
1614    c|=1;
1615  else
1616    if (position.x>(mMax.x+Limits::Small))
1617      c|=2;
1618 
1619  if (position.y<(mMin.y-Limits::Small))
1620    c|=4;
1621  else
1622    if (position.y>(mMax.y+Limits::Small))
1623      c|=8;
1624 
1625  if (position.z<(mMin.z-Limits::Small))
1626    c|=16;
1627  else
1628    if (position.z>(mMax.z+Limits::Small))
1629      c|=32;
1630 
1631  return c;
1632}
1633
1634
1635Rectangle3
1636AxisAlignedBox3::GetFace(const int face) const
1637{
1638  Vector3 v[4];
1639  switch (face) {
1640   
1641  case 0:
1642    v[3].SetValue(mMin.x,mMin.y,mMin.z);
1643    v[2].SetValue(mMin.x,mMax.y,mMin.z);
1644    v[1].SetValue(mMin.x,mMax.y,mMax.z);
1645    v[0].SetValue(mMin.x,mMin.y,mMax.z);
1646    break;
1647   
1648  case 1:
1649    v[0].SetValue(mMax.x,mMin.y,mMin.z);
1650    v[1].SetValue(mMax.x,mMax.y,mMin.z);
1651    v[2].SetValue(mMax.x,mMax.y,mMax.z);
1652    v[3].SetValue(mMax.x,mMin.y,mMax.z);
1653    break;
1654   
1655  case 2:
1656    v[3].SetValue(mMin.x,mMin.y,mMin.z);
1657    v[2].SetValue(mMin.x,mMin.y,mMax.z);
1658    v[1].SetValue(mMax.x,mMin.y,mMax.z);
1659    v[0].SetValue(mMax.x,mMin.y,mMin.z);
1660    break;
1661 
1662  case 3:
1663    v[0].SetValue(mMin.x,mMax.y,mMin.z);
1664    v[1].SetValue(mMin.x,mMax.y,mMax.z);
1665    v[2].SetValue(mMax.x,mMax.y,mMax.z);
1666    v[3].SetValue(mMax.x,mMax.y,mMin.z);
1667    break;
1668
1669  case 4:
1670    v[3].SetValue(mMin.x,mMin.y,mMin.z);
1671    v[2].SetValue(mMax.x,mMin.y,mMin.z);
1672    v[1].SetValue(mMax.x,mMax.y,mMin.z);
1673    v[0].SetValue(mMin.x,mMax.y,mMin.z);
1674    break;
1675 
1676  case 5:
1677    v[0].SetValue(mMin.x,mMin.y,mMax.z);
1678    v[1].SetValue(mMax.x,mMin.y,mMax.z);
1679    v[2].SetValue(mMax.x,mMax.y,mMax.z);
1680    v[3].SetValue(mMin.x,mMax.y,mMax.z);
1681    break;
1682  }
1683 
1684  return Rectangle3(v[0], v[1], v[2], v[3]);
1685}
1686
1687struct VertexData
1688{
1689        Vector3 mVertex;
1690        float mAngle;
1691       
1692        VertexData(Vector3 vtx, float angle): mVertex(vtx), mAngle(angle)
1693        {}
1694
1695        bool operator<(const VertexData &b) const
1696        {
1697                return mAngle > b.mAngle;
1698        }
1699};
1700
1701// TODO: use a table to avoid normal and distance computations
1702Polygon3 *AxisAlignedBox3::CrossSection(const Plane3 &plane) const
1703{
1704        Polygon3 *planePoly = new Polygon3();
1705       
1706        int side[8];
1707        bool onFrontSide = false, onBackSide = false;
1708       
1709        Vector3 vtx;
1710        //-- compute classification of vertices
1711        for (int i = 0; i < 8; ++i)
1712        {
1713                GetVertex(i, vtx);
1714                side[i] = plane.Side(vtx);
1715                if (side[i] > 0)
1716                        onFrontSide = true;
1717                else if (side[i] < 0)
1718                        onBackSide = true;
1719                else // vertex coincident => push_back
1720                        planePoly->mVertices.push_back(vtx);
1721        }
1722
1723        //-- find intersections
1724        if (onFrontSide && onBackSide)
1725        {
1726                Vector3 ptA, ptB;
1727                for (int i = 0; i < 12; ++ i)
1728                {
1729                        int aIdx, bIdx;
1730                        GetEdge(i, aIdx, bIdx);
1731
1732                        ptA = GetVertex(aIdx);
1733                        ptB = GetVertex(bIdx);
1734
1735                        int sideA = side[aIdx];
1736                        int sideB = side[bIdx];
1737
1738                        if (((sideA > 0) && (sideB < 0)) || (sideA < 0) && (sideB > 0))
1739                                planePoly->mVertices.push_back(plane.FindIntersection(ptA, ptB));       
1740                }
1741        }
1742
1743        // order intersections 
1744        if (planePoly->mVertices.size() > 3)
1745        {
1746                Vector3 centerOfMass(0);
1747                int i;
1748                // compute center of mass
1749                for (i = 0; i < (int)planePoly->mVertices.size(); ++ i)
1750                        centerOfMass += planePoly->mVertices[i];
1751               
1752                centerOfMass /= (float)planePoly->mVertices.size();
1753
1754                vector<VertexData> vertexData;
1755               
1756                Vector3 refVec = Normalize(centerOfMass - planePoly->mVertices[0]);
1757
1758                // compute angle to reference point
1759                for (i = 1; i < (int)planePoly->mVertices.size(); ++ i)
1760                  {
1761                    float angle =
1762                      Angle(refVec, centerOfMass - planePoly->mVertices[i], plane.mNormal);
1763                   
1764                    vertexData.push_back(VertexData(planePoly->mVertices[i], angle));
1765                  }
1766               
1767                std::stable_sort(vertexData.begin(), vertexData.end());
1768
1769                // update vertices
1770                for (i = 1; i < (int)planePoly->mVertices.size(); ++ i)
1771                        planePoly->mVertices[i] = vertexData[i - 1].mVertex;
1772        }
1773        else if (planePoly->mVertices.size() == 3)
1774        {
1775                // fix orientation if needed
1776                if (DotProd(planePoly->GetNormal(), plane.mNormal) < 0)
1777                {
1778                        Vector3 v = planePoly->mVertices[1];
1779                        planePoly->mVertices[1] = planePoly->mVertices[2];
1780                        planePoly->mVertices[2] = v;
1781                }
1782        }
1783       
1784        return planePoly;
1785}
1786
1787
1788bool AxisAlignedBox3::GetRaySegment(const Ray &ray,
1789                                                                        float &minT,
1790                                                                        float &maxT) const
1791{
1792        maxT = 1e15f;
1793        minT = 0;
1794       
1795        // test with bounding box
1796        if (!GetMinMaxT(ray, &minT, &maxT))
1797                return false;
1798
1799        if (minT < 0)
1800                minT = 0;
1801
1802        // bound ray or line segment
1803        if (//(ray.GetType() == Ray::LOCAL_RAY) &&
1804            !ray.intersections.empty() &&
1805                (ray.intersections[0].mT <= maxT))
1806        {
1807                maxT = ray.intersections[0].mT;
1808        }
1809
1810        return true;
1811}
1812
1813
1814  // Compute tmin and tmax for a ray, whenever required .. need not pierce box
1815int
1816AxisAlignedBox3::ComputeMinMaxT(const Vector3 &origin,
1817                                                                const Vector3 &direction,
1818                                                                float *tmin,
1819                                                                float *tmax) const
1820{
1821
1822  register float minx, maxx;
1823
1824 
1825  Vector3 invDirection;
1826  const float eps = 1e-6f;
1827  const float invEps = 1e6f;
1828 
1829  // it does change the ray direction very slightly,
1830  // but the size direction vector is not practically changed
1831 
1832  if (fabs(direction.x) < eps) {
1833    if (direction.x < 0.0f)
1834      invDirection.x = -invEps;
1835    else
1836      invDirection.x = invEps;
1837  }
1838  else
1839    invDirection.x = 1.0f / direction.x;
1840 
1841  if (fabs(direction.y) < eps) {
1842    if (direction.y < 0.0)
1843      invDirection.y = -invEps;
1844    else
1845      invDirection.y = invEps;
1846  }
1847  else
1848    invDirection.y = 1.0f / direction.y;
1849 
1850  if (fabs(direction.z) < eps) {
1851    if (direction.z < 0.0f)
1852      invDirection.z = -invEps;
1853    else
1854      invDirection.z = invEps;
1855  }
1856  else
1857    invDirection.z = 1.0f / direction.z;
1858
1859
1860 
1861  if (fabs(direction.x) < 0.001f) {
1862    if (mMin.x < origin.x && mMax.x > origin.x) {
1863      minx = -MAXFLOAT;
1864      maxx = MAXFLOAT;
1865    }
1866    else
1867      return 0;
1868  }
1869  else {
1870    float t1 = (mMin.x - origin.x) * invDirection.x;
1871    float t2 = (mMax.x - origin.x) * invDirection.x;
1872    if (t1 < t2) {
1873      minx = t1;
1874      maxx = t2;
1875    }
1876    else {
1877      minx = t2;
1878      maxx = t1;
1879    }
1880        //    if (maxx < 0.0)
1881        //      return 0;
1882  }
1883
1884  *tmin = minx;
1885  *tmax = maxx;
1886 
1887  if (fabs(direction.y) < 0.001) {
1888    if (mMin.y < origin.y && mMax.y > origin.y) {
1889      minx = -MAXFLOAT;
1890      maxx = MAXFLOAT;
1891    }
1892    else
1893      return 0;
1894  }
1895  else {
1896    float t1 = (mMin.y - origin.y) * invDirection.y;
1897    float t2 = (mMax.y - origin.y) * invDirection.y;
1898    if (t1 < t2) {
1899      minx = t1;
1900      maxx = t2;
1901    }
1902    else {
1903      minx = t2;
1904      maxx = t1;
1905    }
1906        //    if (maxx < 0.0)
1907        //      return 0;
1908  }
1909
1910  if (minx > *tmin)
1911    *tmin = minx;
1912  if (maxx < *tmax)
1913    *tmax = maxx;
1914 
1915  if (fabs(direction.z) < 0.001) {
1916    if (mMin.z < origin.z && mMax.z > origin.z) {
1917      minx = -MAXFLOAT;
1918      maxx = MAXFLOAT;
1919    }
1920    else
1921      return 0;
1922  }
1923  else {
1924    float t1 = (mMin.z - origin.z) * invDirection.z;
1925    float t2 = (mMax.z - origin.z) * invDirection.z;
1926    if (t1 < t2) {
1927      minx = t1;
1928      maxx = t2;
1929    }
1930    else {
1931      minx = t2;
1932      maxx = t1;
1933    }
1934        //    if (maxx < 0.0)
1935        //      return 0;
1936  }
1937
1938  if (minx > *tmin)
1939    *tmin = minx;
1940  if (maxx < *tmax)
1941    *tmax = maxx;
1942
1943  return 1; // yes, intersection was found
1944}
1945
1946
1947bool AxisAlignedBox3::GetIntersectionFace(Rectangle3 &face,
1948                                                                                  const AxisAlignedBox3 &neighbour) const
1949                                                                                   
1950{
1951        if (EpsilonEqual(mMin[0], neighbour.Max(0)))
1952        {
1953                float maxy = min(mMax.y, neighbour.mMax.y);
1954                float maxz = min(mMax.z, neighbour.mMax.z);
1955                float miny = max(mMin.y, neighbour.mMin.y);
1956                float minz = max(mMin.z, neighbour.mMin.z);
1957
1958                face.mVertices[3].SetValue(mMin.x, miny, minz);
1959                face.mVertices[2].SetValue(mMin.x, maxy, minz);
1960                face.mVertices[1].SetValue(mMin.x, maxy, maxz);
1961                face.mVertices[0].SetValue(mMin.x, miny, maxz);
1962
1963        return true;
1964        }
1965    if (EpsilonEqual(mMax[0], neighbour.Min(0)))
1966        {
1967                float maxy = min(mMax.y, neighbour.mMax.y);
1968                float maxz = min(mMax.z, neighbour.mMax.z);
1969                float miny = max(mMin.y, neighbour.mMin.y);
1970                float minz = max(mMin.z, neighbour.mMin.z);
1971
1972                face.mVertices[0].SetValue(mMax.x, miny, minz);
1973                face.mVertices[1].SetValue(mMax.x, maxy, minz);
1974                face.mVertices[2].SetValue(mMax.x, maxy, maxz);
1975                face.mVertices[3].SetValue(mMax.x, miny, maxz);
1976
1977                return true;
1978        }
1979    if (EpsilonEqual(mMin[1], neighbour.Max(1)))
1980        {
1981                float maxx = min(mMax.x, neighbour.mMax.x);
1982                float maxz = min(mMax.z, neighbour.mMax.z);
1983                float minx = max(mMin.x, neighbour.mMin.x);
1984                float minz = max(mMin.z, neighbour.mMin.z);
1985
1986                face.mVertices[3].SetValue(minx, mMin.y, minz);
1987                face.mVertices[2].SetValue(minx, mMin.y, maxz);
1988                face.mVertices[1].SetValue(maxx, mMin.y, maxz);
1989                face.mVertices[0].SetValue(maxx, mMin.y, minz);
1990               
1991                return true;
1992        }
1993        if (EpsilonEqual(mMax[1], neighbour.Min(1)))
1994        {               
1995                float maxx = min(mMax.x, neighbour.mMax.x);
1996                float maxz = min(mMax.z, neighbour.mMax.z);
1997                float minx = max(mMin.x, neighbour.mMin.x);
1998                float minz = max(mMin.z, neighbour.mMin.z);
1999
2000                face.mVertices[0].SetValue(minx, mMax.y, minz);
2001                face.mVertices[1].SetValue(minx, mMax.y, maxz);
2002                face.mVertices[2].SetValue(maxx, mMax.y, maxz);
2003                face.mVertices[3].SetValue(maxx, mMax.y, minz);
2004
2005                return true;
2006        }
2007        if (EpsilonEqual(mMin[2], neighbour.Max(2)))
2008        {
2009                float maxx = min(mMax.x, neighbour.mMax.x);
2010                float maxy = min(mMax.y, neighbour.mMax.y);
2011                float minx = max(mMin.x, neighbour.mMin.x);
2012                float miny = max(mMin.y, neighbour.mMin.y);
2013
2014                face.mVertices[3].SetValue(minx, miny, mMin.z);
2015                face.mVertices[2].SetValue(maxx, miny, mMin.z);
2016                face.mVertices[1].SetValue(maxx, maxy, mMin.z);
2017                face.mVertices[0].SetValue(minx, maxy, mMin.z);
2018       
2019                return true;
2020        }
2021        if (EpsilonEqual(mMax[2], neighbour.Min(2)))
2022        {
2023                float maxx = min(mMax.x, neighbour.mMax.x);
2024                float maxy = min(mMax.y, neighbour.mMax.y);
2025                float minx = max(mMin.x, neighbour.mMin.x);
2026                float miny = max(mMin.y, neighbour.mMin.y);
2027
2028                face.mVertices[0].SetValue(minx, miny, mMax.z);
2029                face.mVertices[1].SetValue(maxx, miny, mMax.z);
2030                face.mVertices[2].SetValue(maxx, maxy, mMax.z);
2031                face.mVertices[3].SetValue(minx, maxy, mMax.z);
2032
2033                return true;
2034        }
2035
2036        return false;
2037}
2038
2039
2040void AxisAlignedBox3::AddBoxToMesh(Mesh *mesh) const
2041{
2042        // add 6 vertices of the box
2043        int index = (int)mesh->mVertices.size();
2044       
2045        for (int i=0; i < 8; i++)
2046        {
2047                Vector3 v;
2048                GetVertex(i, v);
2049                mesh->mVertices.push_back(v);
2050        }
2051       
2052        mesh->AddFace(new Face(index + 0, index + 1, index + 3, index + 2) );
2053        mesh->AddFace(new Face(index + 0, index + 2, index + 6, index + 4) );
2054        mesh->AddFace(new Face(index + 4, index + 6, index + 7, index + 5) );
2055       
2056        mesh->AddFace(new Face(index + 3, index + 1, index + 5, index + 7) );
2057        mesh->AddFace(new Face(index + 0, index + 4, index + 5, index + 1) );
2058        mesh->AddFace(new Face(index + 2, index + 3, index + 7, index + 6) );
2059}
2060
2061
2062void AxisAlignedBox3::ExtractPolys(PolygonContainer &polys) const
2063{
2064        Polygon3 *face1 = new Polygon3();
2065        polys.push_back(face1);   
2066
2067    face1->mVertices.push_back(Vector3(mMin.x, mMin.y, mMax.z));
2068        face1->mVertices.push_back(Vector3(mMin.x, mMax.y, mMax.z));
2069        face1->mVertices.push_back(Vector3(mMin.x, mMax.y ,mMin.z));
2070        face1->mVertices.push_back(Vector3(mMin.x, mMin.y, mMin.z));
2071
2072        Polygon3 *face2 = new Polygon3(); 
2073        polys.push_back(face2);
2074   
2075    face2->mVertices.push_back(Vector3(mMax.x, mMin.y, mMin.z));
2076    face2->mVertices.push_back(Vector3(mMax.x, mMax.y, mMin.z));
2077    face2->mVertices.push_back(Vector3(mMax.x, mMax.y, mMax.z));
2078    face2->mVertices.push_back(Vector3(mMax.x, mMin.y, mMax.z));
2079 
2080        Polygon3 *face3 = new Polygon3(); 
2081        polys.push_back(face3);
2082
2083    face3->mVertices.push_back(Vector3(mMax.x, mMin.y ,mMin.z));
2084        face3->mVertices.push_back(Vector3(mMax.x, mMin.y, mMax.z));
2085        face3->mVertices.push_back(Vector3(mMin.x, mMin.y, mMax.z));
2086        face3->mVertices.push_back(Vector3(mMin.x, mMin.y, mMin.z));
2087
2088        Polygon3 *face4 = new Polygon3(); 
2089        polys.push_back(face4);
2090
2091        face4->mVertices.push_back(Vector3(mMin.x, mMax.y, mMin.z));
2092        face4->mVertices.push_back(Vector3(mMin.x, mMax.y, mMax.z));
2093        face4->mVertices.push_back(Vector3(mMax.x, mMax.y, mMax.z));
2094        face4->mVertices.push_back(Vector3(mMax.x, mMax.y, mMin.z));
2095   
2096        Polygon3 *face5 = new Polygon3(); 
2097        polys.push_back(face5);
2098
2099        face5->mVertices.push_back(Vector3(mMin.x, mMax.y, mMin.z));
2100    face5->mVertices.push_back(Vector3(mMax.x, mMax.y, mMin.z));
2101    face5->mVertices.push_back(Vector3(mMax.x, mMin.y, mMin.z));
2102        face5->mVertices.push_back(Vector3(mMin.x, mMin.y, mMin.z));
2103
2104        Polygon3 *face6 = new Polygon3(); 
2105        polys.push_back(face6); 
2106 
2107    face6->mVertices.push_back(Vector3(mMin.x, mMin.y, mMax.z));
2108    face6->mVertices.push_back(Vector3(mMax.x, mMin.y, mMax.z));
2109    face6->mVertices.push_back(Vector3(mMax.x, mMax.y, mMax.z));
2110    face6->mVertices.push_back(Vector3(mMin.x, mMax.y, mMax.z));
2111
2112}
Note: See TracBrowser for help on using the repository browser.