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

Revision 1951, 61.6 KB checked in by mattausch, 18 years ago (diff)

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