FreeFOAM The Cross-Platform CFD Toolkit
treeBoundBox.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "treeBoundBox.H"
27 #include <OpenFOAM/ListOps.H>
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 const Foam::scalar Foam::treeBoundBox::great(GREAT);
32 
34 (
35  vector(-GREAT, -GREAT, -GREAT),
36  vector(GREAT, GREAT, GREAT)
37 );
38 
39 
41 (
42  vector(GREAT, GREAT, GREAT),
43  vector(-GREAT, -GREAT, -GREAT)
44 );
45 
46 
48 const Foam::label facesArray[6][4] =
49 {
50  {0, 4, 6, 2}, // left
51  {1, 3, 7, 5}, // right
52  {0, 1, 5, 4}, // bottom
53  {2, 6, 7, 3}, // top
54  {0, 2, 3, 1}, // back
55  {4, 5, 7, 6} // front
56 };
58 
59 
61 (
62  initListList<face, label, 6, 4>(facesArray)
63 );
64 
65 
67 const Foam::label edgesArray[12][2] =
68 {
69  {0, 1}, // 0
70  {1, 3},
71  {2, 3}, // 2
72  {0, 2},
73  {4, 5}, // 4
74  {5, 7},
75  {6, 7}, // 6
76  {4, 6},
77  {0, 4}, // 8
78  {1, 5},
79  {3, 7}, // 10
80  {2, 6}
81 };
83 
84 
86 (
87  initListList<edge, label, 12, 2>(edgesArray)
88 );
89 
90 
92 (
93  calcFaceNormals()
94 );
95 
96 
97 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
98 
99 Foam::FixedList<Foam::vector, 6> Foam::treeBoundBox::calcFaceNormals()
100 {
101  FixedList<vector, 6> normals;
102  normals[LEFT] = vector(-1, 0, 0);
103  normals[RIGHT] = vector( 1, 0, 0);
104  normals[BOTTOM] = vector( 0, -1, 0);
105  normals[TOP] = vector( 0, 1, 0);
106  normals[BACK] = vector( 0, 0, -1);
107  normals[FRONT] = vector( 0, 0, 1);
108  return normals;
109 }
110 
111 
112 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
113 
114 // Construct as the bounding box of the given pointField
116 :
117  boundBox()
118 {
119  if (points.empty())
120  {
121  WarningIn
122  (
123  "treeBoundBox::treeBoundBox(const UList<point>&)"
124  ) << "cannot find bounding box for zero-sized pointField"
125  << "returning zero" << endl;
126 
127  return;
128  }
129 
130  min() = points[0];
131  max() = points[0];
132 
133  for (label i = 1; i < points.size(); i++)
134  {
135  min() = ::Foam::min(min(), points[i]);
136  max() = ::Foam::max(max(), points[i]);
137  }
138 }
139 
140 
141 // Construct as the bounding box of the given pointField
143 (
144  const UList<point>& points,
145  const UList<label>& meshPoints
146 )
147 :
148  boundBox()
149 {
150  if (points.empty() || meshPoints.empty())
151  {
152  WarningIn
153  (
154  "treeBoundBox::treeBoundBox"
155  "(const UList<point>&, const UList<label>&)"
156  ) << "cannot find bounding box for zero-sized pointField"
157  << "returning zero" << endl;
158 
159  return;
160  }
161 
162  min() = points[meshPoints[0]];
163  max() = points[meshPoints[0]];
164 
165  for (label i = 1; i < meshPoints.size(); i++)
166  {
167  min() = ::Foam::min(min(), points[meshPoints[i]]);
168  max() = ::Foam::max(max(), points[meshPoints[i]]);
169  }
170 }
171 
172 
173 // Construct from Istream
175 :
176  boundBox(is)
177 {}
178 
179 
180 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
181 
183 {
184  pointField points(8);
185  forAll(points, octant)
186  {
187  points[octant] = corner(octant);
188 
189  }
190  return points;
191 }
192 
193 
195 {
196  return subBbox(midpoint(), octant);
197 }
198 
199 
200 // Octant to bounding box using permutation only.
202 (
203  const point& mid,
204  const direction octant
205 ) const
206 {
207  if (octant > 7)
208  {
210  (
211  "treeBoundBox::subBbox(const point&, const direction)"
212  ) << "octant should be [0..7]"
213  << abort(FatalError);
214  }
215 
216  // start with a copy of this bounding box and adjust limits accordingly
217  treeBoundBox subBb(*this);
218  point& bbMin = subBb.min();
219  point& bbMax = subBb.max();
220 
221  if (octant & treeBoundBox::RIGHTHALF)
222  {
223  bbMin.x() = mid.x(); // mid -> max
224  }
225  else
226  {
227  bbMax.x() = mid.x(); // min -> mid
228  }
229 
230  if (octant & treeBoundBox::TOPHALF)
231  {
232  bbMin.y() = mid.y(); // mid -> max
233  }
234  else
235  {
236  bbMax.y() = mid.y(); // min -> mid
237  }
238 
239  if (octant & treeBoundBox::FRONTHALF)
240  {
241  bbMin.z() = mid.z(); // mid -> max
242  }
243  else
244  {
245  bbMax.z() = mid.z(); // min -> mid
246  }
247 
248  return subBb;
249 }
250 
251 
253 (
254  const point& centre,
255  const scalar radiusSqr
256 ) const
257 {
258  // Find out where centre is in relation to bb.
259  // Find nearest point on bb.
260  scalar distSqr = 0;
261 
262  for (direction dir = 0; dir < vector::nComponents; dir++)
263  {
264  scalar d0 = min()[dir] - centre[dir];
265  scalar d1 = max()[dir] - centre[dir];
266 
267  if ((d0 > 0) != (d1 > 0))
268  {
269  // centre inside both extrema. This component does not add any
270  // distance.
271  }
272  else if (Foam::mag(d0) < Foam::mag(d1))
273  {
274  distSqr += d0*d0;
275  }
276  else
277  {
278  distSqr += d1*d1;
279  }
280 
281  if (distSqr > radiusSqr)
282  {
283  return false;
284  }
285  }
286 
287  return true;
288 }
289 
290 
291 // line intersection. Returns true if line (start to end) inside
292 // bb or intersects bb. Sets pt to intersection.
293 //
294 // Sutherlands algorithm:
295 // loop
296 // - start = intersection of line with one of the planes bounding
297 // the bounding box
298 // - stop if start inside bb (return true)
299 // - stop if start and end in same 'half' (e.g. both above bb)
300 // (return false)
301 //
302 // Uses posBits to efficiently determine 'half' in which start and end
303 // point are.
304 //
305 // Note:
306 // - sets coordinate to exact position: e.g. pt.x() = min().x()
307 // since plane intersect routine might have truncation error.
308 // This makes sure that posBits tests 'inside'
310 (
311  const point& overallStart,
312  const vector& overallVec,
313  const point& start,
314  const point& end,
315  point& pt,
316  direction& ptOnFaces
317 ) const
318 {
319  const direction endBits = posBits(end);
320  pt = start;
321 
322  while (true)
323  {
324  direction ptBits = posBits(pt);
325 
326  if (ptBits == 0)
327  {
328  // pt inside bb
329  ptOnFaces = faceBits(pt);
330  return true;
331  }
332 
333  if ((ptBits & endBits) != 0)
334  {
335  // pt and end in same block outside of bb
336  ptOnFaces = faceBits(pt);
337  return false;
338  }
339 
340  if (ptBits & LEFTBIT)
341  {
342  // Intersect with plane V=min, n=-1,0,0
343  if (Foam::mag(overallVec.x()) > VSMALL)
344  {
345  scalar s = (min().x() - overallStart.x())/overallVec.x();
346  pt.x() = min().x();
347  pt.y() = overallStart.y() + overallVec.y()*s;
348  pt.z() = overallStart.z() + overallVec.z()*s;
349  }
350  else
351  {
352  // Vector not in x-direction. But still intersecting bb planes.
353  // So must be close - just snap to bb.
354  pt.x() = min().x();
355  }
356  }
357  else if (ptBits & RIGHTBIT)
358  {
359  // Intersect with plane V=max, n=1,0,0
360  if (Foam::mag(overallVec.x()) > VSMALL)
361  {
362  scalar s = (max().x() - overallStart.x())/overallVec.x();
363  pt.x() = max().x();
364  pt.y() = overallStart.y() + overallVec.y()*s;
365  pt.z() = overallStart.z() + overallVec.z()*s;
366  }
367  else
368  {
369  pt.x() = max().x();
370  }
371  }
372  else if (ptBits & BOTTOMBIT)
373  {
374  // Intersect with plane V=min, n=0,-1,0
375  if (Foam::mag(overallVec.y()) > VSMALL)
376  {
377  scalar s = (min().y() - overallStart.y())/overallVec.y();
378  pt.x() = overallStart.x() + overallVec.x()*s;
379  pt.y() = min().y();
380  pt.z() = overallStart.z() + overallVec.z()*s;
381  }
382  else
383  {
384  pt.x() = min().y();
385  }
386  }
387  else if (ptBits & TOPBIT)
388  {
389  // Intersect with plane V=max, n=0,1,0
390  if (Foam::mag(overallVec.y()) > VSMALL)
391  {
392  scalar s = (max().y() - overallStart.y())/overallVec.y();
393  pt.x() = overallStart.x() + overallVec.x()*s;
394  pt.y() = max().y();
395  pt.z() = overallStart.z() + overallVec.z()*s;
396  }
397  else
398  {
399  pt.y() = max().y();
400  }
401  }
402  else if (ptBits & BACKBIT)
403  {
404  // Intersect with plane V=min, n=0,0,-1
405  if (Foam::mag(overallVec.z()) > VSMALL)
406  {
407  scalar s = (min().z() - overallStart.z())/overallVec.z();
408  pt.x() = overallStart.x() + overallVec.x()*s;
409  pt.y() = overallStart.y() + overallVec.y()*s;
410  pt.z() = min().z();
411  }
412  else
413  {
414  pt.z() = min().z();
415  }
416  }
417  else if (ptBits & FRONTBIT)
418  {
419  // Intersect with plane V=max, n=0,0,1
420  if (Foam::mag(overallVec.z()) > VSMALL)
421  {
422  scalar s = (max().z() - overallStart.z())/overallVec.z();
423  pt.x() = overallStart.x() + overallVec.x()*s;
424  pt.y() = overallStart.y() + overallVec.y()*s;
425  pt.z() = max().z();
426  }
427  else
428  {
429  pt.z() = max().z();
430  }
431  }
432  }
433 }
434 
435 
437 (
438  const point& start,
439  const point& end,
440  point& pt
441 ) const
442 {
443  direction ptBits;
444  return intersects(start, end-start, start, end, pt, ptBits);
445 }
446 
447 
448 // this.bb fully contains bb
450 {
451  return contains(bb.min()) && contains(bb.max());
452 }
453 
454 
455 bool Foam::treeBoundBox::contains(const vector& dir, const point& pt) const
456 {
457  //
458  // Compare all components against min and max of bb
459  //
460  for (direction cmpt=0; cmpt<3; cmpt++)
461  {
462  if (pt[cmpt] < min()[cmpt])
463  {
464  return false;
465  }
466  else if (pt[cmpt] == min()[cmpt])
467  {
468  // On edge. Outside if direction points outwards.
469  if (dir[cmpt] < 0)
470  {
471  return false;
472  }
473  }
474 
475  if (pt[cmpt] > max()[cmpt])
476  {
477  return false;
478  }
479  else if (pt[cmpt] == max()[cmpt])
480  {
481  // On edge. Outside if direction points outwards.
482  if (dir[cmpt] > 0)
483  {
484  return false;
485  }
486  }
487  }
488 
489  // All components inside bb
490  return true;
491 }
492 
493 
494 // Code position of pt on bounding box faces
496 {
497  direction faceBits = 0;
498  if (pt.x() == min().x())
499  {
500  faceBits |= LEFTBIT;
501  }
502  else if (pt.x() == max().x())
503  {
504  faceBits |= RIGHTBIT;
505  }
506 
507  if (pt.y() == min().y())
508  {
509  faceBits |= BOTTOMBIT;
510  }
511  else if (pt.y() == max().y())
512  {
513  faceBits |= TOPBIT;
514  }
515 
516  if (pt.z() == min().z())
517  {
518  faceBits |= BACKBIT;
519  }
520  else if (pt.z() == max().z())
521  {
522  faceBits |= FRONTBIT;
523  }
524  return faceBits;
525 }
526 
527 
528 // Code position of point relative to box
530 {
531  direction posBits = 0;
532 
533  if (pt.x() < min().x())
534  {
535  posBits |= LEFTBIT;
536  }
537  else if (pt.x() > max().x())
538  {
539  posBits |= RIGHTBIT;
540  }
541 
542  if (pt.y() < min().y())
543  {
544  posBits |= BOTTOMBIT;
545  }
546  else if (pt.y() > max().y())
547  {
548  posBits |= TOPBIT;
549  }
550 
551  if (pt.z() < min().z())
552  {
553  posBits |= BACKBIT;
554  }
555  else if (pt.z() > max().z())
556  {
557  posBits |= FRONTBIT;
558  }
559  return posBits;
560 }
561 
562 
563 // nearest and furthest corner coordinate.
564 // !names of treeBoundBox::min() and treeBoundBox::max() are confusing!
566 (
567  const point& pt,
568  point& nearest,
569  point& furthest
570 ) const
571 {
572  scalar nearX, nearY, nearZ;
573  scalar farX, farY, farZ;
574 
575  if (Foam::mag(min().x() - pt.x()) < Foam::mag(max().x() - pt.x()))
576  {
577  nearX = min().x();
578  farX = max().x();
579  }
580  else
581  {
582  nearX = max().x();
583  farX = min().x();
584  }
585 
586  if (Foam::mag(min().y() - pt.y()) < Foam::mag(max().y() - pt.y()))
587  {
588  nearY = min().y();
589  farY = max().y();
590  }
591  else
592  {
593  nearY = max().y();
594  farY = min().y();
595  }
596 
597  if (Foam::mag(min().z() - pt.z()) < Foam::mag(max().z() - pt.z()))
598  {
599  nearZ = min().z();
600  farZ = max().z();
601  }
602  else
603  {
604  nearZ = max().z();
605  farZ = min().z();
606  }
607 
608  nearest = point(nearX, nearY, nearZ);
609  furthest = point(farX, farY, farZ);
610 }
611 
612 
613 Foam::scalar Foam::treeBoundBox::maxDist(const point& pt) const
614 {
615  point near, far;
616  calcExtremities(pt, near, far);
617 
618  return Foam::mag(far - pt);
619 }
620 
621 
622 // Distance comparator
623 // Compare all vertices of bounding box against all of other bounding
624 // box to see if all vertices of one are nearer
626 (
627  const point& pt,
628  const treeBoundBox& other
629 ) const
630 {
631  //
632  // Distance point <-> nearest and furthest away vertex of this
633  //
634 
635  point nearThis, farThis;
636 
637  // get nearest and furthest away vertex
638  calcExtremities(pt, nearThis, farThis);
639 
640  const scalar minDistThis =
641  sqr(nearThis.x() - pt.x())
642  + sqr(nearThis.y() - pt.y())
643  + sqr(nearThis.z() - pt.z());
644  const scalar maxDistThis =
645  sqr(farThis.x() - pt.x())
646  + sqr(farThis.y() - pt.y())
647  + sqr(farThis.z() - pt.z());
648 
649  //
650  // Distance point <-> other
651  //
652 
653  point nearOther, farOther;
654 
655  // get nearest and furthest away vertex
656  other.calcExtremities(pt, nearOther, farOther);
657 
658  const scalar minDistOther =
659  sqr(nearOther.x() - pt.x())
660  + sqr(nearOther.y() - pt.y())
661  + sqr(nearOther.z() - pt.z());
662  const scalar maxDistOther =
663  sqr(farOther.x() - pt.x())
664  + sqr(farOther.y() - pt.y())
665  + sqr(farOther.z() - pt.z());
666 
667  //
668  // Categorize
669  //
670  if (maxDistThis < minDistOther)
671  {
672  // All vertices of this are nearer to point than any vertex of other
673  return -1;
674  }
675  else if (minDistThis > maxDistOther)
676  {
677  // All vertices of this are further from point than any vertex of other
678  return 1;
679  }
680  else
681  {
682  // Mixed bag
683  return 0;
684  }
685 }
686 
687 
688 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
689 
690 bool Foam::operator==(const treeBoundBox& a, const treeBoundBox& b)
691 {
692  return operator==
693  (
694  static_cast<const boundBox&>(a),
695  static_cast<const boundBox&>(b)
696  );
697 }
698 
699 
700 bool Foam::operator!=(const treeBoundBox& a, const treeBoundBox& b)
701 {
702  return !(a == b);
703 }
704 
705 
706 // * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
707 
708 Foam::Ostream& Foam::operator<<(Ostream& os, const treeBoundBox& bb)
709 {
710  return os << static_cast<const boundBox&>(bb);
711 }
712 
713 
714 Foam::Istream& Foam::operator>>(Istream& is, treeBoundBox& bb)
715 {
716  return is >> static_cast<boundBox&>(bb);
717 }
718 
719 
720 // ************************ vim: set sw=4 sts=4 et: ************************ //