0030790: [REGRESSION] Modeling Algorithms - Crash when loading specific step file
[occt.git] / src / IntCurvesFace / IntCurvesFace_Intersector.cxx
1 // Created on: 1996-06-03
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 #define OPTIMISATION 1 
18
19 #include <Adaptor3d_HCurve.hxx>
20 #include <Adaptor3d_HSurfaceTool.hxx>
21 #include <Bnd_BoundSortBox.hxx>
22 #include <Bnd_Box.hxx>
23 #include <BRepAdaptor_HSurface.hxx>
24 #include <BRepClass_FaceClassifier.hxx>
25 #include <BRepTopAdaptor_TopolTool.hxx>
26 #include <Geom_Line.hxx>
27 #include <GeomAdaptor_Curve.hxx>
28 #include <GeomAdaptor_HCurve.hxx>
29 #include <gp_Lin.hxx>
30 #include <gp_Pnt.hxx>
31 #include <gp_Pnt2d.hxx>
32 #include <IntCurvesFace_Intersector.hxx>
33 #include <IntCurveSurface_HInter.hxx>
34 #include <IntCurveSurface_IntersectionPoint.hxx>
35 #include <IntCurveSurface_SequenceOfPnt.hxx>
36 #include <IntCurveSurface_TheHCurveTool.hxx>
37 #include <IntCurveSurface_ThePolygonOfHInter.hxx>
38 #include <IntCurveSurface_ThePolyhedronOfHInter.hxx>
39 #include <IntCurveSurface_ThePolyhedronToolOfHInter.hxx>
40 #include <Intf_Tool.hxx>
41 #include <TopAbs.hxx>
42 #include <TopoDS_Face.hxx>
43 #include <BRep_Tool.hxx>
44 #include <TopoDS.hxx>
45 #include <Geom2dAPI_ProjectPointOnCurve.hxx>
46 //
47 static void ComputeSamplePars(const Handle(Adaptor3d_HSurface)& Hsurface, 
48                               const Standard_Integer nbsu,
49                               const Standard_Integer nbsv,
50                               Handle(TColStd_HArray1OfReal)& UPars, 
51                               Handle(TColStd_HArray1OfReal)& VPars)
52 {
53   Standard_Integer NbUInts = Hsurface->NbUIntervals(GeomAbs_C2);
54   Standard_Integer NbVInts = Hsurface->NbVIntervals(GeomAbs_C2);
55   TColStd_Array1OfReal UInts(1, NbUInts + 1);
56   TColStd_Array1OfReal VInts(1, NbVInts + 1);
57   Hsurface->UIntervals(UInts, GeomAbs_C2);
58   Hsurface->VIntervals(VInts, GeomAbs_C2);
59   //
60   TColStd_Array1OfInteger NbUSubInts(1, NbUInts);
61   TColStd_Array1OfInteger NbVSubInts(1, NbVInts);
62   //
63   Standard_Integer i, j, ind, NbU, NbV;
64   Standard_Real t, dt;
65   t = UInts(NbUInts + 1) - UInts(1);
66   t = 1. / t;
67   NbU = 0;
68   for(i = 1; i <= NbUInts; ++i)
69   {
70     dt = (UInts(i+1) - UInts(i));
71     NbUSubInts(i) = RealToInt(nbsu * dt * t) + 1;
72     NbU += NbUSubInts(i);
73   }
74   t = VInts(NbVInts + 1) - VInts(1);
75   t = 1. / t;
76   NbV = 0;
77   for(i = 1; i <= NbVInts; ++i)
78   {
79     dt = (VInts(i+1) - VInts(i));
80     NbVSubInts(i) = RealToInt(nbsv * dt * t) + 1;
81     NbV += NbVSubInts(i);
82   }
83   UPars = new TColStd_HArray1OfReal(1, NbU + 1);
84   VPars = new TColStd_HArray1OfReal(1, NbV + 1);
85   //
86   ind = 1;
87   for(i = 1; i <= NbUInts; ++i)
88   {
89     UPars->SetValue(ind++, UInts(i));
90     dt = (UInts(i+1) - UInts(i)) / NbUSubInts(i);
91     t = UInts(i);
92     for(j = 1; j < NbUSubInts(i); ++j)
93     {
94       t += dt;
95       UPars->SetValue(ind++, t);
96     }
97   }
98   UPars->SetValue(ind, UInts(NbUInts + 1));
99   //
100   ind = 1;
101   for(i = 1; i <= NbVInts; ++i)
102   {
103     VPars->SetValue(ind++, VInts(i));
104     dt = (VInts(i+1) - VInts(i)) / NbVSubInts(i);
105     t = VInts(i);
106     for(j = 1; j < NbVSubInts(i); ++j)
107     {
108       t += dt;
109       VPars->SetValue(ind++, t);
110     }
111   }
112   VPars->SetValue(ind, VInts(NbVInts + 1));
113 }
114 //
115 //=======================================================================
116 //function : SurfaceType
117 //purpose  : 
118 //=======================================================================
119 GeomAbs_SurfaceType IntCurvesFace_Intersector::SurfaceType() const 
120
121   return(Adaptor3d_HSurfaceTool::GetType(Hsurface));
122 }
123 //=======================================================================
124 //function : IntCurvesFace_Intersector
125 //purpose  : 
126 //=======================================================================
127 IntCurvesFace_Intersector::IntCurvesFace_Intersector(const TopoDS_Face& Face,
128                                                      const Standard_Real aTol,
129                                                      const Standard_Boolean aRestr,
130                                                      const Standard_Boolean UseBToler)
131                                                     
132
133   Tol(aTol),
134   done(Standard_False),
135   myReady(Standard_False),
136   nbpnt(0),
137   PtrOnPolyhedron(NULL),
138   PtrOnBndBounding(NULL),
139   myUseBoundTol (UseBToler),
140   myIsParallel(Standard_False)
141
142   BRepAdaptor_Surface surface;
143   face = Face;
144   surface.Initialize(Face, aRestr);
145   Hsurface = new BRepAdaptor_HSurface(surface);
146   myTopolTool = new BRepTopAdaptor_TopolTool(Hsurface);
147   
148   GeomAbs_SurfaceType SurfaceType = Adaptor3d_HSurfaceTool::GetType(Hsurface);
149   if(   (SurfaceType != GeomAbs_Plane) 
150      && (SurfaceType != GeomAbs_Cylinder)
151      && (SurfaceType != GeomAbs_Cone) 
152      && (SurfaceType != GeomAbs_Sphere)
153      && (SurfaceType != GeomAbs_Torus)) {
154     Standard_Integer nbsu,nbsv;
155     Standard_Real U0,V0,U1,V1;
156     U0 = Hsurface->FirstUParameter();
157     U1 = Hsurface->LastUParameter();
158     V0 = Hsurface->FirstVParameter();
159     V1 = Hsurface->LastVParameter();
160     //
161     nbsu = myTopolTool->NbSamplesU();
162     nbsv = myTopolTool->NbSamplesV();
163     //
164     Standard_Real aURes = Hsurface->UResolution(1.0);
165     Standard_Real aVRes = Hsurface->VResolution(1.0);
166
167     // Checking correlation between number of samples and length of the face along each axis
168     const Standard_Real aTresh = 100.0;
169     Standard_Integer aMinSamples = 20;
170     const Standard_Integer aMaxSamples = 40;
171     const Standard_Integer aMaxSamples2 = aMaxSamples * aMaxSamples;
172     Standard_Real dU = (U1 - U0) / aURes;
173     Standard_Real dV = (V1 - V0) / aVRes;
174     if (nbsu < aMinSamples) nbsu = aMinSamples;
175     if (nbsv < aMinSamples) nbsv = aMinSamples;
176     if (nbsu > aMaxSamples) nbsu = aMaxSamples;
177     if (nbsv > aMaxSamples) nbsv = aMaxSamples;
178
179     if (dU > Precision::Confusion() && dV > Precision::Confusion()) {
180       if (Max(dU, dV) > Min(dU, dV) * aTresh)
181       {
182         aMinSamples = 10;
183         nbsu = (Standard_Integer)(Sqrt(dU / dV) * aMaxSamples);
184         if (nbsu < aMinSamples) nbsu = aMinSamples;
185         nbsv = aMaxSamples2 / nbsu;
186         if (nbsv < aMinSamples)
187         {
188           nbsv = aMinSamples;
189           nbsu = aMaxSamples2 / aMinSamples;
190         }
191       }
192     }
193     else {
194       return; // surface has no extension along one of directions
195     }
196
197     Standard_Integer NbUOnS = Hsurface->NbUIntervals(GeomAbs_C2);
198     Standard_Integer NbVOnS = Hsurface->NbVIntervals(GeomAbs_C2);
199
200     if(NbUOnS > 1 || NbVOnS > 1)
201     {
202       Handle(TColStd_HArray1OfReal) UPars, VPars;
203       ComputeSamplePars(Hsurface, nbsu, nbsv, UPars, VPars);
204       PtrOnPolyhedron = (IntCurveSurface_ThePolyhedronOfHInter *)
205         new IntCurveSurface_ThePolyhedronOfHInter(Hsurface, UPars->ChangeArray1(), 
206                                                             VPars->ChangeArray1());
207     }
208     else 
209     {
210       PtrOnPolyhedron = (IntCurveSurface_ThePolyhedronOfHInter *)
211         new IntCurveSurface_ThePolyhedronOfHInter(Hsurface,nbsu,nbsv,U0,V0,U1,V1);
212     }
213   }
214   myReady = Standard_True;
215 }
216 //=======================================================================
217 //function : InternalCall
218 //purpose  : 
219 //=======================================================================
220 void IntCurvesFace_Intersector::InternalCall(const IntCurveSurface_HInter &HICS,
221                                              const Standard_Real parinf,
222                                              const Standard_Real parsup) 
223 {
224   if(HICS.IsDone() && HICS.NbPoints() > 0) {
225     //Calculate tolerance for 2d classifier
226     Standard_Real mintol3d = BRep_Tool::Tolerance(face);
227     Standard_Real maxtol3d = mintol3d;
228     Standard_Real mintol2d = Tol, maxtol2d = Tol;
229     TopExp_Explorer anExp(face, TopAbs_EDGE);
230     for(; anExp.More(); anExp.Next())
231     {
232       Standard_Real curtol = BRep_Tool::Tolerance(TopoDS::Edge(anExp.Current()));
233       mintol3d = Min(mintol3d, curtol);
234       maxtol3d = Max(maxtol3d, curtol);
235     }
236     Standard_Real minres = Max(Hsurface->UResolution(mintol3d), Hsurface->VResolution(mintol3d));
237     Standard_Real maxres = Max(Hsurface->UResolution(maxtol3d), Hsurface->VResolution(maxtol3d));
238     mintol2d = Max(minres, Tol); 
239     maxtol2d = Max(maxres, Tol);
240     //
241     Handle(BRepTopAdaptor_TopolTool) anAdditionalTool;
242     for(Standard_Integer index=HICS.NbPoints(); index>=1; index--) {  
243       const IntCurveSurface_IntersectionPoint& HICSPointindex = HICS.Point(index);
244       gp_Pnt2d Puv(HICSPointindex.U(),HICSPointindex.V());
245
246       //TopAbs_State currentstate = myTopolTool->Classify(Puv,Tol);
247       TopAbs_State currentstate = myTopolTool->Classify(Puv, !myUseBoundTol ? 0 : mintol2d);
248       if(myUseBoundTol && currentstate == TopAbs_OUT && maxtol2d > mintol2d) {
249         if(anAdditionalTool.IsNull())
250         {
251           anAdditionalTool = new BRepTopAdaptor_TopolTool(Hsurface);
252         }
253         currentstate = anAdditionalTool->Classify(Puv,maxtol2d);
254         if(currentstate == TopAbs_ON)
255         {
256           currentstate = TopAbs_OUT;
257           //Find out nearest edge and it's tolerance
258           anExp.Init(face, TopAbs_EDGE);
259           for(; anExp.More(); anExp.Next())
260           {
261             TopoDS_Edge anE = TopoDS::Edge(anExp.Current());
262             Standard_Real curtol = BRep_Tool::Tolerance(anE);
263             Standard_Real tol2d = Max(Hsurface->UResolution(curtol), Hsurface->VResolution(curtol));
264             tol2d = Max(tol2d, Tol);
265             Standard_Real f, l;
266             Handle(Geom2d_Curve) aPC = BRep_Tool::CurveOnSurface(anE, face, f, l);
267             Geom2dAPI_ProjectPointOnCurve aProj(Puv, aPC, f, l);
268             if(aProj.NbPoints() > 0)
269             {
270               Standard_Real d = aProj.LowerDistance();
271               if(d <= tol2d)
272               {
273                 //Nearest edge is found, state is really ON
274                 currentstate = TopAbs_ON;
275                 break;
276               }
277             }
278           }
279         }
280       }
281       if(currentstate==TopAbs_IN || currentstate==TopAbs_ON) { 
282         Standard_Real HICSW = HICSPointindex.W();
283         if(HICSW >= parinf && HICSW <= parsup ) { 
284           Standard_Real U          = HICSPointindex.U();
285           Standard_Real V          = HICSPointindex.V();
286           Standard_Real W          = HICSW; 
287           IntCurveSurface_TransitionOnCurve transition = HICSPointindex.Transition();
288           gp_Pnt pnt        = HICSPointindex.Pnt();
289           //      state      = currentstate;
290           //  Modified by skv - Wed Sep  3 16:14:10 2003 OCC578 Begin
291           Standard_Integer anIntState = (currentstate == TopAbs_IN) ? 0 : 1;
292           //  Modified by skv - Wed Sep  3 16:14:11 2003 OCC578 End
293
294           if(transition != IntCurveSurface_Tangent && face.Orientation()==TopAbs_REVERSED) { 
295             if(transition == IntCurveSurface_In) 
296               transition = IntCurveSurface_Out;
297             else 
298               transition = IntCurveSurface_In;
299           }
300           //----- Insertion du point 
301           if(nbpnt==0) { 
302             IntCurveSurface_IntersectionPoint PPP(pnt,U,V,W,transition);
303             SeqPnt.Append(PPP);
304             //  Modified by skv - Wed Sep  3 16:14:10 2003 OCC578 Begin
305             mySeqState.Append(anIntState);
306             //  Modified by skv - Wed Sep  3 16:14:11 2003 OCC578 End
307           }
308           else { 
309             Standard_Integer i = 1;
310             Standard_Integer b = nbpnt+1;                    
311             while(i<=nbpnt) {
312               const IntCurveSurface_IntersectionPoint& Pnti=SeqPnt.Value(i);
313               Standard_Real wi = Pnti.W();
314               if(wi >= W) { b=i; i=nbpnt; }
315               i++;
316             }
317             IntCurveSurface_IntersectionPoint PPP(pnt,U,V,W,transition);
318             //  Modified by skv - Wed Sep  3 16:14:10 2003 OCC578 Begin
319             //      if(b>nbpnt)          { SeqPnt.Append(PPP); } 
320             //      else if(b>0)             { SeqPnt.InsertBefore(b,PPP); } 
321             if(b>nbpnt) {
322               SeqPnt.Append(PPP);
323               mySeqState.Append(anIntState);
324             } else if(b>0) {
325               SeqPnt.InsertBefore(b,PPP);
326               mySeqState.InsertBefore(b, anIntState);
327             }
328             //  Modified by skv - Wed Sep  3 16:14:11 2003 OCC578 End
329           }
330
331
332           nbpnt++;
333         } 
334       } //-- classifier state is IN or ON
335     } //-- Loop on Intersection points.
336   } //-- HICS.IsDone()
337   else if (HICS.IsDone())
338   {
339     myIsParallel = HICS.IsParallel();
340   }
341 }
342 //=======================================================================
343 //function : Perform
344 //purpose  : 
345 //=======================================================================
346 void IntCurvesFace_Intersector::Perform(const gp_Lin& L,
347                                         const Standard_Real ParMin,
348                                         const Standard_Real ParMax)
349
350   done = Standard_False;
351   if (!myReady)
352   {
353     return;
354   }
355   done = Standard_True;
356   SeqPnt.Clear();
357   mySeqState.Clear();
358   nbpnt = 0;
359   
360   IntCurveSurface_HInter HICS; 
361   Handle(Geom_Line) geomline = new Geom_Line(L);
362   GeomAdaptor_Curve LL(geomline);
363   Handle(GeomAdaptor_HCurve) HLL  = new GeomAdaptor_HCurve(LL);
364   Standard_Real parinf=ParMin;
365   Standard_Real parsup=ParMax;
366   //
367   if(PtrOnPolyhedron == NULL) { 
368     HICS.Perform(HLL,Hsurface);
369   }
370   else { 
371     Intf_Tool bndTool;
372     Bnd_Box   boxLine;
373     bndTool.LinBox
374       (L,
375        ((IntCurveSurface_ThePolyhedronOfHInter *)PtrOnPolyhedron)->Bounding(),
376        boxLine);
377     if(bndTool.NbSegments() == 0) 
378       return;
379     for(Standard_Integer nbseg=1; nbseg<= bndTool.NbSegments(); nbseg++) { 
380       Standard_Real pinf = bndTool.BeginParam(nbseg);
381       Standard_Real psup = bndTool.EndParam(nbseg);
382       Standard_Real pppp = 0.05*(psup-pinf);
383       pinf-=pppp;
384       psup+=pppp;
385       if((psup - pinf)<1e-10) { pinf-=1e-10; psup+=1e-10; } 
386       if(nbseg==1) { parinf=pinf; parsup=psup; }
387       else { 
388         if(parinf>pinf) parinf = pinf;
389         if(parsup<psup) parsup = psup;
390       }
391     }
392     if(parinf>ParMax) { return; } 
393     if(parsup<ParMin) { return; }
394     if(parinf<ParMin) parinf=ParMin;
395     if(parsup>ParMax) parsup=ParMax;
396     if(parinf>(parsup-1e-9)) return; 
397     IntCurveSurface_ThePolygonOfHInter polygon(HLL,
398                                                parinf,
399                                                parsup,
400                                                2);
401 #if OPTIMISATION
402     if(PtrOnBndBounding==NULL) { 
403       PtrOnBndBounding = (Bnd_BoundSortBox *) new Bnd_BoundSortBox();
404       IntCurveSurface_ThePolyhedronOfHInter *thePolyh=
405         (IntCurveSurface_ThePolyhedronOfHInter *)PtrOnPolyhedron;
406       ((Bnd_BoundSortBox *)(PtrOnBndBounding))->
407         Initialize(IntCurveSurface_ThePolyhedronToolOfHInter::Bounding(*thePolyh),
408                    IntCurveSurface_ThePolyhedronToolOfHInter::ComponentsBounding(*thePolyh));
409     }
410     HICS.Perform(HLL,
411                  polygon,
412                  Hsurface,
413                  *((IntCurveSurface_ThePolyhedronOfHInter *)PtrOnPolyhedron),
414                  *((Bnd_BoundSortBox *)PtrOnBndBounding));
415 #else
416     HICS.Perform(HLL,
417                  polygon,
418                  Hsurface,
419                  *((IntCurveSurface_ThePolyhedronOfHInter *)PtrOnPolyhedron));
420 #endif
421   }
422   
423   InternalCall(HICS,parinf,parsup);
424 }
425 //=======================================================================
426 //function : Perform
427 //purpose  : 
428 //=======================================================================
429 void IntCurvesFace_Intersector::Perform(const Handle(Adaptor3d_HCurve)& HCu,
430                                         const Standard_Real ParMin,
431                                         const Standard_Real ParMax) 
432
433   done = Standard_False;
434   if (!myReady)
435   {
436     return;
437   }
438   done = Standard_True;
439   SeqPnt.Clear();
440   //  Modified by skv - Wed Sep  3 16:14:10 2003 OCC578 Begin
441   mySeqState.Clear();
442   //  Modified by skv - Wed Sep  3 16:14:11 2003 OCC578 End
443   nbpnt = 0;
444   IntCurveSurface_HInter            HICS; 
445   
446   //-- 
447   Standard_Real parinf=ParMin;
448   Standard_Real parsup=ParMax;
449
450   if(PtrOnPolyhedron == NULL) { 
451     HICS.Perform(HCu,Hsurface);
452   }
453   else { 
454     parinf = IntCurveSurface_TheHCurveTool::FirstParameter(HCu);
455     parsup = IntCurveSurface_TheHCurveTool::LastParameter(HCu);
456     if(parinf<ParMin) parinf = ParMin;
457     if(parsup>ParMax) parsup = ParMax;
458     if(parinf>(parsup-1e-9)) return; 
459     Standard_Integer nbs;
460     nbs = IntCurveSurface_TheHCurveTool::NbSamples(HCu,parinf,parsup);
461     
462     IntCurveSurface_ThePolygonOfHInter polygon(HCu,
463                                                parinf,
464                                                parsup,
465                                                nbs);
466 #if OPTIMISATION
467     if(PtrOnBndBounding==NULL) { 
468       PtrOnBndBounding = (Bnd_BoundSortBox *) new Bnd_BoundSortBox();
469       IntCurveSurface_ThePolyhedronOfHInter *thePolyh=(IntCurveSurface_ThePolyhedronOfHInter *)PtrOnPolyhedron;
470       ((Bnd_BoundSortBox *)(PtrOnBndBounding))->Initialize(IntCurveSurface_ThePolyhedronToolOfHInter::Bounding(*thePolyh),
471                                                            IntCurveSurface_ThePolyhedronToolOfHInter::ComponentsBounding(*thePolyh));
472     }
473     HICS.Perform(HCu,
474                  polygon,
475                  Hsurface,
476                  *((IntCurveSurface_ThePolyhedronOfHInter *)PtrOnPolyhedron),
477                  *((Bnd_BoundSortBox *)PtrOnBndBounding));
478 #else
479     HICS.Perform(HCu,
480                  polygon,
481                  Hsurface,
482                  *((IntCurveSurface_ThePolyhedronOfHInter *)PtrOnPolyhedron));
483 #endif
484   }
485   InternalCall(HICS,parinf,parsup);
486 }
487
488 //============================================================================
489 Bnd_Box IntCurvesFace_Intersector::Bounding() const {
490   if(PtrOnPolyhedron !=NULL) {
491     return(((IntCurveSurface_ThePolyhedronOfHInter *)PtrOnPolyhedron)->Bounding());
492   }
493   else { 
494     Bnd_Box B;
495     return(B);
496   }
497 }
498 TopAbs_State IntCurvesFace_Intersector::ClassifyUVPoint(const gp_Pnt2d& Puv) const { 
499   TopAbs_State state = myTopolTool->Classify(Puv,1e-7);
500   return(state);
501 }
502 //============================================================================
503 void IntCurvesFace_Intersector::Destroy() { 
504   if(PtrOnPolyhedron !=NULL) { 
505     delete (IntCurveSurface_ThePolyhedronOfHInter *)PtrOnPolyhedron;
506     PtrOnPolyhedron = NULL;
507   }
508   if(PtrOnBndBounding !=NULL) { 
509     delete (Bnd_BoundSortBox *)PtrOnBndBounding;
510     PtrOnBndBounding=NULL;
511   }
512 }
513
514  void IntCurvesFace_Intersector::SetUseBoundToler(Standard_Boolean UseBToler)
515  {
516    myUseBoundTol = UseBToler;
517  }
518
519  Standard_Boolean IntCurvesFace_Intersector::GetUseBoundToler() const
520  {
521    return myUseBoundTol;
522  }
523