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-2012 OPEN CASCADE SAS
   5 //
   6 // The content of this file is subject to the Open CASCADE Technology Public
   7 // License Version 6.5 (the "License"). You may not use the content of this file
   8 // except in compliance with the License. Please obtain a copy of the License
   9 // at http://www.opencascade.org and read it completely before using this file.
  10 //
  11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
  12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
  13 //
  14 // The Original Code and all software distributed under the License is
  15 // distributed on an "AS IS" basis, without warranty of any kind, and the
  16 // Initial Developer hereby disclaims all such warranties, including without
  17 // limitation, any warranties of merchantability, fitness for a particular
  18 // purpose or non-infringement. Please see the License for the specific terms
  19 // and conditions governing the rights and limitations under the License.
  20
  21
  22 #include <BRepLib_FindSurface.ixx>
  23
  24 #include <Precision.hxx>
  25 #include <math_Matrix.hxx>
  26 #include <math_Vector.hxx>
  27 #include <math_Gauss.hxx>
  28
  29 #include <gp_Lin.hxx>
  30 #include <gp_Circ.hxx>
  31 #include <gp_Elips.hxx>
  32 #include <gp_Hypr.hxx>
  33 #include <gp_Parab.hxx>
  34 #include <gp_Ax2.hxx>
  35 #include <gp_Ax3.hxx>
  36 #include <gp_Vec.hxx>
  37 #include <TColgp_SequenceOfPnt.hxx>
  38 #include <TColStd_SequenceOfReal.hxx>
  39 #include <TColgp_Array1OfPnt.hxx>
  40 #include <TColgp_HArray1OfPnt.hxx>
  41 #include <Geom_Plane.hxx>
  42
  43 #include <TopoDS.hxx>
  44 #include <TopExp_Explorer.hxx>
  45 #include <BRep_Tool.hxx>
  46 #include <BRepAdaptor_Curve.hxx>
  47
  48 #include <GeomLib.hxx>
  49 #include <Geom2d_Curve.hxx>
  50 #include <Geom_BezierCurve.hxx>
  51 #include <Geom_BSplineCurve.hxx>
  52
  53 //=======================================================================
  54 //function : Controle
  55 //purpose  :
  56 //=======================================================================
  57 static Standard_Real Controle(const TColgp_SequenceOfPnt& thePoints,
  58                               const Handle(Geom_Plane)& thePlane)
  59 {
  60   Standard_Real dfMaxDist=0.;
  61   Standard_Real a,b,c,d, dist;
  62   Standard_Integer ii;
  63   thePlane->Coefficients(a,b,c,d);
  64   for (ii=1; ii<=thePoints.Length(); ii++) {
  65     const gp_XYZ& xyz = thePoints(ii).XYZ();
  66     dist = Abs(a*xyz.X() + b*xyz.Y() + c*xyz.Z() + d);
  67     if (dist > dfMaxDist)
  68       dfMaxDist = dist;
  69   }
  70
  71   return dfMaxDist;
  72 }
  73
  74 //=======================================================================
  75 //function : BRepLib_FindSurface
  76 //purpose  :
  77 //=======================================================================
  78 BRepLib_FindSurface::BRepLib_FindSurface()
  79 {
  80 }
  81 //=======================================================================
  82 //function : BRepLib_FindSurface
  83 //purpose  :
  84 //=======================================================================
  85 BRepLib_FindSurface::BRepLib_FindSurface(const TopoDS_Shape&    S,
  86                                          const Standard_Real    Tol,
  87                                          const Standard_Boolean OnlyPlane)
  88 {
  89   Init(S,Tol,OnlyPlane);
  90 }
  91 //=======================================================================
  92 //function : Init
  93 //purpose  :
  94 //=======================================================================
  95 void BRepLib_FindSurface::Init(const TopoDS_Shape&    S,
  96                                const Standard_Real    Tol,
  97                                const Standard_Boolean OnlyPlane)
  98 {
  99   myTolerance = Tol;
 100   myTolReached = 0.;
 101   isExisted = Standard_False;
 102   myLocation.Identity();
 103   mySurface.Nullify();
 104
 105   // compute the tolerance
 106   TopExp_Explorer ex;
 107
 108   for (ex.Init(S,TopAbs_EDGE); ex.More(); ex.Next()) {
 109     Standard_Real t = BRep_Tool::Tolerance(TopoDS::Edge(ex.Current()));
 110     if (t > myTolerance) myTolerance = t;
 111   }
 112
 113   // search an existing surface
 114   ex.Init(S,TopAbs_EDGE);
 115   if (!ex.More()) return;    // no edges ....
 116
 117   TopoDS_Edge E = TopoDS::Edge(ex.Current());
 118   Standard_Real f,l,ff,ll;
 119   Handle(Geom2d_Curve) PC,PPC;
 120   Handle(Geom_Surface) SS;
 121   TopLoc_Location L;
 122   Standard_Integer i = 0,j;
 123
 124   // iterate on the surfaces of the first edge
 125   while ( Standard_True) {
 126     i++;
 127     BRep_Tool::CurveOnSurface(E,PC,mySurface,myLocation,f,l,i);
 128     if (mySurface.IsNull()) {
 129       break;
 130     }
 131     // check the other edges
 132     for (ex.Init(S,TopAbs_EDGE); ex.More(); ex.Next()) {
 133       if (!E.IsSame(ex.Current())) {
 134         j = 0;
 135         while (Standard_True) {
 136           j++;
 137           BRep_Tool::CurveOnSurface(TopoDS::Edge(ex.Current()),
 138                                     PPC,SS,L,ff,ll,j);
 139           if (SS.IsNull()) {
 140             break;
 141           }
 142           if (SS == mySurface) {
 143             break;
 144           }
 145           SS.Nullify();
 146         }
 147
 148         if (SS.IsNull()) {
 149           mySurface.Nullify();
 150           break;
 151         }
 152       }
 153     }
 154
 155     // if OnlyPlane, eval if mySurface is a plane.
 156     if ( OnlyPlane && !mySurface.IsNull() )
 157       mySurface = Handle(Geom_Plane)::DownCast(mySurface);
 158
 159     if (!mySurface.IsNull()) break;
 160   }
 161
 162   if (!mySurface.IsNull()) {
 163     isExisted = Standard_True;
 164     return;
 165   }
 166   //
 167   // no existing surface, search a plane
 168   // 07/02/02 akm vvv : (OCC157) changed algorithm
 169   //                    1. Collect the points along all edges of the shape
 170   //                       For each point calculate the WEIGHT = sum of
 171   //                       distances from neighboring points (_only_ same edge)
 172   //                    2. Minimizing the weighed sum of squared deviations
 173   //                       compute coefficients of the sought plane.
 174
 175   TColgp_SequenceOfPnt aPoints;
 176   TColStd_SequenceOfReal aWeight;
 177
 178   // ======================= Step #1
 179   for (ex.Init(S,TopAbs_EDGE); ex.More(); ex.Next())
 180   {
 181     BRepAdaptor_Curve c(TopoDS::Edge(ex.Current()));
 182
 183     Standard_Real dfUf = c.FirstParameter();
 184     Standard_Real dfUl = c.LastParameter();
 185     if (IsEqual(dfUf,dfUl)) {
 186       // Degenerate
 187       continue;
 188     }
 189     Standard_Integer iNbPoints=0;
 190
 191     // Add the points with weights to the sequences
 192     switch (c.GetType())
 193     {
 194     case GeomAbs_BezierCurve:
 195       {
 196         // Put all poles for bezier
 197         Handle(Geom_BezierCurve) GC = c.Bezier();
 198         Standard_Integer iNbPol = GC->NbPoles();
 199         if ( iNbPol < 2)
 200           // Degenerate
 201           continue;
 202         else
 203         {
 204           Handle(TColgp_HArray1OfPnt) aPoles = new (TColgp_HArray1OfPnt) (1, iNbPol);
 205           GC->Poles(aPoles->ChangeArray1());
 206           gp_Pnt aPolePrev = aPoles->Value(1), aPoleNext;
 207           Standard_Real dfDistPrev = 0., dfDistNext;
 208           for (Standard_Integer iPol=1; iPol<=iNbPol; iPol++)
 209           {
 210             if (iPol<iNbPol)
 211             {
 212               aPoleNext = aPoles->Value(iPol+1);
 213               dfDistNext = aPolePrev.Distance(aPoleNext);
 214             }
 215             else
 216               dfDistNext = 0.;
 217             aPoints.Append (aPolePrev);
 218             aWeight.Append (dfDistPrev+dfDistNext);
 219             dfDistPrev = dfDistNext;
 220             aPolePrev = aPoleNext;
 221           }
 222         }
 223       }
 224       break;
 225     case GeomAbs_BSplineCurve:
 226       {
 227         // Put all poles for bspline
 228         Handle(Geom_BSplineCurve) GC = c.BSpline();
 229         Standard_Integer iNbPol = GC->NbPoles();
 230         if ( iNbPol < 2)
 231           // Degenerate
 232           continue;
 233         else
 234         {
 235           Handle(TColgp_HArray1OfPnt) aPoles = new (TColgp_HArray1OfPnt) (1, iNbPol);
 236           GC->Poles(aPoles->ChangeArray1());
 237           gp_Pnt aPolePrev = aPoles->Value(1), aPoleNext;
 238           Standard_Real dfDistPrev = 0., dfDistNext;
 239           for (Standard_Integer iPol=1; iPol<=iNbPol; iPol++)
 240           {
 241             if (iPol<iNbPol)
 242             {
 243               aPoleNext = aPoles->Value(iPol+1);
 244               dfDistNext = aPolePrev.Distance(aPoleNext);
 245             }
 246             else
 247               dfDistNext = 0.;
 248             aPoints.Append (aPolePrev);
 249             aWeight.Append (dfDistPrev+dfDistNext);
 250             dfDistPrev = dfDistNext;
 251             aPolePrev = aPoleNext;
 252           }
 253         }
 254       }
 255       break;
 256
 257     case GeomAbs_Line:
 258     case GeomAbs_Circle:
 259     case GeomAbs_Ellipse:
 260     case GeomAbs_Hyperbola:
 261     case GeomAbs_Parabola:
 262       if (c.GetType() == GeomAbs_Line)
 263         // Two points on straight segment
 264         iNbPoints=2;
 265       else
 266         // Four points on otheranalitical curves
 267         iNbPoints=4;
 268     default:
 269       {
 270         // Put some points on other curves
 271         if (iNbPoints==0)
 272           iNbPoints = 15 + c.NbIntervals(GeomAbs_C3);
 273         Standard_Real dfDelta = (dfUl-dfUf)/(iNbPoints-1);
 274         Standard_Integer iPoint;
 275         Standard_Real dfU;
 276         gp_Pnt aPointPrev = c.Value(dfUf), aPointNext;
 277         Standard_Real dfDistPrev = 0., dfDistNext;
 278         for (iPoint=1, dfU=dfUf+dfDelta;
 279              iPoint<=iNbPoints;
 280              iPoint++, dfU+=dfDelta)
 281         {
 282           if (iPoint<iNbPoints)
 283           {
 284             aPointNext = c.Value(dfU);
 285             dfDistNext = aPointPrev.Distance(aPointNext);
 286           }
 287           else
 288             dfDistNext = 0.;
 289           aPoints.Append (aPointPrev);
 290           aWeight.Append (dfDistPrev+dfDistNext);
 291           dfDistPrev = dfDistNext;
 292           aPointPrev = aPointNext;
 293         }
 294       } // default:
 295     } // switch (c.GetType()) ...
 296   } // for (ex.Init(S,TopAbs_EDGE); ex.More() && control; ex.Next()) ...
 297
 298   if (aPoints.Length() < 3) {
 299     return;
 300   }
 301
 302   // ======================= Step #2
 303   myLocation.Identity();
 304   Standard_Integer iPoint;
 305   math_Matrix aMat (1,3,1,3, 0.);
 306   math_Vector aVec (1,3, 0.);
 307   // Find the barycenter and normalize weights
 308   Standard_Real dfMaxWeight=0.;
 309   gp_XYZ aBaryCenter(0.,0.,0.);
 310   Standard_Real dfSumWeight=0.;
 311   for (iPoint=1; iPoint<=aPoints.Length(); iPoint++)  {
 312     Standard_Real dfW = aWeight(iPoint);
 313     aBaryCenter += dfW*aPoints(iPoint).XYZ();
 314     dfSumWeight += dfW;
 315     if (dfW > dfMaxWeight) {
 316       dfMaxWeight = dfW;
 317     }
 318   }
 319   aBaryCenter /= dfSumWeight;
 320
 321   // Fill the matrix and the right vector
 322   for (iPoint=1; iPoint<=aPoints.Length(); iPoint++)  {
 323     gp_XYZ p=aPoints(iPoint).XYZ()-aBaryCenter;
 324     Standard_Real w=aWeight(iPoint)/dfMaxWeight;
 325     aMat(1,1)+=w*p.X()*p.X();
 326         aMat(1,2)+=w*p.X()*p.Y();
 327             aMat(1,3)+=w*p.X()*p.Z();
 328     aMat(2,1)+=w*p.Y()*p.X();
 329         aMat(2,2)+=w*p.Y()*p.Y();
 330             aMat(2,3)+=w*p.Y()*p.Z();
 331     aMat(3,1)+=w*p.Z()*p.X();
 332         aMat(3,2)+=w*p.Z()*p.Y();
 333             aMat(3,3)+=w*p.Z()*p.Z();
 334     aVec(1) -= w*p.X();
 335     aVec(2) -= w*p.Y();
 336     aVec(3) -= w*p.Z();
 337   }
 338
 339   // Solve the system of equations to get plane coefficients
 340   math_Gauss aSolver(aMat);
 341   Standard_Boolean isSolved = aSolver.IsDone();
 342   //
 343   //  let us be more tolerant (occ415)
 344   Standard_Real dfDist = RealLast();
 345   Handle(Geom_Plane) aPlane;
 346   //
 347   if (isSolved)  {
 348     aSolver.Solve(aVec);
 349     if (aVec.Norm2()<gp::Resolution()) {
 350       isSolved = Standard_False;
 351     }
 352   }
 353   //
 354   if (isSolved)  {
 355     aPlane = new Geom_Plane(aBaryCenter,gp_Dir(aVec(1),aVec(2),aVec(3)));
 356     dfDist = Controle (aPoints, aPlane);
 357   }
 358   //
 359   if (!isSolved || myTolerance < dfDist)  {
 360     gp_Pnt aFirstPnt=aPoints(1);
 361     for (iPoint=2; iPoint<=aPoints.Length(); iPoint++)  {
 362       gp_Vec aDir(aFirstPnt,aPoints(iPoint));
 363       Standard_Real dfSide=aDir.Magnitude();
 364       if (dfSide<myTolerance) {
 365         continue; // degeneration
 366       }
 367       for (Standard_Integer iP1=iPoint+1; iP1<=aPoints.Length(); iP1++) {
 368         gp_Vec aCross = gp_Vec(aFirstPnt,aPoints(iP1)) ^ aDir ;
 369         if (aCross.Magnitude() > dfSide*myTolerance) {
 370           Handle(Geom_Plane) aPlane2 = new Geom_Plane(aFirstPnt, aCross);
 371           Standard_Real dfDist2 = Controle (aPoints, aPlane2);
 372           if (dfDist2 < myTolerance)  {
 373             myTolReached = dfDist2;
 374             mySurface = aPlane2;
 375             return;
 376           }
 377           if (dfDist2 < dfDist)  {
 378             dfDist = dfDist2;
 379             aPlane = aPlane2;
 380           }
 381         }
 382       }
 383     }
 384   }
 385   //
 386   //XXf
 387   //static Standard_Real weakness = 5.0;
 388   Standard_Real weakness = 5.0;
 389   //XXf
 390   if(dfDist <= myTolerance || dfDist < myTolerance*weakness && Tol<0) {
 391     //XXf
 392     //myTolReached = dfDist;
 393     //XXt
 394     mySurface = aPlane;
 395   }
 396   //XXf
 397   myTolReached = dfDist;
 398   //XXt
 399 }
 400 //=======================================================================
 401 //function : Found
 402 //purpose  :
 403 //=======================================================================
 404 Standard_Boolean BRepLib_FindSurface::Found() const
 405 {
 406   return !mySurface.IsNull();
 407 }
 408 //=======================================================================
 409 //function : Surface
 410 //purpose  :
 411 //=======================================================================
 412 Handle(Geom_Surface) BRepLib_FindSurface::Surface() const
 413 {
 414   return mySurface;
 415 }
 416 //=======================================================================
 417 //function : Tolerance
 418 //purpose  :
 419 //=======================================================================
 420 Standard_Real BRepLib_FindSurface::Tolerance() const
 421 {
 422   return myTolerance;
 423 }
 424 //=======================================================================
 425 //function : ToleranceReached
 426 //purpose  :
 427 //=======================================================================
 428 Standard_Real BRepLib_FindSurface::ToleranceReached() const
 429 {
 430   return myTolReached;
 431 }
 432 //=======================================================================
 433 //function : Existed
 434 //purpose  :
 435 //=======================================================================
 436 Standard_Boolean BRepLib_FindSurface::Existed() const
 437 {
 438   return isExisted;
 439 }
 440 //=======================================================================
 441 //function : Location
 442 //purpose  :
 443 //=======================================================================
 444 TopLoc_Location BRepLib_FindSurface::Location() const
 445 {
 446   return myLocation;
 447 }