0027117: BRepClass3d_SolidClassifier doesn't take into account vertex/edge/face toler...
[occt.git] / src / BRepClass3d / BRepClass3d_SClassifier.cxx
1 // Created on: 1996-07-15
2 // Created by: Laurent BUCHARD
3 // Copyright (c) 1996-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 //  Modified by skv - Thu Sep  4 11:22:05 2003 OCC578
18
19 #include <BRep_Tool.hxx>
20 #include <BRepClass3d_Intersector3d.hxx>
21 #include <BRepClass3d_SClassifier.hxx>
22 #include <BRepClass3d_SolidExplorer.hxx>
23 #include <BRepTopAdaptor_FClass2d.hxx>
24 #include <ElCLib.hxx>
25 #include <Geom_Surface.hxx>
26 #include <gp_Lin.hxx>
27 #include <gp_Pnt.hxx>
28 #include <gp_Vec.hxx>
29 #include <IntCurvesFace_Intersector.hxx>
30 #include <math_BullardGenerator.hxx>
31 #include <Standard_DomainError.hxx>
32 #include <TopoDS.hxx>
33 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
34 #include <TopTools_MapOfShape.hxx>
35 #include <TopExp.hxx>
36
37 #include <vector>
38
39 // modified by NIZHNY-MKK  Mon Jun 21 15:13:40 2004
40 static
41   Standard_Boolean FaceNormal (const TopoDS_Face& aF,
42                                const Standard_Real U,
43                                const Standard_Real V,
44                                gp_Dir& aDN);
45
46 static 
47   Standard_Real GetAddToParam(const gp_Lin& L,const Standard_Real P,const Bnd_Box& B);
48
49 //gets transition of line <L> passing through/near the edge <e> of faces <f1>, <f2>. <param> is
50 // a parameter on the edge where the minimum distance between <l> and <e> was found 
51 static Standard_Integer GetTransi(const TopoDS_Face& f1, const TopoDS_Face& f2, const TopoDS_Edge e, 
52                      Standard_Real param, const Geom_Line& L, IntCurveSurface_TransitionOnCurve& trans);
53
54 static Standard_Boolean GetNormalOnFaceBound(const TopoDS_Edge& E, const TopoDS_Face& F, Standard_Real param, gp_Dir& OutDir);
55
56 static void Trans(Standard_Real parmin, IntCurveSurface_TransitionOnCurve& tran, int& state);
57
58 //=======================================================================
59 //function : BRepClass3d_SClassifier
60 //purpose  : 
61 //=======================================================================
62 BRepClass3d_SClassifier::BRepClass3d_SClassifier() 
63
64 }
65
66
67 //=======================================================================
68 //function : BRepClass3d_SClassifier
69 //purpose  : 
70 //=======================================================================
71 BRepClass3d_SClassifier::BRepClass3d_SClassifier(BRepClass3d_SolidExplorer& S,
72                                                  const gp_Pnt&  P,
73                                                  const Standard_Real Tol) { 
74   if(S.Reject(P)) { 
75     myState=3; //-- in ds solid case without face 
76   }
77   else { 
78     Perform(S,P,Tol);
79   }
80 }
81
82
83 //=======================================================================
84 //function : PerformInfinitePoint
85 //purpose  : 
86 //=======================================================================
87 void BRepClass3d_SClassifier::PerformInfinitePoint(BRepClass3d_SolidExplorer& aSE,
88                                                    const Standard_Real /*Tol*/) {
89
90   //Take a normal to the first extracted face in its random inner point
91   //and intersect this reversed normal with the faces of the solid.
92   //If the min.par.-intersection point is
93   // a) inner point of a face
94   // b) transition is not TANGENT
95   //    (the line does not touch the face but pierces it)
96   //then set <myState> to IN or OUT according to transition
97   //else take the next random point inside the min.par.-intersected face
98   //and continue
99   
100   if(aSE.Reject(gp_Pnt(0,0,0))) { 
101     myState=3; //-- in ds solid case without face 
102     return;
103   }
104   //
105   //------------------------------------------------------------
106   // 1
107   Standard_Boolean bFound;
108   Standard_Real aParam, aU = 0., aV = 0.;
109   gp_Pnt aPoint;
110   gp_Dir aDN;
111
112   math_BullardGenerator aRandomGenerator;
113   myFace.Nullify();
114   myState=2;
115
116   // Collect faces in sequence to iterate
117   std::vector<TopoDS_Face> aFaces;
118   for (aSE.InitShell(); aSE.MoreShell(); aSE.NextShell())
119   {
120     for (aSE.InitFace(); aSE.MoreFace(); aSE.NextFace())
121     {
122       aFaces.push_back (aSE.CurrentFace());
123     }
124   }
125
126   // iteratively try up to 10 probing points from each face
127   const Standard_Integer NB_MAX_POINTS_PER_FACE = 10;
128   for (Standard_Integer itry = 0; itry < NB_MAX_POINTS_PER_FACE; itry++)
129   {
130     for (std::vector<TopoDS_Face>::iterator iFace = aFaces.begin(); iFace != aFaces.end(); ++iFace)
131     {
132       TopoDS_Face aF = *iFace;
133
134       TopAbs_State aState = TopAbs_OUT;
135       IntCurveSurface_TransitionOnCurve aTransition = IntCurveSurface_Tangent;
136
137       aParam = 0.1 + 0.8 * aRandomGenerator.NextReal(); // random number in range [0.1, 0.9]
138       bFound = aSE.FindAPointInTheFace(aF, aPoint, aU, aV, aParam);
139       if (!bFound || !FaceNormal(aF, aU, aV, aDN))
140         continue;
141
142       gp_Lin aLin(aPoint, -aDN);
143       Standard_Real parmin = RealLast();
144       for (aSE.InitShell();aSE.MoreShell();aSE.NextShell()) { 
145         if (aSE.RejectShell(aLin) == Standard_False) { 
146           for (aSE.InitFace();aSE.MoreFace(); aSE.NextFace()) {
147             if (aSE.RejectFace(aLin) == Standard_False) { 
148               TopoDS_Shape aLocalShape = aSE.CurrentFace();
149               TopoDS_Face CurFace = TopoDS::Face(aLocalShape);
150               IntCurvesFace_Intersector& Intersector3d = aSE.Intersector(CurFace);
151               Intersector3d.Perform(aLin,-RealLast(),parmin); 
152
153               if(Intersector3d.IsDone()) {
154                 if(Intersector3d.NbPnt()) { 
155                   Standard_Integer imin = 1;
156                   for (Standard_Integer i = 2; i <= Intersector3d.NbPnt(); i++)
157                     if (Intersector3d.WParameter(i) < Intersector3d.WParameter(imin))
158                       imin = i;
159                   parmin = Intersector3d.WParameter(imin);
160                   aState = Intersector3d.State(imin);
161                   aTransition = Intersector3d.Transition(imin);
162                 }
163               }
164             }
165           }
166         }
167         else
168           myState = 1;
169       } //end of loop on the whole solid
170         
171       if (aState == TopAbs_IN)
172       {
173         if (aTransition == IntCurveSurface_Out) { 
174           //-- The line is going from inside the solid to outside 
175           //-- the solid.
176           myState = 3; //-- IN --
177           return;
178         }
179         else if (aTransition == IntCurveSurface_In) { 
180           myState = 4; //-- OUT --
181           return;
182         }
183       }
184     } // iteration by faces
185   } // iteration by points
186 }
187
188 //=======================================================================
189 //function : Perform
190 //purpose  : 
191 //=======================================================================
192 void BRepClass3d_SClassifier::Perform(BRepClass3d_SolidExplorer& SolidExplorer,
193                                       const gp_Pnt&  P,
194                                       const Standard_Real Tol)
195
196   if(SolidExplorer.Reject(P))
197   {
198     // Solid without faces => the whole space. Always in.
199     myState = 3; // IN
200     return;
201   }
202
203   const BRepClass3d_BndBoxTree & aTree = SolidExplorer.GetTree();
204   const TopTools_IndexedMapOfShape & aMapEV = SolidExplorer.GetMapEV();
205
206   // Vertices/Edges vs Point.
207   BRepClass3d_BndBoxTreeSelectorPoint aSelectorPoint(aMapEV);
208   aSelectorPoint.SetCurrentPoint(P);
209
210   Standard_Integer SelsVE = 0;
211   SelsVE = aTree.Select(aSelectorPoint); 
212
213   if (SelsVE > 0)
214   {
215     // The point P lays inside the tolerance area of vertices/edges => return ON state.
216     myState = 2; // ON.
217     return;
218   }
219
220   TopTools_IndexedDataMapOfShapeListOfShape mapEF;
221   TopExp::MapShapesAndAncestors(SolidExplorer.GetShape(), TopAbs_EDGE, TopAbs_FACE, mapEF);
222   
223   BRepClass3d_BndBoxTreeSelectorLine aSelectorLine(aMapEV);
224
225   myFace.Nullify();
226   myState = 0;
227
228   gp_Lin L;
229   Standard_Real Par;
230
231   // We compute the intersection between the line built in the Solid Explorer and the shape.
232   //-- --------------------------------------------------------------------------------
233   // Calculate intersection with the face closest to the direction of bounding boxes 
234   // by priority so that to have the smallest possible parmin.
235   // Optimization to produce as much as possible rejections with other faces.
236   Standard_Integer iFlag;
237
238   // If found line passes through a bound of any face, it means that the line
239   // is not found properly and it is necessary to repeat whole procedure.
240   // That's why the main loop while is added.
241   Standard_Boolean isFaultyLine = Standard_True;
242   Standard_Integer anIndFace    = 0;
243   Standard_Real    parmin = 0.0;
244   while (isFaultyLine)
245   {
246     if (anIndFace == 0)
247       iFlag = SolidExplorer.Segment(P,L,Par);
248     else
249       iFlag = SolidExplorer.OtherSegment(P,L,Par);
250
251     Standard_Integer aCurInd = SolidExplorer.GetFaceSegmentIndex();
252
253     if (aCurInd > anIndFace)
254       anIndFace = aCurInd;
255     else
256     {
257       myState = 1; // Faulty.
258       return;
259     }
260
261     if (iFlag==1)
262     {
263       // IsOnFace iFlag==1 i.e face is Infinite
264       myState = 2; // ON.
265       return;
266     }
267     if (iFlag == 2) 
268     {
269       myState = 4; // OUT.
270       return;
271     }
272
273     // Check if the point is ON surface but OUT of the face.
274     // Just skip this face because it is bad for classification.
275     if (iFlag == 3)
276       continue;
277
278     isFaultyLine = Standard_False;
279     parmin = RealLast();
280
281     Standard_Real NearFaultPar = RealLast(); // Parameter on line.
282     aSelectorLine.ClearResults();
283     aSelectorLine.SetCurrentLine(L, Par);
284     Standard_Integer SelsEVL = 0;
285     SelsEVL = aTree.Select(aSelectorLine); //SelsEE > 0 => Line/Edges & Line/Vertex intersection
286     if (SelsEVL > 0 )
287     {    
288       // Line and edges / vertices interference.
289       Standard_Integer VLInterNb = aSelectorLine.GetNbVertParam();
290       TopoDS_Vertex NearIntVert;
291       TopTools_MapOfShape LVInts;
292       for (Standard_Integer i = 1; i <= VLInterNb; i++)
293       {
294         // Line and vertex.
295         Standard_Real LP = 0;
296         TopoDS_Vertex V;
297         aSelectorLine.GetVertParam(i, V, LP);
298
299         LVInts.Add(V);
300         if (Abs(LP) < Abs(NearFaultPar))
301           NearFaultPar = LP;
302       }
303
304       Standard_Real param = 0.0;
305       TopoDS_Edge EE;
306       Standard_Real Lpar = RealLast();
307       for (Standard_Integer i = 1; i <= aSelectorLine.GetNbEdgeParam(); i++)
308       {
309         // Line and edge.
310         aSelectorLine.GetEdgeParam(i, EE, param, Lpar);
311         const TopTools_ListOfShape& ffs = mapEF.FindFromKey(EE); //ffs size == 2
312         if (ffs.Extent() != 2)
313           continue;
314         TopoDS_Face f1 = TopoDS::Face(ffs.First());
315         TopoDS_Face f2 = TopoDS::Face(ffs.Last());
316         IntCurveSurface_TransitionOnCurve tran;
317         TopoDS_Vertex V1, V2;
318         TopExp::Vertices(EE, V1, V2);
319         if (LVInts.Contains(V1) || LVInts.Contains(V2))
320           continue;
321         Standard_Integer Tst = GetTransi(f1, f2, EE, param, L, tran);
322         if (Tst == 1 && Abs(Lpar) < Abs(parmin))
323         {
324           parmin = Lpar;
325           Trans(parmin, tran, myState);
326         }
327         else if (Abs(Lpar) < Abs(NearFaultPar))
328           NearFaultPar = Lpar;
329       }
330     }
331
332     for(SolidExplorer.InitShell();
333         SolidExplorer.MoreShell() && !isFaultyLine;
334         SolidExplorer.NextShell())
335     {
336         if(SolidExplorer.RejectShell(L) == Standard_False)
337         {
338           for(SolidExplorer.InitFace(); 
339               SolidExplorer.MoreFace() && !isFaultyLine; 
340               SolidExplorer.NextFace())
341           {
342               if(SolidExplorer.RejectFace(L) == Standard_False)
343               {
344                 TopoDS_Shape aLocalShape = SolidExplorer.CurrentFace();
345                 TopoDS_Face f = TopoDS::Face(aLocalShape);
346                 IntCurvesFace_Intersector& Intersector3d = SolidExplorer.Intersector(f);
347
348                 // Prolong segment, since there are cases when
349                 // the intersector does not find intersection points with the original
350                 // segment due to rough triangulation of a parameterized surface.
351                 Standard_Real addW = Max(10*Tol, 0.01*Par);
352                 Standard_Real AddW = addW;
353
354                 Bnd_Box aBoxF = Intersector3d.Bounding();
355
356                 // The box must be finite in order to correctly prolong the segment to its bounds.
357                 if (!aBoxF.IsVoid() && !aBoxF.IsWhole())
358                 {
359                   Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
360                   aBoxF.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
361
362                   Standard_Real boxaddW = GetAddToParam(L,Par,aBoxF);
363                   addW = Max(addW,boxaddW);
364                 }
365
366                 Standard_Real minW = -AddW;
367                 Standard_Real maxW = Min(Par*10,Par+addW);
368                 Intersector3d.Perform(L,minW,maxW);
369                 if(Intersector3d.IsDone())
370                 {
371                   for (Standard_Integer i = 1; i <= Intersector3d.NbPnt(); i++)
372                   {
373                     if(Abs(Intersector3d.WParameter(i)) < Abs(parmin) - Precision::PConfusion())
374                     {
375                       parmin = Intersector3d.WParameter(i);
376                       TopAbs_State aState = Intersector3d.State(i);
377                       if(Abs(parmin)<=Tol) 
378                       { 
379                         myState = 2;
380                         myFace  = f;
381                       }
382                       //  Treatment of case TopAbs_ON separately.
383                       else if(aState==TopAbs_IN)
384                       {
385                         //-- The intersection point between the line and a face F 
386                         // -- of the solid is in the face F 
387
388                         IntCurveSurface_TransitionOnCurve tran = Intersector3d.Transition(i);
389                         if (tran == IntCurveSurface_Tangent)
390                         {
391                           #ifdef OCCT_DEBUG
392                             cout<<"*Problem ds BRepClass3d_SClassifier.cxx"<<endl;
393                           #endif
394                           continue; // ignore this point
395                         }
396
397                         Trans(parmin, tran, myState);
398                         myFace  = f;
399                       }
400                       // If the state is TopAbs_ON, it is necessary to chose
401                       // another line and to repeat the whole procedure.
402                       else if(aState==TopAbs_ON)
403                       {
404                         isFaultyLine = Standard_True;
405                         break;
406                       }
407                     }
408                     else
409                     {
410                       //-- No point has been found by the Intersector3d.
411                       //-- Or a Point has been found with a greater parameter.
412                     }
413                   } //-- loop by intersection points
414                 } //-- Face has not been rejected
415                 else
416                   myState = 1;
417               }
418           } //-- Exploration of the faces
419         } //-- Shell has not been rejected
420         else
421           myState = 1;
422     } //-- Exploration of the shells
423
424     if (NearFaultPar != RealLast() &&
425         Abs(parmin) >= Abs(NearFaultPar) - Precision::PConfusion())
426     {
427       isFaultyLine = Standard_True;
428     }
429   }
430
431 #ifdef OCCT_DEBUG
432   //#################################################
433   SolidExplorer.DumpSegment(P,L,parmin,State());
434   //#################################################
435 #endif
436 }
437
438
439 TopAbs_State BRepClass3d_SClassifier::State() const
440 {
441   if(myState == 2)
442     return(TopAbs_ON);
443   else if(myState == 3)
444     return(TopAbs_IN);
445   else if(myState == 4)
446     return(TopAbs_OUT);
447
448   // return OUT state when there is an error during execution.
449   return(TopAbs_OUT);
450 }
451
452 TopoDS_Face BRepClass3d_SClassifier::Face() const {  
453   return(myFace);
454 }
455
456 Standard_Boolean BRepClass3d_SClassifier::Rejected() const { 
457   return(myState==1); 
458 }
459
460   
461 Standard_Boolean BRepClass3d_SClassifier::IsOnAFace() const { 
462   return(myState==2);
463 }
464
465
466 void BRepClass3d_SClassifier::ForceIn() {
467   myState=3;
468 }
469
470 void BRepClass3d_SClassifier::ForceOut() { 
471   myState=4;
472 }
473
474 Standard_Real GetAddToParam(const gp_Lin&       L,
475                             const Standard_Real P,
476                             const Bnd_Box&      B)
477 {
478   Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
479   B.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
480   Standard_Real x[2] = {aXmin,aXmax}, y[2] = {aYmin,aYmax}, z[2] = {aZmin,aZmax};
481   Standard_Integer i = 0, j = 0, k = 0;
482   Standard_Real Par = P;
483   for(i = 0 ; i < 2; i++) {
484     for(j = 0; j < 2; j++) {
485       for(k = 0; k < 2; k++) {
486         Standard_Real X = fabs(x[i]-L.Location().X());
487         Standard_Real Y = fabs(y[j]-L.Location().Y());
488         Standard_Real Z = fabs(z[k]-L.Location().Z());
489         if(X < 1.e+20 && Y < 1.e+20 && Z < 1.e+20) {
490           gp_Pnt aP(x[i],y[j],z[k]);
491           Standard_Real par = ElCLib::Parameter(L,aP);
492           if(par > Par)
493             Par = par;
494         }
495         else
496           return 1.e+20;
497       }
498     }
499   }
500   return Par - P;
501 }
502 //=======================================================================
503 //function : FaceNormal
504 //purpose  : 
505 //=======================================================================
506 Standard_Boolean FaceNormal (const TopoDS_Face& aF,
507                              const Standard_Real U,
508                              const Standard_Real V,
509                              gp_Dir& aDN)
510 {
511   gp_Pnt aPnt ;
512   gp_Vec aD1U, aD1V, aN;
513   Handle(Geom_Surface) aS;
514
515   aS=BRep_Tool::Surface(aF);
516   aS->D1 (U, V, aPnt, aD1U, aD1V);
517   aN=aD1U.Crossed(aD1V);
518   if (aN.Magnitude() <= gp::Resolution())
519     return Standard_False;
520   
521   aN.Normalize();  
522   aDN.SetXYZ(aN.XYZ());
523   if (aF.Orientation() == TopAbs_REVERSED){
524     aDN.Reverse();
525   }
526   return Standard_True;
527 }
528
529 //=======================================================================
530 //function : GetNormalOnFaceBound
531 //purpose  : 
532 //=======================================================================
533 static Standard_Boolean GetNormalOnFaceBound(const TopoDS_Edge& E,
534                                              const TopoDS_Face& F,
535                                              const Standard_Real param,
536                                              gp_Dir& OutDir)
537
538   Standard_Real f = 0, l = 0;
539
540   gp_Pnt2d P2d;
541   Handle(Geom2d_Curve) c2d = BRep_Tool::CurveOnSurface(E, F, f, l);
542   if (c2d.IsNull())
543     return Standard_False;
544   if (param < f || param > l)
545     return Standard_False;
546   c2d->D0(param, P2d);
547   if (!FaceNormal(F, P2d.X(), P2d.Y(), OutDir))
548     return Standard_False;
549   return Standard_True;
550 }
551
552 //=======================================================================
553 //function : GetTransi
554 //purpose  : 
555 //=======================================================================
556 static Standard_Integer GetTransi(const TopoDS_Face& f1,
557                                   const TopoDS_Face& f2,
558                                   const TopoDS_Edge e,
559                                   const Standard_Real param,
560                                   const Geom_Line& L,
561                                   IntCurveSurface_TransitionOnCurve& trans)
562 {
563   //return statuses:
564   //1 => OK
565   //0 => skip
566   //-1 => probably a faulty line
567   gp_Dir nf1, nf2;
568   if (!GetNormalOnFaceBound(e, f1, param, nf1))
569     return -1;
570   if (!GetNormalOnFaceBound(e, f2, param, nf2))
571     return -1;
572
573   const gp_Dir& LDir = L.Lin().Direction();
574
575   if(Abs(LDir.Dot(nf1)) < Precision::Angular() || Abs(LDir.Dot(nf2)) < Precision::Angular())
576   {
577     //line is orthogonal to normal(s)
578     //trans = IntCurveSurface_Tangent;
579     return -1;
580   }
581
582   if (nf1.IsEqual(nf2, Precision::Angular()))
583   {
584     Standard_Real angD = nf1.Dot(LDir);
585     if (Abs(angD) < Precision::Angular())
586       return -1;
587     else if (angD > 0)
588       trans = IntCurveSurface_Out;
589     else //angD < -Precision::Angular())
590       trans = IntCurveSurface_In;
591     return 1;
592   }
593
594   gp_Vec N = nf1^nf2;
595   gp_Dir ProjL = N.XYZ() ^ LDir.XYZ() ^ N.XYZ(); //proj LDir on the plane defined by nf1/nf2 directions
596
597   Standard_Real fAD = nf1.Dot(ProjL); 
598   Standard_Real sAD = nf2.Dot(ProjL);  
599
600   if (fAD < -Precision::Angular() && sAD < -Precision::Angular())
601     trans = IntCurveSurface_In;
602   else if (fAD > Precision::Angular() && sAD > Precision::Angular())
603     trans = IntCurveSurface_Out;
604   else
605     return 0;
606   return 1;
607 }
608
609 //=======================================================================
610 //function : Trans
611 //purpose  : 
612 //=======================================================================
613 static void Trans(const Standard_Real parmin, 
614                   IntCurveSurface_TransitionOnCurve& tran,
615                   int& state)
616 {
617   // if parmin is negative we should reverse transition
618   if (parmin < 0)
619     tran = (tran == IntCurveSurface_Out ? IntCurveSurface_In : IntCurveSurface_Out);
620
621   if(tran == IntCurveSurface_Out)
622     //-- The line is going from inside the solid to outside 
623     //-- the solid.
624     state = 3; // IN
625   else
626     state = 4; // OUT
627 }