source: trunk/VUT/GtpVisibilityPreprocessor/src/AxisAlignedBox3.cpp @ 448

Revision 448, 48.5 KB checked in by mattausch, 19 years ago (diff)

fixed bug in VspBspTree?
view cells in VssPreprocessor?
bounding rays for vspkdtree

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