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

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