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

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