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

Revision 1047, 58.0 KB checked in by mattausch, 19 years ago (diff)

fixed vsp part of vsp osp tree
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
11namespace GtpVisibilityPreprocessor {
12
13#define FATAL Debug
14#define FATAL_ABORT exit(1)
15
16
17// AxisAlignedBox3 implementations
18
19// Overload << operator for C++-style output
20ostream&
21operator<< (ostream &s, const AxisAlignedBox3 &A)
22{
23  return s << '[' << A.mMin.x << ", " << A.mMin.y << ", " << A.mMin.z << "]["
24           << A.mMax.x << ", " << A.mMax.y << ", " << A.mMax.z << ']';
25}
26
27// Overload >> operator for C++-style input
28istream&
29operator>> (istream &s, AxisAlignedBox3 &A)
30{
31  char a;
32  // read "[min.x, min.y, min.z][mMax.x, mMax.y, mMax.z]"
33  return s >> a >> A.mMin.x >> a >> A.mMin.y >> a >> A.mMin.z >> a >> a
34    >> A.mMax.x >> a >> A.mMax.y >> a >> A.mMax.z >> a;
35}
36
37bool
38AxisAlignedBox3::Unbounded() const
39{
40  return (mMin == Vector3(-MAXFLOAT)) ||
41         (mMax == Vector3(-MAXFLOAT));
42}
43
44void
45AxisAlignedBox3::Include(const Vector3 &newpt)
46{
47  Minimize(mMin, newpt);
48  Maximize(mMax, newpt);
49}
50
51void
52AxisAlignedBox3::Include(const Polygon3 &newpoly)
53{
54        VertexContainer::const_iterator it, it_end = newpoly.mVertices.end();
55
56        for (it = newpoly.mVertices.begin(); it != it_end; ++ it)
57                Include(*it);
58}
59
60void
61AxisAlignedBox3::Include(const PolygonContainer &polys)
62{
63        PolygonContainer::const_iterator it, it_end = polys.end();
64
65        for (it = polys.begin(); it != it_end; ++ it)
66                Include(*(*it));
67}
68
69void
70AxisAlignedBox3::Include(Mesh *mesh)
71{
72        VertexContainer::const_iterator it, it_end = mesh->mVertices.end();
73
74        for (it = mesh->mVertices.begin(); it != it_end; ++ it)
75                Include(*it);
76}
77
78
79void
80AxisAlignedBox3::Include(const AxisAlignedBox3 &bbox)
81{
82  Minimize(mMin, bbox.mMin);
83  Maximize(mMax, bbox.mMax);
84}
85
86bool
87AxisAlignedBox3::IsCorrect()
88{
89  if ( (mMin.x > mMax.x) ||
90       (mMin.y > mMax.y) ||
91       (mMin.z > mMax.z) )
92    return false; // box is not formed
93  return true;
94}
95
96void
97AxisAlignedBox3::GetEdge(const int edge, Vector3 *a, Vector3 *b) const
98{
99  switch(edge) {
100  case 0:
101    a->SetValue(mMin.x, mMin.y, mMin.z);
102    b->SetValue(mMin.x, mMin.y, mMax.z);
103    break;
104  case 1:
105    a->SetValue(mMin.x, mMin.y, mMin.z);
106    b->SetValue(mMin.x, mMax.y, mMin.z);
107    break;
108  case 2:
109    a->SetValue(mMin.x, mMin.y, mMin.z);
110    b->SetValue(mMax.x, mMin.y, mMin.z);
111    break;
112  case 3:
113    a->SetValue(mMax.x, mMax.y, mMax.z);
114    b->SetValue(mMax.x, mMax.y, mMin.z);
115    break;
116  case 4:
117    a->SetValue(mMax.x, mMax.y, mMax.z);
118    b->SetValue(mMax.x, mMin.y, mMax.z);
119    break;
120  case 5:
121    a->SetValue(mMax.x, mMax.y, mMax.z);
122    b->SetValue(mMin.x, mMax.y, mMax.z);
123    break;
124   
125  case 6:
126    a->SetValue(mMin.x, mMin.y, mMax.z);
127    b->SetValue(mMin.x, mMax.y, mMax.z);
128    break;
129  case 7:
130    a->SetValue(mMin.x, mMin.y, mMax.z);
131    b->SetValue(mMax.x, mMin.y, mMax.z);
132    break;
133  case 8:
134    a->SetValue(mMin.x, mMax.y, mMin.z);
135    b->SetValue(mMin.x, mMax.y, mMax.z);
136    break;
137  case 9:
138    a->SetValue(mMin.x, mMax.y, mMin.z);
139    b->SetValue(mMax.x, mMax.y, mMin.z);
140    break;
141  case 10:
142    a->SetValue(mMax.x, mMin.y, mMin.z);
143    b->SetValue(mMax.x, mMax.y, mMin.z);
144    break;
145  case 11:
146    a->SetValue(mMax.x, mMin.y, mMin.z);
147    b->SetValue(mMax.x, mMin.y, mMax.z);
148    break;
149  }
150}
151
152// returns the vertex indices in the range <0..7>, v = 4.x + 2.y + z, where
153// x,y,z are either 0 or 1; (0 .. min coordinate, 1 .. max coordinate)
154void
155AxisAlignedBox3::GetEdge(const int edge, int  &aIdx, int &bIdx) const
156{
157  switch(edge) {
158  case 0: aIdx = 0; bIdx = 1; break;
159  case 1: aIdx = 0; bIdx = 2; break;
160  case 2: aIdx = 0; bIdx = 4; break;
161  case 3: aIdx = 7; bIdx = 6; break;
162  case 4: aIdx = 7; bIdx = 5; break;
163  case 5: aIdx = 7; bIdx = 3; break;
164  case 6: aIdx = 1; bIdx = 3; break;
165  case 7: aIdx = 1; bIdx = 5; break;
166  case 8: aIdx = 2; bIdx = 3; break;
167  case 9: aIdx = 2; bIdx = 6; break;
168  case 10: aIdx = 4; bIdx = 6; break;
169  case 11: aIdx = 4; bIdx = 5; break;
170  }
171}
172
173void
174AxisAlignedBox3::Include(const int &axis, const float &newBound)
175{
176  switch (axis) {
177    case 0: { // x-axis
178      if (mMin.x > newBound)
179        mMin.x = newBound;
180      if (mMax.x < newBound)
181        mMax.x = newBound;
182      break;
183    }
184    case 1: { // y-axis
185      if (mMin.y > newBound)
186        mMin.y = newBound;
187      if (mMax.y < newBound)
188        mMax.y = newBound;
189      break;
190    }
191    case 2: { // z-axis
192      if (mMin.z > newBound)
193        mMin.z = newBound;
194      if (mMax.z < newBound)
195        mMax.z = newBound;
196      break;
197    }
198  }
199}
200
201#if 0
202// ComputeMinMaxT computes the minimum and maximum signed distances
203// of intersection with the ray; it returns 1 if the ray hits
204// the bounding box and 0 if it does not.
205int
206AxisAlignedBox3::ComputeMinMaxT(const Ray &ray,
207                                float *tmin,
208                                float *tmax) const
209{
210  float minx, maxx, miny, maxy, minz, maxz;
211
212  if (fabs(ray.dir.x) < 0.001) {
213    if (mMin.x < ray.loc.x && mMax.x > ray.loc.x) {
214      minx = -MAXFLOAT;
215      maxx = MAXFLOAT;
216    }
217    else
218      return 0;
219  }
220  else {
221    float t1 = (mMin.x - ray.loc.x) / ray.dir.x;
222    float t2 = (mMax.x - ray.loc.x) / ray.dir.x;
223    if (t1 < t2) {
224      minx = t1;
225      maxx = t2;
226    }
227    else {
228      minx = t2;
229      maxx = t1;
230    }
231    if (maxx < 0)
232      return 0;
233  }
234
235  if (fabs(ray.dir.y) < 0.001) {
236    if (mMin.y < ray.loc.y && mMax.y > ray.loc.y) {
237      miny = -MAXFLOAT;
238      maxy = MAXFLOAT;
239    }
240    else
241      return 0;
242  }
243  else {
244    float t1 = (mMin.y - ray.loc.y) / ray.dir.y;
245    float t2 = (mMax.y - ray.loc.y) / ray.dir.y;
246    if (t1 < t2) {
247      miny = t1;
248      maxy = t2;
249    }
250    else {
251      miny = t2;
252      maxy = t1;
253    }
254    if (maxy < 0.0)
255      return 0;
256  }
257
258  if (fabs(ray.dir.z) < 0.001) {
259    if (mMin.z < ray.loc.z && mMax.z > ray.loc.z) {
260      minz = -MAXFLOAT;
261      maxz = MAXFLOAT;
262    }
263    else
264      return 0;
265  }
266  else {
267    float t1 = (mMin.z - ray.loc.z) / ray.dir.z;
268    float t2 = (mMax.z - ray.loc.z) / ray.dir.z;
269    if (t1 < t2) {
270      minz = t1;
271      maxz = t2;
272    }
273    else {
274      minz = t2;
275      maxz = t1;
276    }
277    if (maxz < 0.0)
278      return 0;
279  }
280
281  *tmin = minx;
282  if (miny > *tmin)
283    *tmin = miny;
284  if (minz > *tmin)
285    *tmin = minz;
286
287  *tmax = maxx;
288  if (maxy < *tmax)
289    *tmax = maxy;
290  if (maxz < *tmax)
291    *tmax = maxz;
292
293  return 1; // yes, intersection was found
294}
295#else
296// another variant of the same, with less variables
297int
298AxisAlignedBox3::ComputeMinMaxT(const Ray &ray,
299                                float *tmin,
300                                float *tmax) const
301{
302  const float dirEps = 1e-8f;
303  register float minx, maxx;
304  ray.ComputeInvertedDir();
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        //-- compute classification of vertices
1722        for (int i = 0; i < 8; ++i)
1723        {
1724                GetVertex(i, vtx);
1725                side[i] = plane.Side(vtx);
1726                if (side[i] > 0)
1727                        onFrontSide = true;
1728                else if (side[i] < 0)
1729                        onBackSide = true;
1730                else // vertex coincident => push_back
1731                        planePoly->mVertices.push_back(vtx);
1732        }
1733
1734        //-- find intersections
1735        if (onFrontSide && onBackSide)
1736        {
1737                Vector3 ptA, ptB;
1738                for (int i = 0; i < 12; ++ i)
1739                {
1740                        int aIdx, bIdx;
1741                        GetEdge(i, aIdx, bIdx);
1742
1743                        ptA = GetVertex(aIdx);
1744                        ptB = GetVertex(bIdx);
1745
1746                        int sideA = side[aIdx];
1747                        int sideB = side[bIdx];
1748
1749                        if (((sideA > 0) && (sideB < 0)) || (sideA < 0) && (sideB > 0))
1750                                planePoly->mVertices.push_back(plane.FindIntersection(ptA, ptB));       
1751                }
1752        }
1753
1754        // order intersections 
1755        if (planePoly->mVertices.size() > 3)
1756        {
1757                Vector3 centerOfMass(0);
1758                int i;
1759                // compute center of mass
1760                for (i = 0; i < (int)planePoly->mVertices.size(); ++ i)
1761                        centerOfMass += planePoly->mVertices[i];
1762               
1763                centerOfMass /= (float)planePoly->mVertices.size();
1764
1765                vector<VertexData> vertexData;
1766               
1767                Vector3 refVec = Normalize(centerOfMass - planePoly->mVertices[0]);
1768
1769                // compute angle to reference point
1770                for (i = 1; i < (int)planePoly->mVertices.size(); ++ i)
1771                  {
1772                    float angle =
1773                      Angle(refVec, centerOfMass - planePoly->mVertices[i], plane.mNormal);
1774                   
1775                    vertexData.push_back(VertexData(planePoly->mVertices[i], angle));
1776                  }
1777               
1778                std::stable_sort(vertexData.begin(), vertexData.end());
1779
1780                // update vertices
1781                for (i = 1; i < (int)planePoly->mVertices.size(); ++ i)
1782                        planePoly->mVertices[i] = vertexData[i - 1].mVertex;
1783        }
1784        else if (planePoly->mVertices.size() == 3)
1785        {
1786                // fix orientation if needed
1787                if (DotProd(planePoly->GetNormal(), plane.mNormal) < 0)
1788                {
1789                        Vector3 v = planePoly->mVertices[1];
1790                        planePoly->mVertices[1] = planePoly->mVertices[2];
1791                        planePoly->mVertices[2] = v;
1792                }
1793        }
1794       
1795        return planePoly;
1796}
1797
1798
1799bool AxisAlignedBox3::GetRaySegment(const Ray &ray,
1800                                                                        float &minT,
1801                                                                        float &maxT) const
1802{
1803        maxT = 1e15f;
1804        minT = 0;
1805       
1806        // test with bounding box
1807        if (!GetMinMaxT(ray, &minT, &maxT))
1808                return false;
1809
1810        if (minT < 0)
1811                minT = 0;
1812
1813        // bound ray or line segment
1814        if (//(ray.GetType() == Ray::LOCAL_RAY) &&
1815            !ray.intersections.empty() &&
1816                (ray.intersections[0].mT <= maxT))
1817        {
1818                maxT = ray.intersections[0].mT;
1819        }
1820
1821        return true;
1822}
1823
1824
1825  // Compute tmin and tmax for a ray, whenever required .. need not pierce box
1826int
1827AxisAlignedBox3::ComputeMinMaxT(const Vector3 &origin,
1828                                                                const Vector3 &direction,
1829                                                                float *tmin,
1830                                                                float *tmax) const
1831{
1832
1833  register float minx, maxx;
1834
1835 
1836  Vector3 invDirection;
1837  const float eps = 1e-6f;
1838  const float invEps = 1e6f;
1839 
1840  // it does change the ray direction very slightly,
1841  // but the size direction vector is not practically changed
1842 
1843  if (fabs(direction.x) < eps) {
1844    if (direction.x < 0.0f)
1845      invDirection.x = -invEps;
1846    else
1847      invDirection.x = invEps;
1848  }
1849  else
1850    invDirection.x = 1.0f / direction.x;
1851 
1852  if (fabs(direction.y) < eps) {
1853    if (direction.y < 0.0)
1854      invDirection.y = -invEps;
1855    else
1856      invDirection.y = invEps;
1857  }
1858  else
1859    invDirection.y = 1.0f / direction.y;
1860 
1861  if (fabs(direction.z) < eps) {
1862    if (direction.z < 0.0f)
1863      invDirection.z = -invEps;
1864    else
1865      invDirection.z = invEps;
1866  }
1867  else
1868    invDirection.z = 1.0f / direction.z;
1869
1870
1871 
1872  if (fabs(direction.x) < 0.001f) {
1873    if (mMin.x < origin.x && mMax.x > origin.x) {
1874      minx = -MAXFLOAT;
1875      maxx = MAXFLOAT;
1876    }
1877    else
1878      return 0;
1879  }
1880  else {
1881    float t1 = (mMin.x - origin.x) * invDirection.x;
1882    float t2 = (mMax.x - origin.x) * invDirection.x;
1883    if (t1 < t2) {
1884      minx = t1;
1885      maxx = t2;
1886    }
1887    else {
1888      minx = t2;
1889      maxx = t1;
1890    }
1891        //    if (maxx < 0.0)
1892        //      return 0;
1893  }
1894
1895  *tmin = minx;
1896  *tmax = maxx;
1897 
1898  if (fabs(direction.y) < 0.001) {
1899    if (mMin.y < origin.y && mMax.y > origin.y) {
1900      minx = -MAXFLOAT;
1901      maxx = MAXFLOAT;
1902    }
1903    else
1904      return 0;
1905  }
1906  else {
1907    float t1 = (mMin.y - origin.y) * invDirection.y;
1908    float t2 = (mMax.y - origin.y) * invDirection.y;
1909    if (t1 < t2) {
1910      minx = t1;
1911      maxx = t2;
1912    }
1913    else {
1914      minx = t2;
1915      maxx = t1;
1916    }
1917        //    if (maxx < 0.0)
1918        //      return 0;
1919  }
1920
1921  if (minx > *tmin)
1922    *tmin = minx;
1923  if (maxx < *tmax)
1924    *tmax = maxx;
1925 
1926  if (fabs(direction.z) < 0.001) {
1927    if (mMin.z < origin.z && mMax.z > origin.z) {
1928      minx = -MAXFLOAT;
1929      maxx = MAXFLOAT;
1930    }
1931    else
1932      return 0;
1933  }
1934  else {
1935    float t1 = (mMin.z - origin.z) * invDirection.z;
1936    float t2 = (mMax.z - origin.z) * invDirection.z;
1937    if (t1 < t2) {
1938      minx = t1;
1939      maxx = t2;
1940    }
1941    else {
1942      minx = t2;
1943      maxx = t1;
1944    }
1945        //    if (maxx < 0.0)
1946        //      return 0;
1947  }
1948
1949  if (minx > *tmin)
1950    *tmin = minx;
1951  if (maxx < *tmax)
1952    *tmax = maxx;
1953
1954  return 1; // yes, intersection was found
1955}
1956
1957
1958bool AxisAlignedBox3::GetIntersectionFace(Rectangle3 &face,
1959                                                                                  const AxisAlignedBox3 &neighbour) const
1960                                                                                   
1961{
1962        if (EpsilonEqual(mMin[0], neighbour.Max(0)))
1963        {
1964                float maxy = min(mMax.y, neighbour.mMax.y);
1965                float maxz = min(mMax.z, neighbour.mMax.z);
1966                float miny = max(mMin.y, neighbour.mMin.y);
1967                float minz = max(mMin.z, neighbour.mMin.z);
1968
1969                face.mVertices[3].SetValue(mMin.x, miny, minz);
1970                face.mVertices[2].SetValue(mMin.x, maxy, minz);
1971                face.mVertices[1].SetValue(mMin.x, maxy, maxz);
1972                face.mVertices[0].SetValue(mMin.x, miny, maxz);
1973
1974        return true;
1975        }
1976    if (EpsilonEqual(mMax[0], neighbour.Min(0)))
1977        {
1978                float maxy = min(mMax.y, neighbour.mMax.y);
1979                float maxz = min(mMax.z, neighbour.mMax.z);
1980                float miny = max(mMin.y, neighbour.mMin.y);
1981                float minz = max(mMin.z, neighbour.mMin.z);
1982
1983                face.mVertices[0].SetValue(mMax.x, miny, minz);
1984                face.mVertices[1].SetValue(mMax.x, maxy, minz);
1985                face.mVertices[2].SetValue(mMax.x, maxy, maxz);
1986                face.mVertices[3].SetValue(mMax.x, miny, maxz);
1987
1988                return true;
1989        }
1990    if (EpsilonEqual(mMin[1], neighbour.Max(1)))
1991        {
1992                float maxx = min(mMax.x, neighbour.mMax.x);
1993                float maxz = min(mMax.z, neighbour.mMax.z);
1994                float minx = max(mMin.x, neighbour.mMin.x);
1995                float minz = max(mMin.z, neighbour.mMin.z);
1996
1997                face.mVertices[3].SetValue(minx, mMin.y, minz);
1998                face.mVertices[2].SetValue(minx, mMin.y, maxz);
1999                face.mVertices[1].SetValue(maxx, mMin.y, maxz);
2000                face.mVertices[0].SetValue(maxx, mMin.y, minz);
2001               
2002                return true;
2003        }
2004        if (EpsilonEqual(mMax[1], neighbour.Min(1)))
2005        {               
2006                float maxx = min(mMax.x, neighbour.mMax.x);
2007                float maxz = min(mMax.z, neighbour.mMax.z);
2008                float minx = max(mMin.x, neighbour.mMin.x);
2009                float minz = max(mMin.z, neighbour.mMin.z);
2010
2011                face.mVertices[0].SetValue(minx, mMax.y, minz);
2012                face.mVertices[1].SetValue(minx, mMax.y, maxz);
2013                face.mVertices[2].SetValue(maxx, mMax.y, maxz);
2014                face.mVertices[3].SetValue(maxx, mMax.y, minz);
2015
2016                return true;
2017        }
2018        if (EpsilonEqual(mMin[2], neighbour.Max(2)))
2019        {
2020                float maxx = min(mMax.x, neighbour.mMax.x);
2021                float maxy = min(mMax.y, neighbour.mMax.y);
2022                float minx = max(mMin.x, neighbour.mMin.x);
2023                float miny = max(mMin.y, neighbour.mMin.y);
2024
2025                face.mVertices[3].SetValue(minx, miny, mMin.z);
2026                face.mVertices[2].SetValue(maxx, miny, mMin.z);
2027                face.mVertices[1].SetValue(maxx, maxy, mMin.z);
2028                face.mVertices[0].SetValue(minx, maxy, mMin.z);
2029       
2030                return true;
2031        }
2032        if (EpsilonEqual(mMax[2], neighbour.Min(2)))
2033        {
2034                float maxx = min(mMax.x, neighbour.mMax.x);
2035                float maxy = min(mMax.y, neighbour.mMax.y);
2036                float minx = max(mMin.x, neighbour.mMin.x);
2037                float miny = max(mMin.y, neighbour.mMin.y);
2038
2039                face.mVertices[0].SetValue(minx, miny, mMax.z);
2040                face.mVertices[1].SetValue(maxx, miny, mMax.z);
2041                face.mVertices[2].SetValue(maxx, maxy, mMax.z);
2042                face.mVertices[3].SetValue(minx, maxy, mMax.z);
2043
2044                return true;
2045        }
2046
2047        return false;
2048}
2049
2050
2051void IncludeBoxInMesh(const AxisAlignedBox3 &box, Mesh &mesh)
2052{
2053        // add 6 vertices of the box
2054        int index = (int)mesh.mVertices.size();
2055       
2056        for (int i=0; i < 8; ++ i)
2057        {
2058                Vector3 v;
2059                box.GetVertex(i, v);
2060                mesh.mVertices.push_back(v);
2061        }
2062       
2063        mesh.AddFace(new Face(index + 0, index + 1, index + 3, index + 2) );
2064        mesh.AddFace(new Face(index + 0, index + 2, index + 6, index + 4) );
2065        mesh.AddFace(new Face(index + 4, index + 6, index + 7, index + 5) );
2066       
2067        mesh.AddFace(new Face(index + 3, index + 1, index + 5, index + 7) );
2068        mesh.AddFace(new Face(index + 0, index + 4, index + 5, index + 1) );
2069        mesh.AddFace(new Face(index + 2, index + 3, index + 7, index + 6) );
2070}
2071
2072
2073void AxisAlignedBox3::ExtractPolys(PolygonContainer &polys) const
2074{
2075        Polygon3 *face1 = new Polygon3();
2076        polys.push_back(face1);   
2077
2078    face1->mVertices.push_back(Vector3(mMin.x, mMin.y, mMax.z));
2079        face1->mVertices.push_back(Vector3(mMin.x, mMax.y, mMax.z));
2080        face1->mVertices.push_back(Vector3(mMin.x, mMax.y ,mMin.z));
2081        face1->mVertices.push_back(Vector3(mMin.x, mMin.y, mMin.z));
2082
2083        Polygon3 *face2 = new Polygon3(); 
2084        polys.push_back(face2);
2085   
2086    face2->mVertices.push_back(Vector3(mMax.x, mMin.y, mMin.z));
2087    face2->mVertices.push_back(Vector3(mMax.x, mMax.y, mMin.z));
2088    face2->mVertices.push_back(Vector3(mMax.x, mMax.y, mMax.z));
2089    face2->mVertices.push_back(Vector3(mMax.x, mMin.y, mMax.z));
2090 
2091        Polygon3 *face3 = new Polygon3(); 
2092        polys.push_back(face3);
2093
2094    face3->mVertices.push_back(Vector3(mMax.x, mMin.y ,mMin.z));
2095        face3->mVertices.push_back(Vector3(mMax.x, mMin.y, mMax.z));
2096        face3->mVertices.push_back(Vector3(mMin.x, mMin.y, mMax.z));
2097        face3->mVertices.push_back(Vector3(mMin.x, mMin.y, mMin.z));
2098
2099        Polygon3 *face4 = new Polygon3(); 
2100        polys.push_back(face4);
2101
2102        face4->mVertices.push_back(Vector3(mMin.x, mMax.y, mMin.z));
2103        face4->mVertices.push_back(Vector3(mMin.x, mMax.y, mMax.z));
2104        face4->mVertices.push_back(Vector3(mMax.x, mMax.y, mMax.z));
2105        face4->mVertices.push_back(Vector3(mMax.x, mMax.y, mMin.z));
2106   
2107        Polygon3 *face5 = new Polygon3(); 
2108        polys.push_back(face5);
2109
2110        face5->mVertices.push_back(Vector3(mMin.x, mMax.y, mMin.z));
2111    face5->mVertices.push_back(Vector3(mMax.x, mMax.y, mMin.z));
2112    face5->mVertices.push_back(Vector3(mMax.x, mMin.y, mMin.z));
2113        face5->mVertices.push_back(Vector3(mMin.x, mMin.y, mMin.z));
2114
2115        Polygon3 *face6 = new Polygon3(); 
2116        polys.push_back(face6); 
2117 
2118    face6->mVertices.push_back(Vector3(mMin.x, mMin.y, mMax.z));
2119    face6->mVertices.push_back(Vector3(mMax.x, mMin.y, mMax.z));
2120    face6->mVertices.push_back(Vector3(mMax.x, mMax.y, mMax.z));
2121    face6->mVertices.push_back(Vector3(mMin.x, mMax.y, mMax.z));
2122
2123}
2124
2125
2126void AxisAlignedBox3::Split(const int axis,
2127                                                        const float value,
2128                                                        AxisAlignedBox3 &front,
2129                                                        AxisAlignedBox3 &back) const
2130{
2131        if ( (value >= mMin[axis]) && (value <= mMax[axis]) )
2132        {
2133                front.mMin = mMin; front.mMax = mMax;
2134                back.mMin = mMin; back.mMax = mMax;
2135
2136                back.mMax[axis] = value;
2137                front.mMin[axis] = value;
2138        }
2139}
2140
2141
2142}
Note: See TracBrowser for help on using the repository browser.