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