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

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