f2600de4c897de1dcc26f70c5a5b290949d44e7b
[occt.git] / src / BRepExtrema / BRepExtrema_OverlapTool.cxx
1 // Created on: 2015-04-26
2 // Created by: Denis BOGOLEPOV
3 // Copyright (c) 2015 OPEN CASCADE SAS
4 //
5 // This file is part of Open CASCADE Technology software library.
6 //
7 // This library is free software; you can redistribute it and/or modify it under
8 // the terms of the GNU Lesser General Public License version 2.1 as published
9 // by the Free Software Foundation, with special exception defined in the file
10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
11 // distribution for complete text of the license and disclaimer of any warranty.
12 //
13 // Alternatively, this file may be used under the terms of Open CASCADE
14 // commercial license or contractual agreement.
15
16 #include <Precision.hxx>
17
18 #include <BRepExtrema_OverlapTool.hxx>
19
20 //=======================================================================
21 //function : BRepExtrema_OverlapTool
22 //purpose  :
23 //=======================================================================
24 BRepExtrema_OverlapTool::BRepExtrema_OverlapTool()
25 : myFilter (NULL)
26 {
27   myIsDone = Standard_False;
28 }
29
30 //=======================================================================
31 //function : BRepExtrema_OverlapTool
32 //purpose  :
33 //=======================================================================
34 BRepExtrema_OverlapTool::BRepExtrema_OverlapTool (const Handle(BRepExtrema_TriangleSet)& theSet1,
35                                                   const Handle(BRepExtrema_TriangleSet)& theSet2)
36 : myFilter (NULL)
37 {
38   LoadTriangleSets (theSet1, theSet2);
39 }
40
41 //=======================================================================
42 //function : LoadTriangleSets
43 //purpose  :
44 //=======================================================================
45 void BRepExtrema_OverlapTool::LoadTriangleSets (const Handle(BRepExtrema_TriangleSet)& theSet1,
46                                                 const Handle(BRepExtrema_TriangleSet)& theSet2)
47 {
48   mySet1 = theSet1;
49   mySet2 = theSet2;
50
51   myIsDone = Standard_False;
52 }
53
54 #ifndef DBL_EPSILON
55   #define DBL_EPSILON std::numeric_limits<Standard_Real>::epsilon()
56 #endif
57
58 namespace
59 {
60   //! Tool class to describe stack item in traverse function.
61   struct BRepExtrema_StackItem
62   {
63     Standard_Integer Node1;
64     Standard_Integer Node2;
65
66     BRepExtrema_StackItem (const Standard_Integer theNode1 = 0,
67                            const Standard_Integer theNode2 = 0)
68     : Node1 (theNode1),
69       Node2 (theNode2)
70     {
71       //
72     }
73   };
74
75   //! Bounding triangular prism for specified triangle.
76   class BRepExtrema_BoundingPrism
77   {
78   public:
79
80     //! Vertices of the prism.
81     BVH_Vec3d Vertices[6];
82
83     //! Edges of the prism.
84     BVH_Vec3d Edges[3];
85
86     //! Normal to prism caps.
87     BVH_Vec3d Normal;
88
89     //! Normals to prism edges.
90     BVH_Vec3d EdgeNormals[3];
91
92     //! Is prism initialized?
93     Standard_Boolean IsInited;
94
95   public:
96
97     //! Creates uninitialized bounding prism.
98     BRepExtrema_BoundingPrism() : IsInited (Standard_False)
99     {
100       //
101     }
102
103     //! Creates new bounding prism for the given triangle.
104     BRepExtrema_BoundingPrism (const BVH_Vec3d&    theVertex0,
105                                const BVH_Vec3d&    theVertex1,
106                                const BVH_Vec3d&    theVertex2,
107                                const Standard_Real theDeflect)
108     {
109       Init (theVertex0,
110             theVertex1,
111             theVertex2,
112             theDeflect);
113     }
114
115     //! Calculates bounding prism for the given triangle.
116     void Init (const BVH_Vec3d&    theVertex0,
117                const BVH_Vec3d&    theVertex1,
118                const BVH_Vec3d&    theVertex2,
119                const Standard_Real theDeflect)
120     {
121       Edges[0] = theVertex1 - theVertex0;
122       Edges[1] = theVertex2 - theVertex0;
123       Edges[2] = theVertex2 - theVertex1;
124
125       Normal = BVH_Vec3d::Cross (Edges[0], Edges[1]);
126
127       EdgeNormals[0] = BVH_Vec3d::Cross (Edges[0], Normal);
128       EdgeNormals[1] = BVH_Vec3d::Cross (Edges[1], Normal);
129       EdgeNormals[2] = BVH_Vec3d::Cross (Edges[2], Normal);
130
131       EdgeNormals[0] *= 1.0 / Max (EdgeNormals[0].Modulus(), Precision::Confusion());
132       EdgeNormals[1] *= 1.0 / Max (EdgeNormals[1].Modulus(), Precision::Confusion());
133       EdgeNormals[2] *= 1.0 / Max (EdgeNormals[2].Modulus(), Precision::Confusion());
134
135       const BVH_Vec3d aDirect01 = EdgeNormals[0] - EdgeNormals[1];
136       const BVH_Vec3d aDirect02 = EdgeNormals[0] + EdgeNormals[2];
137       const BVH_Vec3d aDirect12 = EdgeNormals[2] - EdgeNormals[1];
138
139       Vertices[0] = Vertices[3] = theVertex0 + aDirect01 * (theDeflect / aDirect01.Dot (EdgeNormals[0]));
140       Vertices[1] = Vertices[4] = theVertex1 + aDirect02 * (theDeflect / aDirect02.Dot (EdgeNormals[2]));
141       Vertices[2] = Vertices[5] = theVertex2 + aDirect12 * (theDeflect / aDirect12.Dot (EdgeNormals[2]));
142
143       const BVH_Vec3d aNormOffset = Normal * (theDeflect / Max (Normal.Modulus(), Precision::Confusion()));
144
145       for (Standard_Integer aVertIdx = 0; aVertIdx < 3; ++aVertIdx)
146       {
147         Vertices[aVertIdx + 0] += aNormOffset;
148         Vertices[aVertIdx + 3] -= aNormOffset;
149       }
150
151       IsInited = Standard_True;
152     }
153
154     //! Checks if two prisms are separated along the given axis.
155     Standard_Boolean Separated (const BRepExtrema_BoundingPrism& thePrism, const BVH_Vec3d& theAxis) const
156     {
157       Standard_Real aMin1 =  DBL_MAX;
158       Standard_Real aMax1 = -DBL_MAX;
159
160       Standard_Real aMin2 =  DBL_MAX;
161       Standard_Real aMax2 = -DBL_MAX;
162
163       for (Standard_Integer aVertIdx = 0; aVertIdx < 6; ++aVertIdx)
164       {
165         const Standard_Real aProj1 = Vertices[aVertIdx].Dot (theAxis);
166
167         aMin1 = Min (aMin1, aProj1);
168         aMax1 = Max (aMax1, aProj1);
169
170         const Standard_Real aProj2 = thePrism.Vertices[aVertIdx].Dot (theAxis);
171
172         aMin2 = Min (aMin2, aProj2);
173         aMax2 = Max (aMax2, aProj2);
174
175         if (aMin1 <= aMax2 && aMax1 >= aMin2)
176         {
177           return Standard_False;
178         }
179       }
180
181       return aMin1 > aMax2 || aMax1 < aMin2;
182     }
183   };
184
185   // =======================================================================
186   // function : sign
187   // purpose  :
188   // =======================================================================
189   Standard_Real sign (const BVH_Vec3d& theVertex0,
190                       const BVH_Vec3d& theVertex1,
191                       const BVH_Vec3d& theVertex2,
192                       const Standard_Integer theX,
193                       const Standard_Integer theY)
194   {
195     return (theVertex0[theX] - theVertex2[theX]) * (theVertex1[theY] - theVertex2[theY]) -
196            (theVertex1[theX] - theVertex2[theX]) * (theVertex0[theY] - theVertex2[theY]);
197   }
198
199   // =======================================================================
200   // function : pointInTriangle
201   // purpose  :
202   // =======================================================================
203   Standard_Boolean pointInTriangle (const BVH_Vec3d& theTestPnt,
204                                     const BVH_Vec3d& theTrgVtx0,
205                                     const BVH_Vec3d& theTrgVtx1,
206                                     const BVH_Vec3d& theTrgVtx2,
207                                     const Standard_Integer theX,
208                                     const Standard_Integer theY)
209   {
210     const Standard_Boolean aSign0 = sign (theTestPnt, theTrgVtx0, theTrgVtx1, theX, theY) <= 0.0;
211     const Standard_Boolean aSign1 = sign (theTestPnt, theTrgVtx1, theTrgVtx2, theX, theY) <= 0.0;
212     const Standard_Boolean aSign2 = sign (theTestPnt, theTrgVtx2, theTrgVtx0, theX, theY) <= 0.0;
213
214     return (aSign0 == aSign1) && (aSign1 == aSign2);
215   }
216
217   // =======================================================================
218   // function : segmentsIntersected
219   // purpose  : Checks if two line segments are intersected
220   // =======================================================================
221   Standard_Boolean segmentsIntersected (const BVH_Vec2d& theOriginSeg0,
222                                         const BVH_Vec2d& theOriginSeg1,
223                                         const BVH_Vec2d& theDirectSeg0,
224                                         const BVH_Vec2d& theDirectSeg1)
225   {
226     const Standard_Real aDet = -theDirectSeg1.x() * theDirectSeg0.y() +
227                                 theDirectSeg0.x() * theDirectSeg1.y();
228
229     if (fabs (aDet) < DBL_EPSILON) // segments are parallel
230     {
231       const BVH_Vec2d aDirect = theDirectSeg0 * (1.0 / theDirectSeg0.Modulus());
232
233       const Standard_Real aEdge0Time0 = theOriginSeg0.Dot (aDirect);
234       const Standard_Real aEdge1Time0 = theOriginSeg1.Dot (aDirect);
235
236       const Standard_Real aEdge0Time1 = aEdge0Time0 + theDirectSeg0.Dot (aDirect);
237       const Standard_Real aEdge1Time1 = aEdge1Time0 + theDirectSeg1.Dot (aDirect);
238
239       const Standard_Real aEdge0Min = Min (aEdge0Time0, aEdge0Time1);
240       const Standard_Real aEdge1Min = Min (aEdge1Time0, aEdge1Time1);
241       const Standard_Real aEdge0Max = Max (aEdge0Time0, aEdge0Time1);
242       const Standard_Real aEdge1Max = Max (aEdge1Time0, aEdge1Time1);
243
244       if (Max (aEdge0Min, aEdge1Min) > Min (aEdge0Max, aEdge1Max))
245       {
246         return Standard_False;
247       }
248
249       const BVH_Vec2d aNormal (-aDirect.y(), aDirect.x());
250
251       return fabs (theOriginSeg0.Dot (aNormal) - theOriginSeg1.Dot (aNormal)) < DBL_EPSILON;
252     }
253
254     const BVH_Vec2d aDelta = theOriginSeg0 - theOriginSeg1;
255
256     const Standard_Real aU = (-theDirectSeg0.y() * aDelta.x() + theDirectSeg0.x() * aDelta.y()) / aDet;
257     const Standard_Real aV = ( theDirectSeg1.x() * aDelta.y() - theDirectSeg1.y() * aDelta.x()) / aDet;
258
259     return aU >= 0.0 && aU <= 1.0 && aV >= 0.0 && aV <= 1.0;
260   }
261
262   // =======================================================================
263   // function : trianglesIntersected
264   // purpose  : Checks if two triangles are intersected
265   //            ("A Fast Triangle-Triangle Intersection Test" by T. Moller)
266   // =======================================================================
267   Standard_Boolean trianglesIntersected (const BVH_Vec3d& theTrng0Vert0,
268                                          const BVH_Vec3d& theTrng0Vert1,
269                                          const BVH_Vec3d& theTrng0Vert2,
270                                          const BVH_Vec3d& theTrng1Vert0,
271                                          const BVH_Vec3d& theTrng1Vert1,
272                                          const BVH_Vec3d& theTrng1Vert2)
273   {
274     const BVH_Vec3d aTrng1Normal = BVH_Vec3d::Cross (theTrng1Vert1 - theTrng1Vert0,
275                                                      theTrng1Vert2 - theTrng1Vert0).Normalized();
276
277     const Standard_Real aTrng1PlaneDist = aTrng1Normal.Dot (-theTrng1Vert0);
278
279     Standard_Real aDistTrng0Vert0 = aTrng1Normal.Dot (theTrng0Vert0) + aTrng1PlaneDist;
280     Standard_Real aDistTrng0Vert1 = aTrng1Normal.Dot (theTrng0Vert1) + aTrng1PlaneDist;
281     Standard_Real aDistTrng0Vert2 = aTrng1Normal.Dot (theTrng0Vert2) + aTrng1PlaneDist;
282
283     if ((aDistTrng0Vert0 < 0.0 && aDistTrng0Vert1 < 0.0 && aDistTrng0Vert2 < 0.0)
284      || (aDistTrng0Vert0 > 0.0 && aDistTrng0Vert1 > 0.0 && aDistTrng0Vert2 > 0.0))
285     {
286       return Standard_False; // 1st triangle lies on one side of the 2nd triangle
287     }
288
289     if (fabs (aDistTrng0Vert0) > Precision::Confusion()
290      || fabs (aDistTrng0Vert1) > Precision::Confusion()
291      || fabs (aDistTrng0Vert2) > Precision::Confusion()) // general 3D case
292     {
293       const BVH_Vec3d aTrng0Normal = BVH_Vec3d::Cross (theTrng0Vert1 - theTrng0Vert0,
294                                                        theTrng0Vert2 - theTrng0Vert0).Normalized();
295
296       const Standard_Real aTrng0PlaneDist = aTrng0Normal.Dot (-theTrng0Vert0);
297
298       Standard_Real aDistTrng1Vert0 = aTrng0Normal.Dot (theTrng1Vert0) + aTrng0PlaneDist;
299       Standard_Real aDistTrng1Vert1 = aTrng0Normal.Dot (theTrng1Vert1) + aTrng0PlaneDist;
300       Standard_Real aDistTrng1Vert2 = aTrng0Normal.Dot (theTrng1Vert2) + aTrng0PlaneDist;
301
302       if ((aDistTrng1Vert0 < 0.0 && aDistTrng1Vert1 < 0.0 && aDistTrng1Vert2 < 0.0)
303        || (aDistTrng1Vert0 > 0.0 && aDistTrng1Vert1 > 0.0 && aDistTrng1Vert2 > 0.0))
304       {
305         return Standard_False; // 2nd triangle lies on one side of the 1st triangle
306       }
307
308       const BVH_Vec3d aCrossLine = BVH_Vec3d::Cross (aTrng0Normal,
309                                                      aTrng1Normal);
310
311       Standard_Real aProjTrng0Vert0 = theTrng0Vert0.Dot (aCrossLine);
312       Standard_Real aProjTrng0Vert1 = theTrng0Vert1.Dot (aCrossLine);
313       Standard_Real aProjTrng0Vert2 = theTrng0Vert2.Dot (aCrossLine);
314
315       if (aDistTrng0Vert0 * aDistTrng0Vert1 > 0.0)
316       {
317         std::swap (aDistTrng0Vert1, aDistTrng0Vert2);
318         std::swap (aProjTrng0Vert1, aProjTrng0Vert2);
319       }
320       else if (aDistTrng0Vert1 * aDistTrng0Vert2 > 0.0)
321       {
322         std::swap (aDistTrng0Vert1, aDistTrng0Vert0);
323         std::swap (aProjTrng0Vert1, aProjTrng0Vert0);
324       }
325
326       Standard_Real aTime1 = fabs (aDistTrng0Vert0) <= DBL_EPSILON ? aProjTrng0Vert0 :
327         aProjTrng0Vert0 + (aProjTrng0Vert1 - aProjTrng0Vert0) * aDistTrng0Vert0 / (aDistTrng0Vert0 - aDistTrng0Vert1);
328       Standard_Real aTime2 = fabs (aDistTrng0Vert2) <= DBL_EPSILON ? aProjTrng0Vert2 :
329         aProjTrng0Vert2 + (aProjTrng0Vert1 - aProjTrng0Vert2) * aDistTrng0Vert2 / (aDistTrng0Vert2 - aDistTrng0Vert1);
330
331       const Standard_Real aTimeMin1 = Min (aTime1, aTime2);
332       const Standard_Real aTimeMax1 = Max (aTime1, aTime2);
333
334       Standard_Real aProjTrng1Vert0 = theTrng1Vert0.Dot (aCrossLine);
335       Standard_Real aProjTrng1Vert1 = theTrng1Vert1.Dot (aCrossLine);
336       Standard_Real aProjTrng1Vert2 = theTrng1Vert2.Dot (aCrossLine);
337
338       if (aDistTrng1Vert0 * aDistTrng1Vert1 > 0.0)
339       {
340         std::swap (aDistTrng1Vert1, aDistTrng1Vert2);
341         std::swap (aProjTrng1Vert1, aProjTrng1Vert2);
342       }
343       else if (aDistTrng1Vert1 * aDistTrng1Vert2 > 0.0)
344       {
345         std::swap (aDistTrng1Vert1, aDistTrng1Vert0);
346         std::swap (aProjTrng1Vert1, aProjTrng1Vert0);
347       }
348
349       aTime1 = fabs (aDistTrng1Vert0) <= DBL_EPSILON ? aProjTrng1Vert0 :
350         aProjTrng1Vert0 + (aProjTrng1Vert1 - aProjTrng1Vert0) * aDistTrng1Vert0 / (aDistTrng1Vert0 - aDistTrng1Vert1);
351       aTime2 = fabs (aDistTrng1Vert2) <= DBL_EPSILON ? aProjTrng1Vert2 :
352         aProjTrng1Vert2 + (aProjTrng1Vert1 - aProjTrng1Vert2) * aDistTrng1Vert2 / (aDistTrng1Vert2 - aDistTrng1Vert1);
353
354       const Standard_Real aTimeMin2 = Min (aTime1, aTime2);
355       const Standard_Real aTimeMax2 = Max (aTime1, aTime2);
356
357       aTime1 = Max (aTimeMin1, aTimeMin2);
358       aTime2 = Min (aTimeMax1, aTimeMax2);
359
360       return aTime1 <= aTime2; // intervals intersected --> triangles overlapped
361     }
362     else // triangles are co-planar
363     {
364       Standard_Integer anX;
365       Standard_Integer anY;
366
367       if (fabs (aTrng1Normal[0]) > fabs (aTrng1Normal[1]))
368       {
369         anX = fabs (aTrng1Normal[0]) > fabs (aTrng1Normal[2]) ? 1 : 0;
370         anY = fabs (aTrng1Normal[0]) > fabs (aTrng1Normal[2]) ? 2 : 1;
371       }
372       else
373       {
374         anX = fabs (aTrng1Normal[1]) > fabs (aTrng1Normal[2]) ? 0 : 0;
375         anY = fabs (aTrng1Normal[1]) > fabs (aTrng1Normal[2]) ? 2 : 1;
376       }
377
378       const BVH_Vec2d aOriginSeg0 [] = {BVH_Vec2d (theTrng0Vert0[anX], theTrng0Vert0[anY]),
379                                         BVH_Vec2d (theTrng0Vert1[anX], theTrng0Vert1[anY]),
380                                         BVH_Vec2d (theTrng0Vert2[anX], theTrng0Vert2[anY]) };
381
382       const BVH_Vec2d aDirectSeg0 [] = {aOriginSeg0[1] - aOriginSeg0[0],
383                                         aOriginSeg0[2] - aOriginSeg0[1],
384                                         aOriginSeg0[0] - aOriginSeg0[2] };
385
386       const BVH_Vec2d aOriginSeg1 [] = {BVH_Vec2d (theTrng1Vert0[anX], theTrng1Vert0[anY]),
387                                         BVH_Vec2d (theTrng1Vert1[anX], theTrng1Vert1[anY]),
388                                         BVH_Vec2d (theTrng1Vert2[anX], theTrng1Vert2[anY]) };
389
390       const BVH_Vec2d aDirectSeg1 [] = {aOriginSeg1[1] - aOriginSeg1[0],
391                                         aOriginSeg1[2] - aOriginSeg1[1],
392                                         aOriginSeg1[0] - aOriginSeg1[2] };
393
394       for (Standard_Integer aTrg0Edge = 0; aTrg0Edge < 3; ++aTrg0Edge)
395       {
396         for (Standard_Integer aTrg1Edge = 0; aTrg1Edge < 3; ++aTrg1Edge)
397         {
398           if (segmentsIntersected (aOriginSeg0[aTrg0Edge],
399                                    aOriginSeg1[aTrg1Edge],
400                                    aDirectSeg0[aTrg0Edge],
401                                    aDirectSeg1[aTrg1Edge]))
402           {
403             return Standard_True; // edges intersected --> triangles overlapped
404           }
405         }
406       }
407
408       if (pointInTriangle (theTrng1Vert0,
409                            theTrng0Vert0,
410                            theTrng0Vert1,
411                            theTrng0Vert2,
412                            anX,
413                            anY))
414       {
415         return Standard_True; // 1st triangle inside 2nd --> triangles overlapped
416       }
417
418       if (pointInTriangle (theTrng0Vert0,
419                            theTrng1Vert0,
420                            theTrng1Vert1,
421                            theTrng1Vert2,
422                            anX,
423                            anY))
424       {
425         return Standard_True; // 2nd triangle inside 1st --> triangles overlapped
426       }
427     }
428
429     return Standard_False;
430   }
431
432   // =======================================================================
433   // function : prismsIntersected
434   // purpose  : Checks if two triangular prisms are intersected
435   //            (test uses SAT - Separating Axis Theorem)
436   // =======================================================================
437   Standard_Boolean prismsIntersected (const BRepExtrema_BoundingPrism& thePrism1,
438                                       const BRepExtrema_BoundingPrism& thePrism2)
439   {
440     if (thePrism1.Separated (thePrism2, thePrism1.Normal))
441     {
442       return Standard_False;
443     }
444
445     if (thePrism1.Separated (thePrism2, thePrism2.Normal))
446     {
447       return Standard_False;
448     }
449
450     for (Standard_Integer anIdx = 0; anIdx < 3; ++anIdx)
451     {
452       if (thePrism1.Separated (thePrism2, thePrism1.EdgeNormals[anIdx]))
453       {
454         return Standard_False;
455       }
456     }
457
458     for (Standard_Integer anIdx = 0; anIdx < 3; ++anIdx)
459     {
460       if (thePrism1.Separated (thePrism2, thePrism2.EdgeNormals[anIdx]))
461       {
462         return Standard_False;
463       }
464     }
465
466     for (Standard_Integer anIdx1 = 0; anIdx1 < 4; ++anIdx1)
467     {
468       const BVH_Vec3d& aEdge1 = (anIdx1 == 3) ? thePrism1.Normal : thePrism1.Edges[anIdx1];
469
470       for (Standard_Integer anIdx2 = 0; anIdx2 < 4; ++anIdx2)
471       {
472         const BVH_Vec3d& aEdge2 = (anIdx2 == 3) ? thePrism2.Normal : thePrism2.Edges[anIdx2];
473
474         if (thePrism1.Separated (thePrism2, BVH_Vec3d::Cross (aEdge1, aEdge2)))
475         {
476           return Standard_False;
477         }
478       }
479     }
480
481     return Standard_True;
482   }
483
484   // =======================================================================
485   // function : overlapBoxes
486   // purpose  : Checks if two boxes (AABBs) are overlapped
487   // =======================================================================
488   inline Standard_Boolean overlapBoxes (const BVH_Vec3d&    theBoxMin1,
489                                         const BVH_Vec3d&    theBoxMax1,
490                                         const BVH_Vec3d&    theBoxMin2,
491                                         const BVH_Vec3d&    theBoxMax2,
492                                         const Standard_Real theTolerance)
493   {
494     // Check for overlap
495     return !(theBoxMin1.x() > theBoxMax2.x() + theTolerance ||
496              theBoxMax1.x() < theBoxMin2.x() - theTolerance ||
497              theBoxMin1.y() > theBoxMax2.y() + theTolerance ||
498              theBoxMax1.y() < theBoxMin2.y() - theTolerance ||
499              theBoxMin1.z() > theBoxMax2.z() + theTolerance ||
500              theBoxMax1.z() < theBoxMin2.z() - theTolerance);
501   }
502
503   //=======================================================================
504   //function : getSetOfFaces
505   //purpose  :
506   //=======================================================================
507   TColStd_PackedMapOfInteger& getSetOfFaces (
508     BRepExtrema_MapOfIntegerPackedMapOfInteger& theFaces, const Standard_Integer theFaceIdx)
509   {
510     if (!theFaces.IsBound (theFaceIdx))
511     {
512       theFaces.Bind (theFaceIdx, TColStd_PackedMapOfInteger());
513     }
514
515     return theFaces.ChangeFind (theFaceIdx);
516   }
517 }
518
519 //=======================================================================
520 //function : intersectTriangleRangesExact
521 //purpose  :
522 //=======================================================================
523 void BRepExtrema_OverlapTool::intersectTriangleRangesExact (const BVH_Vec4i& theLeaf1,
524                                                             const BVH_Vec4i& theLeaf2)
525 {
526   for (Standard_Integer aTrgIdx1 = theLeaf1.y(); aTrgIdx1 <= theLeaf1.z(); ++aTrgIdx1)
527   {
528     const Standard_Integer aFaceIdx1 = mySet1->GetFaceID (aTrgIdx1);
529
530     BVH_Vec3d aTrg1Vert1;
531     BVH_Vec3d aTrg1Vert2;
532     BVH_Vec3d aTrg1Vert3;
533
534     mySet1->GetVertices (aTrgIdx1,
535                          aTrg1Vert1,
536                          aTrg1Vert2,
537                          aTrg1Vert3);
538
539     const Standard_Boolean aIsInSet = myOverlapSubShapes1.IsBound (aFaceIdx1);
540
541     for (Standard_Integer aTrgIdx2 = theLeaf2.y(); aTrgIdx2 <= theLeaf2.z(); ++aTrgIdx2)
542     {
543       const Standard_Integer aFaceIdx2 = mySet2->GetFaceID (aTrgIdx2);
544
545       if (aIsInSet && myOverlapSubShapes1.Find (aFaceIdx1).Contains (aFaceIdx2))
546       {
547         continue;
548       }
549
550       BRepExtrema_ElementFilter::FilterResult aResult = myFilter == NULL ?
551         BRepExtrema_ElementFilter::DoCheck : myFilter->PreCheckElements (aTrgIdx1, aTrgIdx2);
552
553       if (aResult == BRepExtrema_ElementFilter::Overlap)
554       {
555         getSetOfFaces (myOverlapSubShapes1, aFaceIdx1).Add (aFaceIdx2);
556         getSetOfFaces (myOverlapSubShapes2, aFaceIdx2).Add (aFaceIdx1);
557
558 #ifdef OVERLAP_TOOL_OUTPUT_TRIANGLES
559         if (mySet1 == mySet2)
560         {
561           myOverlapTriangles1.Add (aTrgIdx1);
562           myOverlapTriangles1.Add (aTrgIdx2);
563         }
564         else
565         {
566           myOverlapTriangles1.Add (aTrgIdx1);
567           myOverlapTriangles2.Add (aTrgIdx2);
568         }
569 #endif
570       }
571       else if (aResult == BRepExtrema_ElementFilter::DoCheck)
572       {
573         BVH_Vec3d aTrg2Vert1;
574         BVH_Vec3d aTrg2Vert2;
575         BVH_Vec3d aTrg2Vert3;
576
577         mySet2->GetVertices (aTrgIdx2, aTrg2Vert1, aTrg2Vert2, aTrg2Vert3);
578
579         if (trianglesIntersected (aTrg1Vert1,
580                                   aTrg1Vert2,
581                                   aTrg1Vert3,
582                                   aTrg2Vert1,
583                                   aTrg2Vert2,
584                                   aTrg2Vert3))
585         {
586           getSetOfFaces (myOverlapSubShapes1, aFaceIdx1).Add (aFaceIdx2);
587           getSetOfFaces (myOverlapSubShapes2, aFaceIdx2).Add (aFaceIdx1);
588
589 #ifdef OVERLAP_TOOL_OUTPUT_TRIANGLES
590         if (mySet1 == mySet2)
591         {
592           myOverlapTriangles1.Add (aTrgIdx1);
593           myOverlapTriangles1.Add (aTrgIdx2);
594         }
595         else
596         {
597           myOverlapTriangles1.Add (aTrgIdx1);
598           myOverlapTriangles2.Add (aTrgIdx2);
599         }
600 #endif
601         }
602       }
603     }
604   }
605 }
606
607 //=======================================================================
608 //function : intersectTriangleRangesToler
609 //purpose  :
610 //=======================================================================
611 void BRepExtrema_OverlapTool::intersectTriangleRangesToler (const BVH_Vec4i&    theLeaf1,
612                                                             const BVH_Vec4i&    theLeaf2,
613                                                             const Standard_Real theToler)
614 {
615   for (Standard_Integer aTrgIdx1 = theLeaf1.y(); aTrgIdx1 <= theLeaf1.z(); ++aTrgIdx1)
616   {
617     const Standard_Integer aFaceIdx1 = mySet1->GetFaceID (aTrgIdx1);
618
619     BVH_Vec3d aTrg1Vert1;
620     BVH_Vec3d aTrg1Vert2;
621     BVH_Vec3d aTrg1Vert3;
622
623     mySet1->GetVertices (aTrgIdx1,
624                          aTrg1Vert1,
625                          aTrg1Vert2,
626                          aTrg1Vert3);
627
628     BRepExtrema_BoundingPrism aPrism1; // not initialized
629
630     const Standard_Boolean aIsInSet = myOverlapSubShapes1.IsBound (aFaceIdx1);
631
632     for (Standard_Integer aTrgIdx2 = theLeaf2.y(); aTrgIdx2 <= theLeaf2.z(); ++aTrgIdx2)
633     {
634       const Standard_Integer aFaceIdx2 = mySet2->GetFaceID (aTrgIdx2);
635
636       if (aIsInSet && myOverlapSubShapes1.Find (aFaceIdx1).Contains (aFaceIdx2))
637       {
638         continue;
639       }
640
641       BRepExtrema_ElementFilter::FilterResult aResult = myFilter == NULL ?
642         BRepExtrema_ElementFilter::DoCheck : myFilter->PreCheckElements (aTrgIdx1, aTrgIdx2);
643
644       if (aResult == BRepExtrema_ElementFilter::Overlap)
645       {
646         getSetOfFaces (myOverlapSubShapes1, aFaceIdx1).Add (aFaceIdx2);
647         getSetOfFaces (myOverlapSubShapes2, aFaceIdx2).Add (aFaceIdx1);
648
649 #ifdef OVERLAP_TOOL_OUTPUT_TRIANGLES
650         if (mySet1 == mySet2)
651         {
652           myOverlapTriangles1.Add (aTrgIdx1);
653           myOverlapTriangles1.Add (aTrgIdx2);
654         }
655         else
656         {
657           myOverlapTriangles1.Add (aTrgIdx1);
658           myOverlapTriangles2.Add (aTrgIdx2);
659         }
660 #endif
661       }
662       else if (aResult == BRepExtrema_ElementFilter::DoCheck)
663       {
664         if (!aPrism1.IsInited)
665         {
666           aPrism1.Init (aTrg1Vert1, aTrg1Vert2, aTrg1Vert3, theToler);
667         }
668
669         BVH_Vec3d aTrg2Vert1;
670         BVH_Vec3d aTrg2Vert2;
671         BVH_Vec3d aTrg2Vert3;
672
673         mySet2->GetVertices (aTrgIdx2,
674                              aTrg2Vert1,
675                              aTrg2Vert2,
676                              aTrg2Vert3);
677
678         BRepExtrema_BoundingPrism aPrism2 (aTrg2Vert1,
679                                            aTrg2Vert2,
680                                            aTrg2Vert3,
681                                            theToler);
682
683         if (prismsIntersected (aPrism1, aPrism2))
684         {
685           getSetOfFaces (myOverlapSubShapes1, aFaceIdx1).Add (aFaceIdx2);
686           getSetOfFaces (myOverlapSubShapes2, aFaceIdx2).Add (aFaceIdx1);
687         }
688       }
689     }
690   }
691 }
692
693 //=======================================================================
694 //function : Perform
695 //purpose  : Performs search for overlapped faces
696 //=======================================================================
697 void BRepExtrema_OverlapTool::Perform (const Standard_Real theTolerance)
698 {
699   if (mySet1.IsNull() || mySet2.IsNull())
700   {
701     return;
702   }
703
704   BRepExtrema_StackItem aStack[96];
705
706   const opencascade::handle<BVH_Tree<Standard_Real, 3> >& aBVH1 = mySet1->BVH();
707   const opencascade::handle<BVH_Tree<Standard_Real, 3> >& aBVH2 = mySet2->BVH();
708
709   if (aBVH1.IsNull() || aBVH2.IsNull())
710   {
711     return;
712   }
713
714   BRepExtrema_StackItem aNodes; // current pair of nodes
715
716   Standard_Integer aHead = -1; // stack head position
717
718   for (;;)
719   {
720     BVH_Vec4i aNodeData1 = aBVH1->NodeInfoBuffer()[aNodes.Node1];
721     BVH_Vec4i aNodeData2 = aBVH2->NodeInfoBuffer()[aNodes.Node2];
722
723     if (aNodeData1.x() != 0 && aNodeData2.x() != 0) // leaves
724     {
725       if (theTolerance == 0.0)
726       {
727         intersectTriangleRangesExact (aNodeData1, aNodeData2);
728       }
729       else
730       {
731         intersectTriangleRangesToler (aNodeData1, aNodeData2, theTolerance);
732       }
733
734       if (aHead < 0)
735         break;
736
737       aNodes = aStack[aHead--];
738     }
739     else
740     {
741       BRepExtrema_StackItem aPairsToProcess[4];
742
743       Standard_Integer aNbPairs = 0;
744
745       if (aNodeData1.x() == 0) // inner node
746       {
747         const BVH_Vec3d& aMinPntLft1 = aBVH1->MinPoint (aNodeData1.y());
748         const BVH_Vec3d& aMaxPntLft1 = aBVH1->MaxPoint (aNodeData1.y());
749         const BVH_Vec3d& aMinPntRgh1 = aBVH1->MinPoint (aNodeData1.z());
750         const BVH_Vec3d& aMaxPntRgh1 = aBVH1->MaxPoint (aNodeData1.z());
751
752         if (aNodeData2.x() == 0) // inner node
753         {
754           const BVH_Vec3d& aMinPntLft2 = aBVH2->MinPoint (aNodeData2.y());
755           const BVH_Vec3d& aMaxPntLft2 = aBVH2->MaxPoint (aNodeData2.y());
756           const BVH_Vec3d& aMinPntRgh2 = aBVH2->MinPoint (aNodeData2.z());
757           const BVH_Vec3d& aMaxPntRgh2 = aBVH2->MaxPoint (aNodeData2.z());
758
759           if (overlapBoxes (aMinPntLft1, aMaxPntLft1, aMinPntLft2, aMaxPntLft2, theTolerance))
760           {
761             aPairsToProcess[aNbPairs++] = BRepExtrema_StackItem (aNodeData1.y(), aNodeData2.y());
762           }
763           if (overlapBoxes (aMinPntLft1, aMaxPntLft1, aMinPntRgh2, aMaxPntRgh2, theTolerance))
764           {
765             aPairsToProcess[aNbPairs++] = BRepExtrema_StackItem (aNodeData1.y(), aNodeData2.z());
766           }
767           if (overlapBoxes (aMinPntRgh1, aMaxPntRgh1, aMinPntLft2, aMaxPntLft2, theTolerance))
768           {
769             aPairsToProcess[aNbPairs++] = BRepExtrema_StackItem (aNodeData1.z(), aNodeData2.y());
770           }
771           if (overlapBoxes (aMinPntRgh1, aMaxPntRgh1, aMinPntRgh2, aMaxPntRgh2, theTolerance))
772           {
773             aPairsToProcess[aNbPairs++] = BRepExtrema_StackItem (aNodeData1.z(), aNodeData2.z());
774           }
775         }
776         else
777         {
778           const BVH_Vec3d& aMinPntLeaf = aBVH2->MinPoint (aNodes.Node2);
779           const BVH_Vec3d& aMaxPntLeaf = aBVH2->MaxPoint (aNodes.Node2);
780
781           if (overlapBoxes (aMinPntLft1, aMaxPntLft1, aMinPntLeaf, aMaxPntLeaf, theTolerance))
782           {
783             aPairsToProcess[aNbPairs++] = BRepExtrema_StackItem (aNodeData1.y(), aNodes.Node2);
784           }
785           if (overlapBoxes (aMinPntRgh1, aMaxPntRgh1, aMinPntLeaf, aMaxPntLeaf, theTolerance))
786           {
787             aPairsToProcess[aNbPairs++] = BRepExtrema_StackItem (aNodeData1.z(), aNodes.Node2);
788           }
789         }
790       }
791       else
792       {
793         const BVH_Vec3d& aMinPntLeaf = aBVH1->MinPoint (aNodes.Node1);
794         const BVH_Vec3d& aMaxPntLeaf = aBVH1->MaxPoint (aNodes.Node1);
795
796         const BVH_Vec3d& aMinPntLft2 = aBVH2->MinPoint (aNodeData2.y());
797         const BVH_Vec3d& aMaxPntLft2 = aBVH2->MaxPoint (aNodeData2.y());
798         const BVH_Vec3d& aMinPntRgh2 = aBVH2->MinPoint (aNodeData2.z());
799         const BVH_Vec3d& aMaxPntRgh2 = aBVH2->MaxPoint (aNodeData2.z());
800
801         if (overlapBoxes (aMinPntLft2, aMaxPntLft2, aMinPntLeaf, aMaxPntLeaf, theTolerance))
802         {
803           aPairsToProcess[aNbPairs++] = BRepExtrema_StackItem (aNodes.Node1, aNodeData2.y());
804         }
805         if (overlapBoxes (aMinPntRgh2, aMaxPntRgh2, aMinPntLeaf, aMaxPntLeaf, theTolerance))
806         {
807           aPairsToProcess[aNbPairs++] = BRepExtrema_StackItem (aNodes.Node1, aNodeData2.z());
808         }
809       }
810
811       if (aNbPairs > 0)
812       {
813         aNodes = aPairsToProcess[0];
814
815         for (Standard_Integer anIdx = 1; anIdx < aNbPairs; ++anIdx)
816         {
817           aStack[++aHead] = aPairsToProcess[anIdx];
818         }
819       }
820       else
821       {
822         if (aHead < 0)
823           break;
824
825         aNodes = aStack[aHead--];
826       }
827     }
828   }
829
830   myIsDone = Standard_True;
831 }