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

Revision 2691, 64.7 KB checked in by mattausch, 16 years ago (diff)

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