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

Revision 2176, 63.3 KB checked in by mattausch, 18 years ago (diff)

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