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

Revision 209, 44.7 KB checked in by bittner, 19 years ago (diff)

data added

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