0030997: Foundation Classes - name correction of dump macros
[occt.git] / src / Bnd / Bnd_OBB.cxx
1 // Created by: Eugeny MALTCHIKOV
2 // Copyright (c) 2017 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
15 #include <Bnd_OBB.hxx>
16
17 #include <Bnd_Tools.hxx>
18 #include <Bnd_Range.hxx>
19
20 #include <BVH_BoxSet.hxx>
21 #include <BVH_LinearBuilder.hxx>
22
23 #include <BVH_Traverse.hxx>
24
25 #include <NCollection_Array1.hxx>
26 #include <Precision.hxx>
27 #include <Standard_Dump.hxx>
28 #include <TColStd_Array1OfReal.hxx>
29
30 //! Auxiliary class to select from the points stored in
31 //! BVH tree the two points giving the extreme projection
32 //! parameters on the axis
33 class OBB_ExtremePointsSelector : 
34   public BVH_Traverse <Standard_Real, 3, BVH_BoxSet <Standard_Real, 3, gp_XYZ>, Bnd_Range>
35 {
36 public:
37
38   //! Constructor
39   OBB_ExtremePointsSelector() :
40     BVH_Traverse <Standard_Real, 3, BVH_BoxSet <Standard_Real, 3, gp_XYZ>, Bnd_Range>(),
41     myPrmMin (RealLast()),
42     myPrmMax (RealFirst())
43   {}
44
45 public: //! @name Set axis for projection
46
47   //! Sets the axis
48   void SetAxis (const gp_XYZ& theAxis) { myAxis = theAxis; }
49
50 public: //! @name Clears the points from previous runs
51
52   //! Clear
53   void Clear()
54   {
55     myPrmMin = RealLast();
56     myPrmMax = RealFirst();
57   }
58
59 public: //! @name Getting the results
60
61   //! Returns the minimal projection parameter
62   Standard_Real MinPrm() const { return myPrmMin; }
63
64   //! Returns the maximal projection parameter
65   Standard_Real MaxPrm() const { return myPrmMax; }
66
67   //! Returns the minimal projection point
68   const gp_XYZ& MinPnt() const { return myPntMin; }
69
70   //! Returns the maximal projection point
71   const gp_XYZ& MaxPnt() const { return myPntMax; }
72
73 public: //! @name Definition of rejection/acceptance rules
74
75   //! Defines the rules for node rejection
76   virtual Standard_Boolean RejectNode (const BVH_Vec3d& theCMin,
77                                        const BVH_Vec3d& theCMax,
78                                        Bnd_Range& theMetric) const Standard_OVERRIDE
79   {
80     if (myPrmMin > myPrmMax)
81       // No parameters computed yet
82       return Standard_False;
83
84     Standard_Real aPrmMin = myPrmMin, aPrmMax = myPrmMax;
85     Standard_Boolean isToReject = Standard_True;
86
87     // Check if the current node is between already found parameters
88     for (Standard_Integer i = 0; i < 2; ++i)
89     {
90       Standard_Real x = !i ? theCMin.x() : theCMax.x();
91       for (Standard_Integer j = 0; j < 2; ++j)
92       {
93         Standard_Real y = !j ? theCMin.y() : theCMax.y();
94         for (Standard_Integer k = 0; k < 2; ++k)
95         {
96           Standard_Real z = !k ? theCMin.z() : theCMax.z();
97
98           Standard_Real aPrm = myAxis.Dot (gp_XYZ (x, y, z));
99           if (aPrm < aPrmMin)
100           {
101             aPrmMin = aPrm;
102             isToReject = Standard_False;
103           }
104           else if (aPrm > aPrmMax)
105           {
106             aPrmMax = aPrm;
107             isToReject = Standard_False;
108           }
109         }
110       }
111     }
112
113     theMetric = Bnd_Range (aPrmMin, aPrmMax);
114
115     return isToReject;
116   }
117
118   //! Rules for node rejection by the metric
119   virtual Standard_Boolean RejectMetric (const Bnd_Range& theMetric) const Standard_OVERRIDE
120   {
121     if (myPrmMin > myPrmMax)
122       // no parameters computed
123       return Standard_False;
124
125     Standard_Real aMin, aMax;
126     if (!theMetric.GetBounds (aMin, aMax))
127       // void metric
128       return Standard_False;
129
130     // Check if the box of the branch is inside of the already computed parameters
131     return aMin > myPrmMin && aMax < myPrmMax;
132   }
133
134   //! Defines the rules for leaf acceptance
135   virtual Standard_Boolean Accept (const Standard_Integer theIndex,
136                                    const Bnd_Range&) Standard_OVERRIDE
137   {
138     const gp_XYZ& theLeaf = myBVHSet->Element (theIndex);
139     Standard_Real aPrm = myAxis.Dot (theLeaf);
140     if (aPrm < myPrmMin)
141     {
142       myPrmMin = aPrm;
143       myPntMin = theLeaf;
144     }
145     if (aPrm > myPrmMax)
146     {
147       myPrmMax = aPrm;
148       myPntMax = theLeaf;
149     }
150     return Standard_True;
151   }
152
153 public: //! @name Choosing the best branch
154
155   //! Returns true if the metric of the left branch is better than the metric of the right
156   virtual Standard_Boolean IsMetricBetter (const Bnd_Range& theLeft,
157                                            const Bnd_Range& theRight) const Standard_OVERRIDE
158   {
159     if (myPrmMin > myPrmMax)
160       // no parameters computed
161       return Standard_True;
162
163     Standard_Real aMin[2], aMax[2];
164     if (!theLeft.GetBounds  (aMin[0], aMax[0]) ||
165         !theRight.GetBounds (aMin[1], aMax[1]))
166       // void metrics
167       return Standard_True;
168
169     // Choose branch with larger extension over computed parameters
170     Standard_Real anExt[2] = {0.0, 0.0};
171     for (int i = 0; i < 2; ++i)
172     {
173       if (aMin[i] < myPrmMin) anExt[i] += myPrmMin - aMin[i];
174       if (aMax[i] > myPrmMax) anExt[i] += aMax[i] - myPrmMax;
175     }
176     return anExt[0] > anExt[1];
177   }
178
179 protected: //! @name Fields
180
181   gp_XYZ myAxis;          //!< Axis to project the points to
182   Standard_Real myPrmMin; //!< Minimal projection parameter
183   Standard_Real myPrmMax; //!< Maximal projection parameter
184   gp_XYZ myPntMin;        //!< Minimal projection point
185   gp_XYZ myPntMax;        //!< Maximal projection point
186 };
187
188 //! Tool for OBB construction
189 class OBBTool
190 {
191 public:
192   //! Constructor. theL - list of points.
193   //! theLT is a pointer to the list of tolerances
194   //! (i-th element of this array is a tolerance
195   //! of i-th point in theL). If theLT is empty
196   //! then the tolerance of every point is equal to 0.
197   //! Attention! The objects, which theL and theLT links on,
198   //! must be available during all time of OBB creation
199   //! (i.e. while the object of OBBTool exists).
200   OBBTool(const TColgp_Array1OfPnt& theL,
201           const TColStd_Array1OfReal *theLT = 0,
202           Standard_Boolean theIsOptimal = Standard_False);
203
204   //! DiTO algorithm for OBB construction
205   //! (http://www.idt.mdh.se/~tla/publ/FastOBBs.pdf)
206   void ProcessDiTetrahedron();
207
208   //! Creates OBB with already computed parameters
209   void BuildBox(Bnd_OBB& theBox);
210
211 protected:
212
213   // Computes the extreme points on the set of Initial axes
214   void ComputeExtremePoints ();
215
216   //! Works with the triangle set by the points in myTriIdx.
217   //! If theIsBuiltTrg == TRUE, new set of triangles will be
218   //! recomputed.
219   void ProcessTriangle(const Standard_Integer theIdx1,
220                        const Standard_Integer theIdx2,
221                        const Standard_Integer theIdx3,
222                        const Standard_Boolean theIsBuiltTrg);
223
224   //! Computes myTriIdx[2]
225   void FillToTriangle3();
226
227   //! Computes myTriIdx[3] and myTriIdx[4]
228   void FillToTriangle5(const gp_XYZ& theNormal,
229                        const gp_XYZ& theBarryCenter);
230
231   //! Returns half of the Surface area of the box
232   static Standard_Real ComputeQuality(const Standard_Real* const thePrmArr)
233   {
234     const Standard_Real aDX = thePrmArr[1] - thePrmArr[0],
235                         aDY = thePrmArr[3] - thePrmArr[2],
236                         aDZ = thePrmArr[5] - thePrmArr[4];
237
238     return (aDX*aDY + aDY*aDZ + aDX*aDZ);
239   }
240
241 protected:
242   //! Assignment operator is forbidden
243   OBBTool& operator=(const OBBTool&);
244
245 private:
246   //! Params structure stores the two values meaning
247   //! min and max parameters on the axis
248   struct Params
249   {
250     Params() :
251       _ParamMin(RealLast()), _ParamMax(RealFirst())
252     {}
253
254     Params(Standard_Real theMin, Standard_Real theMax)
255       : _ParamMin(theMin), _ParamMax(theMax)
256     {}
257
258     Standard_Real _ParamMin;
259     Standard_Real _ParamMax;
260   };
261
262   //! Computes the Minimal and maximal parameters on the vector
263   //! connecting the points myLExtremalPoints[theId1] and myLExtremalPoints[theId2]
264   void ComputeParams (const Standard_Integer theId1,
265                       const Standard_Integer theId2,
266                       Standard_Real &theMin,
267                       Standard_Real &theMax)
268   {
269     theMin = myParams[theId1][theId2]._ParamMin;
270     theMax = myParams[theId1][theId2]._ParamMax;
271
272     if (theMin > theMax)
273     {
274       FindMinMax ((myLExtremalPoints[theId1] - myLExtremalPoints[theId2]).Normalized(), theMin, theMax);
275       myParams[theId1][theId2]._ParamMin = myParams[theId2][theId1]._ParamMin = theMin;
276       myParams[theId1][theId2]._ParamMax = myParams[theId2][theId1]._ParamMax = theMax;
277     }
278   }
279
280   //! Looks for the min-max parameters on the axis.
281   //! For optimal case projects all the points on the axis,
282   //! for not optimal - only the set of extreme points.
283   void FindMinMax (const gp_XYZ& theAxis,
284                    Standard_Real &theMin,
285                    Standard_Real &theMax)
286   {
287     theMin = RealLast(), theMax = RealFirst();
288
289     if (myOptimal)
290       Project (theAxis, theMin, theMax);
291     else
292     {
293       for (Standard_Integer i = 0; i < myNbExtremalPoints; ++i)
294       {
295         Standard_Real aPrm = theAxis.Dot (myLExtremalPoints[i]);
296         if (aPrm < theMin) theMin = aPrm;
297         if (aPrm > theMax) theMax = aPrm;
298       }
299     }
300   }
301
302   //! Projects the set of points on the axis
303   void Project (const gp_XYZ& theAxis,
304                 Standard_Real& theMin, Standard_Real& theMax,
305                 gp_XYZ* thePntMin = 0, gp_XYZ* thePntMax = 0)
306   {
307     theMin = RealLast(), theMax = RealFirst();
308
309     if (myOptimal)
310     {
311       // Project BVH
312       OBB_ExtremePointsSelector anExtremePointsSelector;
313       anExtremePointsSelector.SetBVHSet (myPointBoxSet.get());
314       anExtremePointsSelector.SetAxis (theAxis);
315       anExtremePointsSelector.Select();
316       theMin = anExtremePointsSelector.MinPrm();
317       theMax = anExtremePointsSelector.MaxPrm();
318       if (thePntMin) *thePntMin = anExtremePointsSelector.MinPnt();
319       if (thePntMax) *thePntMax = anExtremePointsSelector.MaxPnt();
320     }
321     else
322     {
323       // Project all points
324       for (Standard_Integer iP = myPntsList.Lower(); iP <= myPntsList.Upper(); ++iP)
325       {
326         const gp_XYZ& aPoint = myPntsList(iP).XYZ();
327         const Standard_Real aPrm = theAxis.Dot (aPoint);
328         if (aPrm < theMin)
329         {
330           theMin = aPrm;
331           if (thePntMin)
332             *thePntMin = aPoint;
333         }
334         if (aPrm > theMax)
335         {
336           theMax = aPrm;
337           if (thePntMax)
338             *thePntMax = aPoint;
339         }
340       }
341     }
342   }
343
344 private:
345
346   //! Number of the initial axes.
347   static const Standard_Integer myNbInitAxes = 7;
348
349   //! Number of extremal points
350   static const Standard_Integer myNbExtremalPoints = 2 * myNbInitAxes;
351
352   //! The source list of points
353   const TColgp_Array1OfPnt& myPntsList;
354
355   //! Pointer to the array of tolerances
356   const TColStd_Array1OfReal *myListOfTolers;
357
358   //! Points of ditetrahedron
359   //! given by their indices in myLExtremalPoints.
360   Standard_Integer myTriIdx[5];
361
362   //! List of extremal points
363   gp_XYZ myLExtremalPoints[myNbExtremalPoints];
364
365   //! The axes of the box (always normalized or
366   //! can be null-vector)
367   gp_XYZ myAxes[3];
368
369   //! The surface area of the OBB
370   Standard_Real myQualityCriterion;
371
372   //! Defines if the OBB should be computed more tight.
373   //! Takes more time, but the volume is less.
374   Standard_Boolean myOptimal;
375
376   //! Point box set organized with BVH
377   opencascade::handle<BVH_BoxSet <Standard_Real, 3, gp_XYZ>> myPointBoxSet;
378
379   //! Stored min/max parameters for the axes between extremal points
380   Params myParams[myNbExtremalPoints][myNbExtremalPoints];
381 };
382
383 //=======================================================================
384 // Function : SetMinMax
385 // purpose : 
386 //    ATTENTION!!! thePrmArr must be initialized before this method calling.
387 //=======================================================================
388 static inline void SetMinMax(Standard_Real* const thePrmArr,
389                              const Standard_Real theNewParam)
390 {
391   if(theNewParam < thePrmArr[0])
392   {
393     thePrmArr[0] = theNewParam;
394   }
395   if(theNewParam > thePrmArr[1])
396   {
397     thePrmArr[1] = theNewParam;
398   }
399 }
400
401 //=======================================================================
402 // Function : Constructor
403 // purpose : 
404 //=======================================================================
405 OBBTool::
406     OBBTool(const TColgp_Array1OfPnt& theL,
407             const TColStd_Array1OfReal *theLT,
408             const Standard_Boolean theIsOptimal) : myPntsList(theL),
409                                                    myListOfTolers(theLT),
410                                                    myQualityCriterion(RealLast()),
411                                                    myOptimal (theIsOptimal)
412 {
413   if (myOptimal)
414   {
415     // Use linear builder for BVH construction with 30 elements in the leaf
416     opencascade::handle<BVH_LinearBuilder<Standard_Real, 3> > aLBuilder =
417       new BVH_LinearBuilder<Standard_Real, 3> (30);
418     myPointBoxSet = new BVH_BoxSet <Standard_Real, 3, gp_XYZ> (aLBuilder);
419     myPointBoxSet->SetSize(myPntsList.Length());
420
421     // Add the points into Set
422     for (Standard_Integer iP = 0; iP < theL.Length(); ++iP)
423     {
424       const gp_Pnt& aP = theL (iP);
425       Standard_Real aTol = theLT ? theLT->Value(iP) : Precision::Confusion();
426       BVH_Box <Standard_Real, 3> aBox (BVH_Vec3d (aP.X() - aTol, aP.Y() - aTol, aP.Z() - aTol),
427                                        BVH_Vec3d (aP.X() + aTol, aP.Y() + aTol, aP.Z() + aTol));
428       myPointBoxSet->Add (aP.XYZ(), aBox);
429     }
430
431     myPointBoxSet->Build();
432   }
433
434   ComputeExtremePoints();
435 }
436
437 //=======================================================================
438 // Function : ComputeExtremePoints
439 // purpose : 
440 //=======================================================================
441 void OBBTool::ComputeExtremePoints()
442 {
443   // Six initial axes show great quality on the Optimal OBB, plus
444   // the performance is better (due to the less number of operations).
445   // But they show worse quality for the not optimal approach.
446   //const Standard_Real a = (sqrt(5) - 1) / 2.;
447   //const gp_XYZ anInitialAxes6[myNbInitAxes] = { gp_XYZ (0, 1, a),
448   //                                              gp_XYZ (0, 1, -a),
449   //                                              gp_XYZ (1, a, 0),
450   //                                              gp_XYZ (1, -a, 0),
451   //                                              gp_XYZ (a, 0, 1),
452   //                                              gp_XYZ (a, 0, -1) };
453   const Standard_Real aSqrt3 = Sqrt(3);
454   const gp_XYZ anInitialAxes7[myNbInitAxes] = { gp_XYZ (1.0, 0.0, 0.0),
455                                                 gp_XYZ (0.0, 1.0, 0.0),
456                                                 gp_XYZ (0.0, 0.0, 1.0),
457                                                 gp_XYZ (1.0, 1.0, 1.0) / aSqrt3,
458                                                 gp_XYZ (1.0, 1.0, -1.0) / aSqrt3,
459                                                 gp_XYZ (1.0, -1.0, 1.0) / aSqrt3,
460                                                 gp_XYZ (1.0, -1.0, -1.0) / aSqrt3 };
461
462   // Set of initial axes
463   const gp_XYZ *anInitialAxesArray = anInitialAxes7;
464
465   // Min and Max parameter
466   Standard_Real aParams[myNbExtremalPoints];
467   // Look for the extremal points (myLExtremalPoints)
468   for (Standard_Integer anAxeInd = 0, aPrmInd = -1; anAxeInd < myNbInitAxes; ++anAxeInd)
469   {
470     Standard_Integer aMinInd = ++aPrmInd, aMaxInd = ++aPrmInd;
471     aParams[aMinInd] = RealLast();
472     aParams[aMaxInd] = -RealLast();
473     Project (anInitialAxesArray[anAxeInd],
474              aParams[aMinInd], aParams[aMaxInd],
475              &myLExtremalPoints[aMinInd], &myLExtremalPoints[aMaxInd]);
476   }
477
478   // For not optimal box it is necessary to compute the max axis
479   // created by the maximally distant extreme points
480   if (!myOptimal)
481   {
482     for(Standard_Integer i = 0; i < 5; i++)
483       myTriIdx[i] = INT_MAX;
484
485     // Compute myTriIdx[0] and myTriIdx[1].
486     Standard_Real aMaxSqDist = -1.0;
487     for (Standard_Integer aPrmInd = 0; aPrmInd < myNbExtremalPoints; aPrmInd += 2)
488     {
489       const gp_Pnt &aP1 = myLExtremalPoints[aPrmInd],
490         &aP2 = myLExtremalPoints[aPrmInd + 1];
491       const Standard_Real aSqDist = aP1.SquareDistance(aP2);
492       if (aSqDist > aMaxSqDist)
493       {
494         aMaxSqDist = aSqDist;
495         myTriIdx[0] = aPrmInd;
496         myTriIdx[1] = aPrmInd + 1;
497       }
498     }
499
500     // Compute the maximal axis orthogonal to the found one
501     FillToTriangle3();
502   }
503 }
504
505 //=======================================================================
506 // Function : FillToTriangle3
507 // purpose : Two value of myTriIdx array is known. Let us find myTriIdx[2].
508 //            It must be in maximal distance from the infinite axis going
509 //            through the points with indexes myTriIdx[0] and myTriIdx[1].
510 //=======================================================================
511 void OBBTool::FillToTriangle3()
512 {
513   const gp_XYZ &aP0 = myLExtremalPoints[myTriIdx[0]];
514   const gp_XYZ anAxis = myLExtremalPoints[myTriIdx[1]] - aP0;
515   Standard_Real aMaxSqDist = -1.0;
516   for(Standard_Integer i = 0; i < myNbExtremalPoints; i++)
517   {
518     if((i == myTriIdx[0]) || (i == myTriIdx[1]))
519       continue;
520
521     const gp_XYZ &aP = myLExtremalPoints[i];
522     const Standard_Real aDistToAxe = anAxis.CrossSquareMagnitude(aP - aP0);
523     if(aDistToAxe > aMaxSqDist)
524     {
525       myTriIdx[2] = i;
526       aMaxSqDist = aDistToAxe;
527     }
528   }
529 }
530
531 //=======================================================================
532 // Function : FillToTriangle5
533 // purpose : Three value of myTriIdx array is known.
534 //            Let us find myTriIdx[3] and myTriIdx[4].
535 //           They must be in the different sides of the plane of
536 //            triangle set by points myTriIdx[0], myTriIdx[1] and
537 //            myTriIdx[2]. Moreover, the distance from these points
538 //            to the triangle plane must be maximal.
539 //=======================================================================
540 void OBBTool::FillToTriangle5(const gp_XYZ& theNormal,
541                               const gp_XYZ& theBarryCenter)
542 {
543   Standard_Real aParams[2] = {0.0, 0.0};
544   Standard_Integer id3 = -1, id4 = -1;
545
546   for(Standard_Integer aPtIdx = 0; aPtIdx < myNbExtremalPoints; aPtIdx++)
547   {
548     if((aPtIdx == myTriIdx[0]) || (aPtIdx == myTriIdx[1]) || (aPtIdx == myTriIdx[2]))
549       continue;
550
551     const gp_XYZ &aCurrPoint = myLExtremalPoints[aPtIdx];
552     const Standard_Real aParam = theNormal.Dot(aCurrPoint - theBarryCenter);
553
554     if (aParam < aParams[0])
555     {
556       id3 = aPtIdx;
557       aParams[0] = aParam;
558     }
559     else if (aParam > aParams[1])
560     {
561       id4 = aPtIdx;
562       aParams[1] = aParam;
563     }
564   }
565
566   // The points must be in the different sides of the triangle plane.
567   if (id3 >= 0 && aParams[0] < -Precision::Confusion())
568     myTriIdx[3] = id3;
569
570   if (id4 >= 0 && aParams[1] > Precision::Confusion())
571     myTriIdx[4] = id4;
572 }
573
574 //=======================================================================
575 // Function : ProcessTriangle
576 // purpose : Choose the optimal box with triple axes containing normal
577 //            to the triangle and some edge of the triangle (3rd axis is
578 //            computed from these two ones).
579 //=======================================================================
580 void OBBTool::ProcessTriangle(const Standard_Integer theIdx1,
581                               const Standard_Integer theIdx2,
582                               const Standard_Integer theIdx3,
583                               const Standard_Boolean theIsBuiltTrg)
584 {
585   const Standard_Integer aNbAxes = 3;
586
587   // All axes must be normalized in order to provide correct area computation
588   // (see ComputeQuality(...) method).
589   int ID1[3] = { theIdx2, theIdx3, theIdx1 },
590       ID2[3] = { theIdx1, theIdx2, theIdx3 };
591   gp_XYZ aYAxis[aNbAxes] = {(myLExtremalPoints[ID1[0]] - myLExtremalPoints[ID2[0]]),
592                             (myLExtremalPoints[ID1[1]] - myLExtremalPoints[ID2[1]]),
593                             (myLExtremalPoints[ID1[2]] - myLExtremalPoints[ID2[2]])};
594
595   // Normal to the triangle plane
596   gp_XYZ aZAxis = aYAxis[0].Crossed(aYAxis[1]);
597
598   Standard_Real aSqMod = aZAxis.SquareModulus();
599
600   if (aSqMod < Precision::SquareConfusion())
601     return;
602
603   aZAxis /= Sqrt(aSqMod);
604
605   gp_XYZ aXAxis[aNbAxes];
606   for (Standard_Integer i = 0; i < aNbAxes; i++)
607     aXAxis[i] = aYAxis[i].Crossed(aZAxis).Normalized();
608
609   if (theIsBuiltTrg)
610     FillToTriangle5 (aZAxis, myLExtremalPoints[theIdx1]);
611
612   // Min and Max parameter
613   const Standard_Integer aNbPoints = 2 * aNbAxes;
614
615   // Compute Min/Max params for ZAxis
616   Standard_Real aParams[aNbPoints];
617   FindMinMax (aZAxis, aParams[4], aParams[5]); // Compute params on ZAxis once
618
619   Standard_Integer aMinIdx = -1;
620   for(Standard_Integer anAxeInd = 0; anAxeInd < aNbAxes; anAxeInd++)
621   {
622     const gp_XYZ &aAX = aXAxis[anAxeInd];
623     // Compute params on XAxis
624     FindMinMax (aAX, aParams[0], aParams[1]);
625     // Compute params on YAxis checking for stored values
626     ComputeParams (ID1[anAxeInd], ID2[anAxeInd], aParams[2], aParams[3]);
627
628     const Standard_Real anArea = ComputeQuality(aParams);
629     if (anArea < myQualityCriterion)
630     {
631       myQualityCriterion = anArea;
632       aMinIdx = anAxeInd;
633     }
634   }
635
636   if (aMinIdx < 0)
637     return;
638
639   myAxes[0] = aXAxis[aMinIdx];
640   myAxes[1] = aYAxis[aMinIdx].Normalized();
641   myAxes[2] = aZAxis;
642 }
643 //=======================================================================
644 // Function : ProcessDiTetrahedron
645 // purpose : DiTo-algorithm (http://www.idt.mdh.se/~tla/publ/FastOBBs.pdf)
646 //=======================================================================
647 void OBBTool::ProcessDiTetrahedron()
648 {
649   // To compute the optimal OBB it is necessary to check all possible
650   // axes created by the extremal points. It is also necessary to project
651   // all the points on the axis, as for each different axis there will be
652   // different extremal points.
653   if (myOptimal)
654   {
655     for (Standard_Integer i = 0; i < myNbExtremalPoints - 2; i++)
656     {
657       for (Standard_Integer j = i + 1; j < myNbExtremalPoints - 1; j++)
658       {
659         for (Standard_Integer k = j + 1; k < myNbExtremalPoints; k++)
660         {
661           ProcessTriangle (i, j, k, Standard_False);
662         }
663       }
664     }
665   }
666   else
667   {
668     // Use the standard DiTo approach
669     ProcessTriangle(myTriIdx[0], myTriIdx[1], myTriIdx[2], Standard_True);
670
671     if (myTriIdx[3] <= myNbExtremalPoints)
672     {
673       ProcessTriangle(myTriIdx[0], myTriIdx[1], myTriIdx[3], Standard_False);
674       ProcessTriangle(myTriIdx[1], myTriIdx[2], myTriIdx[3], Standard_False);
675       ProcessTriangle(myTriIdx[0], myTriIdx[2], myTriIdx[3], Standard_False);
676     }
677
678     if (myTriIdx[4] <= myNbExtremalPoints)
679     {
680       ProcessTriangle(myTriIdx[0], myTriIdx[1], myTriIdx[4], Standard_False);
681       ProcessTriangle(myTriIdx[1], myTriIdx[2], myTriIdx[4], Standard_False);
682       ProcessTriangle(myTriIdx[0], myTriIdx[2], myTriIdx[4], Standard_False);
683     }
684   }
685 }
686
687 //=======================================================================
688 // Function : BuildBox
689 // purpose : 
690 //=======================================================================
691 void OBBTool::BuildBox(Bnd_OBB& theBox)
692 {
693   theBox.SetVoid();
694
695   // In fact, use Precision::SquareConfusion().
696   const Standard_Boolean isOBB = myAxes[0].SquareModulus()*
697                                  myAxes[1].SquareModulus()*
698                                  myAxes[2].SquareModulus() > 1.0e-14;
699
700   const gp_Dir aXDir = isOBB ? myAxes[0] : gp_Dir(1, 0, 0);
701   const gp_Dir aYDir = isOBB ? myAxes[1] : gp_Dir(0, 1, 0);
702   const gp_Dir aZDir = isOBB ? myAxes[2] : gp_Dir(0, 0, 1);
703
704   const Standard_Integer aNbPoints = 6;
705   Standard_Real aParams[aNbPoints];
706
707   gp_XYZ aFCurrPoint = myPntsList.First().XYZ();
708   
709   aParams[0] = aParams[1] = aFCurrPoint.Dot(aXDir.XYZ());
710   aParams[2] = aParams[3] = aFCurrPoint.Dot(aYDir.XYZ());
711   aParams[4] = aParams[5] = aFCurrPoint.Dot(aZDir.XYZ());
712
713   if(myListOfTolers != 0)
714   {
715     const Standard_Real aTol = myListOfTolers->First();
716     aParams[0] -= aTol;
717     aParams[1] += aTol;
718     aParams[2] -= aTol;
719     aParams[3] += aTol;
720     aParams[4] -= aTol;
721     aParams[5] += aTol;
722   }
723
724   for(Standard_Integer i = myPntsList.Lower() + 1; i <= myPntsList.Upper(); i++)
725   {
726     const gp_XYZ &aCurrPoint = myPntsList(i).XYZ();
727     const Standard_Real aDx = aCurrPoint.Dot(aXDir.XYZ()),
728                         aDy = aCurrPoint.Dot(aYDir.XYZ()),
729                         aDz = aCurrPoint.Dot(aZDir.XYZ());
730         
731     if(myListOfTolers == 0)
732     {
733       SetMinMax(&aParams[0], aDx);
734       SetMinMax(&aParams[2], aDy);
735       SetMinMax(&aParams[4], aDz);
736     }
737     else
738     {
739       const Standard_Real aTol = myListOfTolers->Value(i);
740       aParams[0] = Min(aParams[0], aDx - aTol);
741       aParams[1] = Max(aParams[1], aDx + aTol);
742       aParams[2] = Min(aParams[2], aDy - aTol);
743       aParams[3] = Max(aParams[3], aDy + aTol);
744       aParams[4] = Min(aParams[4], aDz - aTol);
745       aParams[5] = Max(aParams[5], aDz + aTol);
746     }
747   }
748
749   //Half-sizes
750   const Standard_Real aHX = 0.5*(aParams[1] - aParams[0]);
751   const Standard_Real aHY = 0.5*(aParams[3] - aParams[2]);
752   const Standard_Real aHZ = 0.5*(aParams[5] - aParams[4]);
753
754   const gp_XYZ aCenter = 0.5*((aParams[1] + aParams[0])*aXDir.XYZ() +
755                               (aParams[3] + aParams[2])*aYDir.XYZ() +
756                               (aParams[5] + aParams[4])*aZDir.XYZ());
757
758   theBox.SetCenter(aCenter);
759   theBox.SetXComponent(aXDir, aHX);
760   theBox.SetYComponent(aYDir, aHY);
761   theBox.SetZComponent(aZDir, aHZ);
762   theBox.SetAABox(!isOBB);
763 }
764
765 // =======================================================================
766 // function : ReBuild
767 // purpose  : http://www.idt.mdh.se/~tla/publ/
768 // =======================================================================
769 void Bnd_OBB::ReBuild(const TColgp_Array1OfPnt& theListOfPoints,
770                       const TColStd_Array1OfReal *theListOfTolerances,
771                       const Standard_Boolean theIsOptimal)
772 {
773   switch(theListOfPoints.Length())
774   {
775     case 1:
776       ProcessOnePoint(theListOfPoints.First());
777       if(theListOfTolerances)
778         Enlarge(theListOfTolerances->First());
779       return;
780     case 2:
781     {
782       const Standard_Real aTol1 = (theListOfTolerances == 0) ? 0.0 :
783                                       theListOfTolerances->First();
784
785       const Standard_Real aTol2 = (theListOfTolerances == 0) ? 0.0 :
786                                       theListOfTolerances->Last();
787       
788       const gp_XYZ &aP1 = theListOfPoints.First().XYZ(),
789                    &aP2 = theListOfPoints.Last().XYZ();
790       const gp_XYZ aDP = aP2 - aP1;
791       const Standard_Real aDPm = aDP.Modulus();
792       myIsAABox = Standard_False;
793       myHDims[1] = myHDims[2] = Max(aTol1, aTol2);
794
795       if(aDPm < Precision::Confusion())
796       {
797         ProcessOnePoint(aP1);
798         Enlarge(myHDims[1] + Precision::Confusion());
799         return;
800       }
801
802       myHDims[0] = 0.5*(aDPm+aTol1+aTol2);
803       myAxes[0] = aDP/aDPm;
804       if(Abs(myAxes[0].X()) > Abs(myAxes[0].Y()))
805       {
806         // Z-coord. is maximal or X-coord. is maximal
807         myAxes[1].SetCoord(-myAxes[0].Z(), 0.0, myAxes[0].X());
808       }
809       else
810       {
811         // Z-coord. is maximal or Y-coord. is maximal
812         myAxes[1].SetCoord(0.0, -myAxes[0].Z(), myAxes[0].Y());
813       }
814
815       myAxes[2] = myAxes[0].Crossed(myAxes[1]).Normalized();
816       myCenter = aP1 + 0.5*(aDPm - aTol1 + aTol2)*myAxes[0];
817     }
818     return;
819     default:
820       break;
821   }
822
823   OBBTool aTool(theListOfPoints, theListOfTolerances, theIsOptimal);
824   aTool.ProcessDiTetrahedron();
825   aTool.BuildBox(*this);
826 }
827
828 // =======================================================================
829 // function : IsOut
830 // purpose  : 
831 // =======================================================================
832 Standard_Boolean Bnd_OBB::IsOut(const Bnd_OBB& theOther) const
833 {
834   if (IsVoid() || theOther.IsVoid())
835     return Standard_True;
836
837   if (myIsAABox && theOther.myIsAABox)
838   {
839     return ((Abs(theOther.myCenter.X() - myCenter.X()) > theOther.myHDims[0] + myHDims[0]) ||
840             (Abs(theOther.myCenter.Y() - myCenter.Y()) > theOther.myHDims[1] + myHDims[1]) ||
841             (Abs(theOther.myCenter.Z() - myCenter.Z()) > theOther.myHDims[2] + myHDims[2]));
842   }
843
844   // According to the Separating Axis Theorem for Oriented Bounding Boxes
845   // it is necessary to check the 15 separating axes (Ls):
846   // - 6 axes of the boxes;
847   // - 9 cross products of the axes of the boxes.
848   // If any of these axes is valid, the boxes do not interfere.
849
850   // The algorithm is following:
851   // 1. Compute the "length" for j-th BndBox (j=1...2) according to the formula:
852   //    L(j)=Sum(myHDims[i]*Abs(myAxes[i].Dot(Ls)))
853   // 2. If (theCenter2 - theCenter1).Dot(Ls) > (L(1) + L(2))
854   //    then the considered OBBs are not interfered in terms of the axis Ls.
855   //
856   // If OBBs are not interfered in terms of at least one axis (of 15) then
857   // they are not interfered at all.
858
859   // Precomputed difference between centers
860   gp_XYZ D = theOther.myCenter - myCenter;
861
862   // Check the axes of the this box, i.e. L is one of myAxes
863   // Since the Dot product of two of these directions is null, it could be skipped:
864   // myXDirection.Dot(myYDirection) = 0
865
866   for(Standard_Integer i = 0; i < 3; ++i)
867   {
868     // Length of the second segment
869     Standard_Real aLSegm2 = 0;
870     for(Standard_Integer j = 0; j < 3; ++j)
871       aLSegm2 += theOther.myHDims[j] * Abs(theOther.myAxes[j].Dot(myAxes[i]));
872
873     // Distance between projected centers
874     Standard_Real aDistCC = Abs(D.Dot(myAxes[i]));
875
876     if(aDistCC > myHDims[i] + aLSegm2)
877       return Standard_True;
878   }
879
880   // Check the axes of the Other box, i.e. L is one of theOther.myAxes
881
882   for(Standard_Integer i = 0; i < 3; ++i)
883   {
884     // Length of the first segment
885     Standard_Real aLSegm1 = 0.;
886     for(Standard_Integer j = 0; j < 3; ++j)
887       aLSegm1 += myHDims[j] * Abs(myAxes[j].Dot(theOther.myAxes[i]));
888
889     // Distance between projected centers
890     Standard_Real aDistCC = Abs(D.Dot(theOther.myAxes[i]));
891
892     if(aDistCC > aLSegm1 + theOther.myHDims[i])
893       return Standard_True;
894   }
895
896   const Standard_Real aTolNull = Epsilon(1.0);
897
898   // Check the axes produced by the cross products
899   for(Standard_Integer i = 0; i < 3; ++i)
900   {
901     for(Standard_Integer j = 0; j < 3; ++j)
902     {
903       // Separating axis
904       gp_XYZ aLAxe = myAxes[i].Crossed(theOther.myAxes[j]);
905
906       const Standard_Real aNorm = aLAxe.Modulus();
907       if(aNorm < aTolNull)
908         continue;
909
910       aLAxe /= aNorm;
911
912       // Length of the first segment
913       Standard_Real aLSegm1 = 0.;
914       for(Standard_Integer k = 0; k < 3; ++k)
915         aLSegm1 += myHDims[k] * Abs(myAxes[k].Dot(aLAxe));
916
917       // Length of the second segment
918       Standard_Real aLSegm2 = 0.;
919       for(Standard_Integer k = 0; k < 3; ++k)
920         aLSegm2 += theOther.myHDims[k] * Abs(theOther.myAxes[k].Dot(aLAxe));
921
922       // Distance between projected centers
923       Standard_Real aDistCC = Abs(D.Dot(aLAxe));
924
925       if(aDistCC > aLSegm1 + aLSegm2)
926         return Standard_True;
927     }
928   }
929
930   return Standard_False;
931 }
932
933 // =======================================================================
934 // function : IsOut
935 // purpose  : 
936 // =======================================================================
937 Standard_Boolean Bnd_OBB::IsOut(const gp_Pnt& theP) const
938 {
939   // 1. Project the point to myAxes[i] (i=0...2).
940   // 2. Check, whether the absolute value of the correspond
941   //    projection parameter is greater than myHDims[i].
942   //    In this case, IsOut method will return TRUE.
943
944   const gp_XYZ aRV = theP.XYZ() - myCenter;
945
946   return ((Abs(myAxes[0].Dot(aRV)) > myHDims[0]) || 
947           (Abs(myAxes[1].Dot(aRV)) > myHDims[1]) ||
948           (Abs(myAxes[2].Dot(aRV)) > myHDims[2]));
949 }
950
951 // =======================================================================
952 // function : IsCompletelyInside
953 // purpose  : Checks if every vertex of theOther is completely inside *this
954 // =======================================================================
955 Standard_Boolean Bnd_OBB::IsCompletelyInside(const Bnd_OBB& theOther) const
956 {
957   if(IsVoid() || theOther.IsVoid())
958     return Standard_False;
959
960   gp_Pnt aVert[8];
961   theOther.GetVertex(aVert);
962   for(Standard_Integer i = 0; i < 8; i++)
963   {
964     if(IsOut(aVert[i]))
965       return Standard_False;
966   }
967
968   return Standard_True;
969 }
970
971 // =======================================================================
972 // function : Add
973 // purpose  : 
974 // =======================================================================
975 void Bnd_OBB::Add(const gp_Pnt& theP)
976 {
977   if (IsVoid())
978   {
979     myCenter = theP.XYZ();
980     myAxes[0] = gp::DX().XYZ();
981     myAxes[1] = gp::DY().XYZ();
982     myAxes[2] = gp::DZ().XYZ();
983     myHDims[0] = 0.0;
984     myHDims[1] = 0.0;
985     myHDims[2] = 0.0;
986     myIsAABox = Standard_True;
987   }
988   else
989   {
990     gp_Pnt aList[9];
991     GetVertex(aList);
992     aList[8] = theP;
993     ReBuild(TColgp_Array1OfPnt(aList[0], 0, 8));
994   }
995 }
996
997 // =======================================================================
998 // function : Add
999 // purpose  : 
1000 // =======================================================================
1001 void Bnd_OBB::Add(const Bnd_OBB& theOther)
1002 {
1003   if (!theOther.IsVoid())
1004   {
1005     if (IsVoid())
1006     {
1007       myCenter = theOther.myCenter;
1008       myAxes[0] = theOther.myAxes[0];
1009       myAxes[1] = theOther.myAxes[1];
1010       myAxes[2] = theOther.myAxes[2];
1011       myHDims[0] = theOther.myHDims[0];
1012       myHDims[1] = theOther.myHDims[1];
1013       myHDims[2] = theOther.myHDims[2];
1014       myIsAABox = theOther.myIsAABox;
1015     }
1016     else
1017     {
1018       gp_Pnt aList[16];
1019       GetVertex(&aList[0]);
1020       theOther.GetVertex(&aList[8]);
1021       ReBuild(TColgp_Array1OfPnt(aList[0], 0, 15));
1022     }
1023   }
1024 }
1025
1026 //=======================================================================
1027 //function : DumpJson
1028 //purpose  : 
1029 //=======================================================================
1030 void Bnd_OBB::DumpJson (Standard_OStream& theOStream, const Standard_Integer theDepth) const
1031 {
1032   OCCT_DUMP_CLASS_BEGIN (theOStream, Bnd_OBB);
1033
1034   OCCT_DUMP_FIELD_VALUES_DUMPED (theOStream, theDepth, &myCenter);
1035   OCCT_DUMP_FIELD_VALUES_DUMPED (theOStream, theDepth, &myAxes[0]);
1036   OCCT_DUMP_FIELD_VALUES_DUMPED (theOStream, theDepth, &myAxes[1]);
1037   OCCT_DUMP_FIELD_VALUES_DUMPED (theOStream, theDepth, &myAxes[2]);
1038
1039   OCCT_DUMP_FIELD_VALUE_NUMERICAL (theOStream, myHDims[0]);
1040   OCCT_DUMP_FIELD_VALUE_NUMERICAL (theOStream, myHDims[1]);
1041   OCCT_DUMP_FIELD_VALUE_NUMERICAL (theOStream, myHDims[2]);
1042   OCCT_DUMP_FIELD_VALUE_NUMERICAL (theOStream, myIsAABox);
1043 }