0024255: Regressions in test cases on OCCT vc9 win64 Release
[occt.git] / src / BRepLib / BRepLib_FindSurface.cxx
1 // Created on: 1994-07-22
2 // Created by: Remi LEQUETTE
3 // Copyright (c) 1994-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 #include <BRepLib_FindSurface.ixx>
18
19 #include <Precision.hxx>
20 #include <math_Matrix.hxx>
21 #include <math_Vector.hxx>
22 #include <math_Jacobi.hxx>
23
24 #include <gp_Lin.hxx>
25 #include <gp_Circ.hxx>
26 #include <gp_Elips.hxx>
27 #include <gp_Hypr.hxx>
28 #include <gp_Parab.hxx>
29 #include <gp_Ax2.hxx>
30 #include <gp_Ax3.hxx>
31 #include <gp_Vec.hxx>
32 #include <TColgp_SequenceOfPnt.hxx>
33 #include <TColStd_SequenceOfReal.hxx>
34 #include <TColgp_Array1OfPnt.hxx>
35 #include <TColgp_HArray1OfPnt.hxx>
36 #include <Geom_Plane.hxx>
37
38 #include <BRepAdaptor_Curve.hxx>
39 #include <BRepAdaptor_Surface.hxx>
40 #include <BRepLib_MakeFace.hxx>
41 #include <BRepTools_WireExplorer.hxx>
42 #include <BRep_Tool.hxx>
43 #include <TopExp.hxx>
44 #include <TopExp_Explorer.hxx>
45 #include <TopoDS_Vertex.hxx>
46 #include <TopoDS_Wire.hxx>
47 #include <TopoDS.hxx>
48 #include <BRep_Builder.hxx>
49 #include <BRepTopAdaptor_FClass2d.hxx>
50
51 #include <GeomLib.hxx>
52 #include <Geom2d_Curve.hxx>
53 #include <Geom_BezierCurve.hxx>
54 #include <Geom_BSplineCurve.hxx>
55 #include <Geom_RectangularTrimmedSurface.hxx>
56 #include <Standard_ErrorHandler.hxx>
57
58 //=======================================================================
59 //function : Controle
60 //purpose  : 
61 //=======================================================================
62 static Standard_Real Controle(const TColgp_SequenceOfPnt& thePoints,
63                               const Handle(Geom_Plane)& thePlane)
64 {
65   Standard_Real dfMaxDist=0.;
66   Standard_Real a,b,c,d, dist;
67   Standard_Integer ii;
68   thePlane->Coefficients(a,b,c,d);
69   for (ii=1; ii<=thePoints.Length(); ii++) {
70     const gp_XYZ& xyz = thePoints(ii).XYZ();
71     dist = Abs(a*xyz.X() + b*xyz.Y() + c*xyz.Z() + d);
72     if (dist > dfMaxDist)
73       dfMaxDist = dist;
74   }
75
76   return dfMaxDist;
77 }
78 //=======================================================================
79 //function : Is2DConnected
80 //purpose  : Return true if the last vertex of theEdge1 coincides with
81 //           the first vertex of theEdge2 in parametric space of theFace
82 //=======================================================================
83 inline static Standard_Boolean Is2DConnected(const TopoDS_Edge& theEdge1,
84                                              const TopoDS_Edge& theEdge2,
85                                              const Handle(Geom_Surface)& theSurface,
86                                              const TopLoc_Location& theLocation)
87 {
88   Standard_Real f,l;
89   //TopLoc_Location aLoc;
90   Handle(Geom2d_Curve) aCurve;
91   gp_Pnt2d p1, p2;
92
93   // get 2D points
94   aCurve=BRep_Tool::CurveOnSurface(theEdge1, theSurface, theLocation, f, l);
95   p1       = aCurve->Value( theEdge1.Orientation() == TopAbs_FORWARD ? l : f );
96   aCurve=BRep_Tool::CurveOnSurface(theEdge2, theSurface, theLocation, f, l);
97   p2       = aCurve->Value(  theEdge2.Orientation() == TopAbs_FORWARD ? f : l );
98
99   // compare 2D points
100   GeomAdaptor_Surface aSurface( theSurface );
101   TopoDS_Vertex    aV = TopExp::FirstVertex( theEdge2, Standard_True );
102   Standard_Real tol3D = BRep_Tool::Tolerance( aV );
103   Standard_Real tol2D = aSurface.UResolution( tol3D ) + aSurface.VResolution( tol3D );
104   Standard_Real dist2 = p1.SquareDistance( p2 );
105   return dist2 < tol2D * tol2D;
106 }
107
108 //=======================================================================
109 //function : Is2DClosed
110 //purpose  : Return true if edges of theShape form a closed wire in
111 //           parametric space of theSurface
112 //=======================================================================
113
114 static Standard_Boolean Is2DClosed(const TopoDS_Shape&         theShape,
115                                    const Handle(Geom_Surface)& theSurface,
116                                    const TopLoc_Location& theLocation)
117 {
118   try
119   {
120     // get a wire theShape 
121     TopExp_Explorer aWireExp( theShape, TopAbs_WIRE );
122     if ( !aWireExp.More() ) {
123       return Standard_False;
124     }
125     TopoDS_Wire aWire = TopoDS::Wire( aWireExp.Current() );
126     // a tmp face
127     TopoDS_Face aTmpFace = BRepLib_MakeFace( theSurface, Precision::PConfusion() );
128
129     // check topological closeness using wire explorer, if the wire is not closed
130     // the 1st and the last vertices of wire are different
131     BRepTools_WireExplorer aWireExplorer( aWire, aTmpFace );
132     if ( !aWireExplorer.More()) {
133       return Standard_False;
134     }
135     // remember the 1st and the last edges of aWire
136     TopoDS_Edge aFisrtEdge = aWireExplorer.Current(), aLastEdge = aFisrtEdge;
137     // check if edges connected topologically (that is assured by BRepTools_WireExplorer)
138     // are connected in 2D
139     TopoDS_Edge aPrevEdge = aFisrtEdge;
140     for ( aWireExplorer.Next(); aWireExplorer.More(); aWireExplorer.Next() )    {
141       aLastEdge = aWireExplorer.Current();
142       if ( !Is2DConnected( aPrevEdge, aLastEdge, theSurface, theLocation)) { 
143         return false;
144       }
145       aPrevEdge = aLastEdge;
146     }
147     // wire is closed if ( 1st vertex of aFisrtEdge ) ==
148     // ( last vertex of aLastEdge ) in 2D
149     TopoDS_Vertex aV1 = TopExp::FirstVertex( aFisrtEdge, Standard_True );
150     TopoDS_Vertex aV2 = TopExp::LastVertex( aLastEdge, Standard_True );
151     return ( aV1.IsSame( aV2 ) && Is2DConnected( aLastEdge, aFisrtEdge, theSurface, theLocation));
152   }
153   catch ( Standard_Failure )  {
154     return Standard_False;
155   }
156 }
157 //=======================================================================
158 //function : BRepLib_FindSurface
159 //purpose  : 
160 //=======================================================================
161 BRepLib_FindSurface::BRepLib_FindSurface() 
162 {
163 }
164 //=======================================================================
165 //function : BRepLib_FindSurface
166 //purpose  : 
167 //=======================================================================
168 BRepLib_FindSurface::BRepLib_FindSurface(const TopoDS_Shape&    S, 
169                                          const Standard_Real    Tol,
170                                          const Standard_Boolean OnlyPlane,
171                                          const Standard_Boolean OnlyClosed)
172 {
173   Init(S,Tol,OnlyPlane,OnlyClosed);
174 }
175 //=======================================================================
176 //function : Init
177 //purpose  : 
178 //=======================================================================
179 void BRepLib_FindSurface::Init(const TopoDS_Shape&    S, 
180                                                  const Standard_Real    Tol,
181                                                  const Standard_Boolean OnlyPlane,
182                                const Standard_Boolean OnlyClosed)
183 {
184   myTolerance = Tol;
185   myTolReached = 0.;
186   isExisted = Standard_False;
187   myLocation.Identity();
188   mySurface.Nullify();
189
190   // compute the tolerance
191   TopExp_Explorer ex;
192
193   for (ex.Init(S,TopAbs_EDGE); ex.More(); ex.Next()) {
194     Standard_Real t = BRep_Tool::Tolerance(TopoDS::Edge(ex.Current()));
195     if (t > myTolerance) myTolerance = t;
196   }
197
198   // search an existing surface
199   ex.Init(S,TopAbs_EDGE);
200   if (!ex.More()) return;    // no edges ....
201
202   TopoDS_Edge E = TopoDS::Edge(ex.Current());
203   Standard_Real f,l,ff,ll;
204   Handle(Geom2d_Curve) PC,PPC;
205   Handle(Geom_Surface) SS;
206   TopLoc_Location L;
207   Standard_Integer i = 0,j;
208
209   // iterate on the surfaces of the first edge
210   for(;;) {
211     i++;
212     BRep_Tool::CurveOnSurface(E,PC,mySurface,myLocation,f,l,i);
213     if (mySurface.IsNull()) {
214       break;
215     }
216     // check the other edges
217     for (ex.Init(S,TopAbs_EDGE); ex.More(); ex.Next()) {
218       if (!E.IsSame(ex.Current())) {
219         j = 0;
220         for(;;) {
221           j++;
222           BRep_Tool::CurveOnSurface(TopoDS::Edge(ex.Current()),
223             PPC,SS,L,ff,ll,j);
224           if (SS.IsNull()) {
225             break;
226           }
227           if (SS == mySurface) {
228             break;
229           }
230           SS.Nullify();
231         }
232
233         if (SS.IsNull()) {
234           mySurface.Nullify();
235           break;
236         }
237       }
238     }
239
240     // if OnlyPlane, eval if mySurface is a plane.
241     if ( OnlyPlane && !mySurface.IsNull() ) 
242     {
243       if ( mySurface->IsKind( STANDARD_TYPE(Geom_RectangularTrimmedSurface)))
244         mySurface = Handle(Geom_RectangularTrimmedSurface)::DownCast(mySurface)->BasisSurface();
245       mySurface = Handle(Geom_Plane)::DownCast(mySurface);
246     }
247
248     if (!mySurface.IsNull())
249       // if S is e.g. the bottom face of a cylinder, mySurface can be the
250       // lateral (cylindrical) face of the cylinder; reject an improper mySurface
251       if ( !OnlyClosed || Is2DClosed( S, mySurface, myLocation ))
252         break;
253   }
254
255   if (!mySurface.IsNull()) {
256     isExisted = Standard_True;
257     return;
258   }
259   //
260   // no existing surface, search a plane
261   // 07/02/02 akm vvv : (OCC157) changed algorithm
262   //                    1. Collect the points along all edges of the shape
263   //                       For each point calculate the WEIGHT = sum of
264   //                       distances from neighboring points (_only_ same edge)
265   //                    2. Minimizing the weighed sum of squared deviations
266   //                       compute coefficients of the sought plane.
267
268   TColgp_SequenceOfPnt aPoints;
269   TColStd_SequenceOfReal aWeight;
270
271   // ======================= Step #1
272   for (ex.Init(S,TopAbs_EDGE); ex.More(); ex.Next()) 
273   {
274     BRepAdaptor_Curve c(TopoDS::Edge(ex.Current()));
275
276     Standard_Real dfUf = c.FirstParameter();
277     Standard_Real dfUl = c.LastParameter();
278     if (IsEqual(dfUf,dfUl)) {
279       // Degenerate
280       continue;
281     }
282     Standard_Integer iNbPoints=0;
283
284     // Add the points with weights to the sequences
285     switch (c.GetType()) 
286     {
287     case GeomAbs_BezierCurve:
288       {
289         // Put all poles for bezier
290         Handle(Geom_BezierCurve) GC = c.Bezier();
291         Standard_Integer iNbPol = GC->NbPoles();
292         Standard_Real tf = GC->FirstParameter();
293         Standard_Real tl = GC->LastParameter();
294         Standard_Real r =  (dfUl - dfUf) / (tl - tf);
295         r *= iNbPol;
296         if ( iNbPol < 2 || r < 1.)
297           // Degenerate
298           continue;
299         else
300         {
301           Handle(TColgp_HArray1OfPnt) aPoles = new (TColgp_HArray1OfPnt) (1, iNbPol);
302           GC->Poles(aPoles->ChangeArray1());
303           gp_Pnt aPolePrev = aPoles->Value(1), aPoleNext;
304           Standard_Real dfDistPrev = 0., dfDistNext;
305           for (Standard_Integer iPol=1; iPol<=iNbPol; iPol++)
306           {
307             if (iPol<iNbPol)
308             {
309               aPoleNext = aPoles->Value(iPol+1);
310               dfDistNext = aPolePrev.Distance(aPoleNext);
311             }
312             else
313               dfDistNext = 0.;
314             aPoints.Append (aPolePrev);
315             aWeight.Append (dfDistPrev+dfDistNext);
316             dfDistPrev = dfDistNext;
317             aPolePrev = aPoleNext;
318           }
319         }
320       }
321       break;
322     case GeomAbs_BSplineCurve:
323       {
324         // Put all poles for bspline
325         Handle(Geom_BSplineCurve) GC = c.BSpline();
326         Standard_Integer iNbPol = GC->NbPoles();
327         Standard_Real tf = GC->FirstParameter();
328         Standard_Real tl = GC->LastParameter();
329         Standard_Real r =  (dfUl - dfUf) / (tl - tf);
330         r *= iNbPol;
331         if ( iNbPol < 2 || r < 1.)
332           // Degenerate
333           continue;
334         else
335         {
336           Handle(TColgp_HArray1OfPnt) aPoles = new (TColgp_HArray1OfPnt) (1, iNbPol);
337           GC->Poles(aPoles->ChangeArray1());
338           gp_Pnt aPolePrev = aPoles->Value(1), aPoleNext;
339           Standard_Real dfDistPrev = 0., dfDistNext;
340           for (Standard_Integer iPol=1; iPol<=iNbPol; iPol++)
341           {
342             if (iPol<iNbPol)
343             {
344               aPoleNext = aPoles->Value(iPol+1);
345               dfDistNext = aPolePrev.Distance(aPoleNext);
346             }
347             else
348               dfDistNext = 0.;
349             aPoints.Append (aPolePrev);
350             aWeight.Append (dfDistPrev+dfDistNext);
351             dfDistPrev = dfDistNext;
352             aPolePrev = aPoleNext;
353           }
354         }
355       }
356       break;
357
358     case GeomAbs_Line:
359     case GeomAbs_Circle:
360     case GeomAbs_Ellipse:
361     case GeomAbs_Hyperbola:
362     case GeomAbs_Parabola:
363       if (c.GetType() == GeomAbs_Line)
364         // Two points on straight segment
365         iNbPoints=2;
366       else
367         // Four points on otheranalitical curves
368         iNbPoints=4;
369     default:
370       { 
371         // Put some points on other curves
372         if (iNbPoints==0)
373           iNbPoints = 15 + c.NbIntervals(GeomAbs_C3);
374         Standard_Real dfDelta = (dfUl-dfUf)/(iNbPoints-1);
375         Standard_Integer iPoint;
376         Standard_Real dfU;
377         gp_Pnt aPointPrev = c.Value(dfUf), aPointNext;
378         Standard_Real dfDistPrev = 0., dfDistNext;
379         for (iPoint=1, dfU=dfUf+dfDelta; 
380           iPoint<=iNbPoints; 
381           iPoint++, dfU+=dfDelta) 
382         {
383           if (iPoint<iNbPoints)
384           {
385             aPointNext = c.Value(dfU);
386             dfDistNext = aPointPrev.Distance(aPointNext);
387           }
388           else
389             dfDistNext = 0.;
390           aPoints.Append (aPointPrev);
391           aWeight.Append (dfDistPrev+dfDistNext);
392           dfDistPrev = dfDistNext;
393           aPointPrev = aPointNext;
394         }
395       } // default:
396     } // switch (c.GetType()) ...
397   } // for (ex.Init(S,TopAbs_EDGE); ex.More() && control; ex.Next()) ...
398
399   if (aPoints.Length() < 3) {
400     return;
401   }
402
403   // ======================= Step #2
404   myLocation.Identity();
405   Standard_Integer iPoint;
406   math_Matrix aMat (1,3,1,3, 0.);
407   math_Vector aVec (1,3, 0.);
408   // Find the barycenter and normalize weights 
409   Standard_Real dfMaxWeight=0.;
410   gp_XYZ aBaryCenter(0.,0.,0.);
411   Standard_Real dfSumWeight=0.;
412   for (iPoint=1; iPoint<=aPoints.Length(); iPoint++)  {
413     Standard_Real dfW = aWeight(iPoint);
414     aBaryCenter += dfW*aPoints(iPoint).XYZ();
415     dfSumWeight += dfW;
416     if (dfW > dfMaxWeight) {
417       dfMaxWeight = dfW;
418     }
419   }
420   aBaryCenter /= dfSumWeight;
421
422   // Fill the matrix and the right vector
423   for (iPoint=1; iPoint<=aPoints.Length(); iPoint++)  {
424     gp_XYZ p=aPoints(iPoint).XYZ()-aBaryCenter;
425     Standard_Real w=aWeight(iPoint)/dfMaxWeight;
426     aMat(1,1)+=w*p.X()*p.X(); 
427     aMat(1,2)+=w*p.X()*p.Y(); 
428     aMat(1,3)+=w*p.X()*p.Z();
429     //  
430     aMat(2,2)+=w*p.Y()*p.Y();  
431     aMat(2,3)+=w*p.Y()*p.Z();
432     //  
433     aMat(3,3)+=w*p.Z()*p.Z();
434   }
435   aMat(2,1) = aMat(1,2);
436   aMat(3,1) = aMat(1,3);
437   aMat(3,2) = aMat(2,3);
438   //
439   math_Jacobi anEignval(aMat);
440   math_Vector anEVals(1,3);
441   Standard_Boolean isSolved = anEignval.IsDone();
442   Standard_Integer isol = 0;
443   if(isSolved)
444   {
445     anEVals = anEignval.Values();
446     //We need vector with eigenvalue ~ 0.
447     Standard_Real anEMin = RealLast();
448     Standard_Real anEMax = -anEMin;
449     for(i = 1; i <= 3; ++i)
450     {
451       Standard_Real anE = Abs(anEVals(i));
452       if(anEMin > anE)
453       {
454         anEMin = anE;
455         isol = i;
456       }
457       if(anEMax < anE)
458       {
459         anEMax = anE;
460       }
461     }
462     
463     if(isol == 0)
464     {
465       isSolved = Standard_False;
466     }
467     else
468     {
469       Standard_Real eps = Epsilon(anEMax);
470       if(anEMin <= eps)
471       {
472         anEignval.Vector(isol, aVec);
473       }
474       else
475       {
476         //try using vector product of other axes
477         Standard_Integer ind[2] = {0,0};
478         for(i = 1; i <= 3; ++i)
479         {
480           if(i == isol)
481           {
482             continue;
483           }
484           if(ind[0] == 0)
485           {
486             ind[0] = i;
487             continue;
488           }
489           if(ind[1] == 0)
490           {
491             ind[1] = i;
492             continue;
493           }
494         }
495         math_Vector aVec1(1, 3, 0.), aVec2(1, 3, 0.);
496         anEignval.Vector(ind[0], aVec1);
497         anEignval.Vector(ind[1], aVec2);
498         gp_Vec aV1(aVec1(1), aVec1(2), aVec1(3));
499         gp_Vec aV2(aVec2(1), aVec2(2), aVec2(3));
500         gp_Vec aN = aV1^ aV2;
501         aVec(1) = aN.X();
502         aVec(2) = aN.Y();
503         aVec(3) = aN.Z();
504       }
505       if (aVec.Norm2() < gp::Resolution()) {
506         isSolved = Standard_False;
507       }
508     }
509   }
510     
511   //
512   //  let us be more tolerant (occ415)
513   Standard_Real dfDist = RealLast();
514   Handle(Geom_Plane) aPlane;
515   //
516   if (isSolved)  {
517     //Plane normal can have two directions, direction is chosen
518     //according to direction of eigenvector
519     gp_Vec anN(aVec(1), aVec(2), aVec(3));
520     aPlane = new Geom_Plane(aBaryCenter,anN);
521     dfDist = Controle (aPoints, aPlane);
522   }
523   //
524   if (!isSolved || myTolerance < dfDist)  {
525     gp_Pnt aFirstPnt=aPoints(1);
526     for (iPoint=2; iPoint<=aPoints.Length(); iPoint++)  {
527       gp_Vec aDir(aFirstPnt,aPoints(iPoint));
528       Standard_Real dfSide=aDir.Magnitude();
529       if (dfSide<myTolerance) {
530         continue; // degeneration
531       }
532       for (Standard_Integer iP1=iPoint+1; iP1<=aPoints.Length(); iP1++) {
533
534         gp_Vec aCross = gp_Vec(aFirstPnt,aPoints(iP1)) ^ aDir ;
535
536         if (aCross.Magnitude() > dfSide*myTolerance) {
537           Handle(Geom_Plane) aPlane2 = new Geom_Plane(aBaryCenter, aCross);
538           Standard_Real dfDist2 = Controle (aPoints, aPlane2);
539           if (dfDist2 < myTolerance)  {
540             myTolReached = dfDist2;
541             mySurface = aPlane2;
542             return;
543           }
544           if (dfDist2 < dfDist)  {
545             dfDist = dfDist2;
546             aPlane = aPlane2;
547           }
548         }
549       }
550     }
551   }
552   //
553   //XXf
554   //static Standard_Real weakness = 5.0;
555   Standard_Real weakness = 5.0;
556   //XXf
557   if(dfDist <= myTolerance || (dfDist < myTolerance*weakness && Tol<0)) { 
558     //XXf 
559     //myTolReached = dfDist;
560     //XXt
561     mySurface = aPlane;
562     //If S is wire, try to orient surface according to orientation of wire.
563     if(S.ShapeType() == TopAbs_WIRE && S.Closed())
564     {
565        //
566       TopoDS_Wire aW = TopoDS::Wire(S);
567       TopoDS_Face aTmpFace = BRepLib_MakeFace(mySurface, Precision::Confusion());
568       BRep_Builder BB;
569       BB.Add(aTmpFace, aW);
570       BRepTopAdaptor_FClass2d FClass(aTmpFace, 0.);
571       if ( FClass.PerformInfinitePoint() == TopAbs_IN ) 
572       {
573         gp_Dir aN = aPlane->Position().Direction();
574         aN.Reverse();
575         mySurface = new Geom_Plane(aPlane->Position().Location(), aN);
576       }
577
578     }
579   }
580   //XXf
581   myTolReached = dfDist;
582   //XXt
583 }
584 //=======================================================================
585 //function : Found
586 //purpose  : 
587 //=======================================================================
588 Standard_Boolean BRepLib_FindSurface::Found() const 
589 {
590   return !mySurface.IsNull();
591 }
592 //=======================================================================
593 //function : Surface
594 //purpose  : 
595 //=======================================================================
596 Handle(Geom_Surface) BRepLib_FindSurface::Surface() const 
597 {
598   return mySurface;
599 }
600 //=======================================================================
601 //function : Tolerance
602 //purpose  : 
603 //=======================================================================
604 Standard_Real BRepLib_FindSurface::Tolerance() const 
605 {
606   return myTolerance;
607 }
608 //=======================================================================
609 //function : ToleranceReached
610 //purpose  : 
611 //=======================================================================
612 Standard_Real BRepLib_FindSurface::ToleranceReached() const 
613 {
614   return myTolReached;
615 }
616 //=======================================================================
617 //function : Existed
618 //purpose  : 
619 //=======================================================================
620 Standard_Boolean BRepLib_FindSurface::Existed() const 
621 {
622   return isExisted;
623 }
624 //=======================================================================
625 //function : Location
626 //purpose  : 
627 //=======================================================================
628 TopLoc_Location BRepLib_FindSurface::Location() const
629 {
630   return myLocation;
631 }
632