src/BRepMesh/BRepMesh_GeomTool.cxx

   1 // Created on: 1993-09-29
   2 // Created by: Isabelle GRIGNON
   3 // Copyright (c) 1993-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 <BRepMesh_GeomTool.hxx>
  18
  19 #include <BRepMesh_DefaultRangeSplitter.hxx>
  20
  21 #include <TopAbs_Orientation.hxx>
  22 #include <CSLib.hxx>
  23 #include <Precision.hxx>
  24 #include <Adaptor3d_IsoCurve.hxx>
  25 #include <Adaptor3d_HCurve.hxx>
  26 #include <BRepAdaptor_Curve.hxx>
  27 #include <BRepAdaptor_HSurface.hxx>
  28 #include <Geom2d_Curve.hxx>
  29 #include <BRep_Tool.hxx>
  30
  31 namespace
  32 {
  33   void ComputeErrFactors (const Standard_Real               theDeflection,
  34                           const Handle(Adaptor3d_HSurface)& theFace,
  35                           Standard_Real&                    theErrFactorU,
  36                           Standard_Real&                    theErrFactorV)
  37   {
  38     theErrFactorU = theDeflection * 10.;
  39     theErrFactorV = theDeflection * 10.;
  40
  41     switch (theFace->GetType ())
  42     {
  43     case GeomAbs_Cylinder:
  44     case GeomAbs_Cone:
  45     case GeomAbs_Sphere:
  46     case GeomAbs_Torus:
  47       break;
  48
  49     case GeomAbs_SurfaceOfExtrusion:
  50     case GeomAbs_SurfaceOfRevolution:
  51       {
  52         Handle(Adaptor3d_HCurve) aCurve = theFace->BasisCurve ();
  53         if (aCurve->GetType () == GeomAbs_BSplineCurve && aCurve->Degree () > 2)
  54         {
  55           theErrFactorV /= (aCurve->Degree () * aCurve->NbKnots ());
  56         }
  57         break;
  58       }
  59     case GeomAbs_BezierSurface:
  60       {
  61         if (theFace->UDegree () > 2)
  62         {
  63           theErrFactorU /= (theFace->UDegree ());
  64         }
  65         if (theFace->VDegree () > 2)
  66         {
  67           theErrFactorV /= (theFace->VDegree ());
  68         }
  69         break;
  70       }
  71     case GeomAbs_BSplineSurface:
  72       {
  73         if (theFace->UDegree () > 2)
  74         {
  75           theErrFactorU /= (theFace->UDegree () * theFace->NbUKnots ());
  76         }
  77         if (theFace->VDegree () > 2)
  78         {
  79           theErrFactorV /= (theFace->VDegree () *  theFace->NbVKnots ());
  80         }
  81         break;
  82       }
  83
  84     case GeomAbs_Plane:
  85     default:
  86       theErrFactorU = theErrFactorV = 1.;
  87     }
  88   }
  89
  90   void AdjustCellsCounts (const Handle(Adaptor3d_HSurface)& theFace,
  91                           const Standard_Integer            theNbVertices,
  92                           Standard_Integer&                 theCellsCountU,
  93                           Standard_Integer&                 theCellsCountV)
  94   {
  95     const GeomAbs_SurfaceType aType = theFace->GetType ();
  96     if (aType == GeomAbs_OtherSurface)
  97     {
  98       // fallback to the default behavior
  99       theCellsCountU = theCellsCountV = -1;
 100       return;
 101     }
 102
 103     Standard_Real aSqNbVert = theNbVertices;
 104     if (aType == GeomAbs_Plane)
 105     {
 106       theCellsCountU = theCellsCountV = (Standard_Integer)Ceiling (Pow (2, Log10 (aSqNbVert)));
 107     }
 108     else if (aType == GeomAbs_Cylinder || aType == GeomAbs_Cone)
 109     {
 110       theCellsCountV = (Standard_Integer)Ceiling (Pow (2, Log10 (aSqNbVert)));
 111     }
 112     else if (aType == GeomAbs_SurfaceOfExtrusion || aType == GeomAbs_SurfaceOfRevolution)
 113     {
 114       Handle (Adaptor3d_HCurve) aCurve = theFace->BasisCurve ();
 115       if (aCurve->GetType () == GeomAbs_Line ||
 116          (aCurve->GetType () == GeomAbs_BSplineCurve && aCurve->Degree () < 2))
 117       {
 118         // planar, cylindrical, conical cases
 119         if (aType == GeomAbs_SurfaceOfExtrusion)
 120           theCellsCountU = (Standard_Integer)Ceiling (Pow (2, Log10 (aSqNbVert)));
 121         else
 122           theCellsCountV = (Standard_Integer)Ceiling (Pow (2, Log10 (aSqNbVert)));
 123       }
 124       if (aType == GeomAbs_SurfaceOfExtrusion)
 125       {
 126         // V is always a line
 127         theCellsCountV = (Standard_Integer)Ceiling (Pow (2, Log10 (aSqNbVert)));
 128       }
 129     }
 130     else if (aType == GeomAbs_BezierSurface || aType == GeomAbs_BSplineSurface)
 131     {
 132       if (theFace->UDegree () < 2)
 133       {
 134         theCellsCountU = (Standard_Integer)Ceiling (Pow (2, Log10 (aSqNbVert)));
 135       }
 136       if (theFace->VDegree () < 2)
 137       {
 138         theCellsCountV = (Standard_Integer)Ceiling (Pow (2, Log10 (aSqNbVert)));
 139       }
 140     }
 141
 142     theCellsCountU = Max (theCellsCountU, 2);
 143     theCellsCountV = Max (theCellsCountV, 2);
 144   }
 145 }
 146
 147 //=======================================================================
 148 //function : Constructor
 149 //purpose  :
 150 //=======================================================================
 151 BRepMesh_GeomTool::BRepMesh_GeomTool(
 152   const BRepAdaptor_Curve& theCurve,
 153   const Standard_Real      theFirstParam,
 154   const Standard_Real      theLastParam,
 155   const Standard_Real      theLinDeflection,
 156   const Standard_Real      theAngDeflection,
 157   const Standard_Integer   theMinPointsNb,
 158   const Standard_Real      theMinSize)
 159   : myEdge(&theCurve.Edge()),
 160     myIsoType(GeomAbs_NoneIso)
 161 {
 162   myDiscretTool.Initialize(theCurve, theFirstParam, theLastParam,
 163     theAngDeflection, theLinDeflection, theMinPointsNb,
 164     Precision::PConfusion(), theMinSize);
 165 }
 166
 167 //=======================================================================
 168 //function : Constructor
 169 //purpose  :
 170 //=======================================================================
 171 BRepMesh_GeomTool::BRepMesh_GeomTool(
 172   const Handle(BRepAdaptor_HSurface)& theSurface,
 173   const GeomAbs_IsoType               theIsoType,
 174   const Standard_Real                 theParamIso,
 175   const Standard_Real                 theFirstParam,
 176   const Standard_Real                 theLastParam,
 177   const Standard_Real                 theLinDeflection,
 178   const Standard_Real                 theAngDeflection,
 179   const Standard_Integer              theMinPointsNb,
 180   const Standard_Real                 theMinSize)
 181   : myEdge(NULL),
 182     myIsoType(theIsoType)
 183 {
 184   Adaptor3d_IsoCurve aIso(theSurface, theIsoType, theParamIso,
 185     theFirstParam, theLastParam);
 186
 187   myDiscretTool.Initialize(aIso, theFirstParam, theLastParam,
 188     theAngDeflection, theLinDeflection, theMinPointsNb,
 189     Precision::PConfusion(), theMinSize);
 190 }
 191
 192 //=======================================================================
 193 //function : Value
 194 //purpose  :
 195 //=======================================================================
 196 Standard_Boolean BRepMesh_GeomTool::Value(
 197   const Standard_Integer              theIndex,
 198   const Handle(BRepAdaptor_HSurface)& theSurface,
 199   Standard_Real&                      theParam,
 200   gp_Pnt&                             thePoint,
 201   gp_Pnt2d&                           theUV) const
 202 {
 203   if (theIndex < 1 || theIndex > NbPoints())
 204     return Standard_False;
 205
 206   if (myEdge == NULL)
 207     return Standard_False;
 208
 209   thePoint = myDiscretTool.Value(theIndex);
 210   theParam = myDiscretTool.Parameter(theIndex);
 211
 212   const TopoDS_Face& aFace = ((BRepAdaptor_Surface*)&(theSurface->Surface()))->Face();
 213
 214   Standard_Real aFirst, aLast;
 215   Handle(Geom2d_Curve) aCurve =
 216     BRep_Tool::CurveOnSurface(*myEdge, aFace, aFirst, aLast);
 217
 218   aCurve->D0(theParam, theUV);
 219
 220   return Standard_True;
 221 }
 222
 223 //=======================================================================
 224 //function : Value
 225 //purpose  :
 226 //=======================================================================
 227 Standard_Boolean BRepMesh_GeomTool::Value(
 228   const Standard_Integer theIndex,
 229   const Standard_Real    theIsoParam,
 230   Standard_Real&         theParam,
 231   gp_Pnt&                thePoint,
 232   gp_Pnt2d&              theUV) const
 233 {
 234   if (theIndex < 1 || theIndex > NbPoints())
 235     return Standard_False;
 236
 237   thePoint = myDiscretTool.Value(theIndex);
 238   theParam = myDiscretTool.Parameter(theIndex);
 239
 240   if (myIsoType == GeomAbs_IsoU)
 241     theUV.SetCoord(theIsoParam, theParam);
 242   else
 243     theUV.SetCoord(theParam, theIsoParam);
 244
 245   return Standard_True;
 246 }
 247
 248 //=======================================================================
 249 //function : Normal
 250 //purpose  :
 251 //=======================================================================
 252 Standard_Boolean BRepMesh_GeomTool::Normal(
 253   const Handle(BRepAdaptor_HSurface)& theSurface,
 254   const Standard_Real                 theParamU,
 255   const Standard_Real                 theParamV,
 256   gp_Pnt&                             thePoint,
 257   gp_Dir&                             theNormal)
 258 {
 259   Standard_Boolean isOK = Standard_True;
 260   gp_Vec aD1U, aD1V;
 261
 262   theSurface->D1(theParamU, theParamV, thePoint, aD1U, aD1V);
 263
 264   CSLib_DerivativeStatus aStatus;
 265   CSLib::Normal(aD1U, aD1V, Precision::Angular(), aStatus, theNormal);
 266   if (aStatus != CSLib_Done)
 267   {
 268     gp_Vec aD2U,aD2V,aD2UV;
 269     theSurface->D2(theParamU, theParamV, thePoint, aD1U, aD1V, aD2U, aD2V, aD2UV);
 270     CSLib_NormalStatus aNormalStatus;
 271     CSLib::Normal(aD1U, aD1V, aD2U, aD2V, aD2UV, Precision::Angular(),
 272       isOK, aNormalStatus, theNormal);
 273   }
 274
 275   if (!isOK)
 276     return Standard_False;
 277
 278   const TopoDS_Face& aFace = ((BRepAdaptor_Surface*)&(theSurface->Surface()))->Face();
 279   TopAbs_Orientation aOri = aFace.Orientation();
 280   if (aOri == TopAbs_REVERSED)
 281     theNormal.Reverse();
 282
 283   return Standard_True;
 284 }
 285
 286 //=============================================================================
 287 //function : IntLinLin
 288 //purpose  :
 289 //=============================================================================
 290 BRepMesh_GeomTool::IntFlag BRepMesh_GeomTool::IntLinLin(
 291   const gp_XY&  theStartPnt1,
 292   const gp_XY&  theEndPnt1,
 293   const gp_XY&  theStartPnt2,
 294   const gp_XY&  theEndPnt2,
 295   gp_XY&        theIntPnt,
 296   Standard_Real (&theParamOnSegment)[2])
 297 {
 298   gp_XY aVec1    = theEndPnt1   - theStartPnt1;
 299   gp_XY aVec2    = theEndPnt2   - theStartPnt2;
 300   gp_XY aVecO1O2 = theStartPnt2 - theStartPnt1;
 301
 302   Standard_Real aCrossD1D2 = aVec1    ^ aVec2;
 303   Standard_Real aCrossD1D3 = aVecO1O2 ^ aVec2;
 304
 305   const Standard_Real aPrec = gp::Resolution();
 306   // Are edgegs codirectional
 307   if ( Abs( aCrossD1D2 ) < aPrec )
 308   {
 309     // Just a parallel case?
 310     if( Abs( aCrossD1D3 ) < aPrec )
 311       return BRepMesh_GeomTool::Same;
 312     else
 313       return BRepMesh_GeomTool::NoIntersection;
 314   }
 315
 316   theParamOnSegment[0] = aCrossD1D3 / aCrossD1D2;
 317   theIntPnt = theStartPnt1 + theParamOnSegment[0] * aVec1;
 318
 319   Standard_Real aCrossD2D3 = aVecO1O2 ^ aVec1;
 320   theParamOnSegment[1] = aCrossD2D3 / aCrossD1D2;
 321
 322   return BRepMesh_GeomTool::Cross;
 323 }
 324
 325 //=============================================================================
 326 //function : IntSegSeg
 327 //purpose  :
 328 //=============================================================================
 329 BRepMesh_GeomTool::IntFlag BRepMesh_GeomTool::IntSegSeg(
 330   const gp_XY&           theStartPnt1,
 331   const gp_XY&           theEndPnt1,
 332   const gp_XY&           theStartPnt2,
 333   const gp_XY&           theEndPnt2,
 334   const Standard_Boolean isConsiderEndPointTouch,
 335   const Standard_Boolean isConsiderPointOnSegment,
 336   gp_Pnt2d&              theIntPnt)
 337 {
 338   Standard_Integer aPointHash[] = {
 339     classifyPoint(theStartPnt1, theEndPnt1, theStartPnt2),
 340     classifyPoint(theStartPnt1, theEndPnt1, theEndPnt2  ),
 341     classifyPoint(theStartPnt2, theEndPnt2, theStartPnt1),
 342     classifyPoint(theStartPnt2, theEndPnt2, theEndPnt1  )
 343   };
 344
 345   // Consider case when edges have shared vertex
 346   if ( aPointHash[0] < 0 || aPointHash[1] < 0 )
 347   {
 348     if ( isConsiderEndPointTouch )
 349       return BRepMesh_GeomTool::EndPointTouch;
 350
 351     return BRepMesh_GeomTool::NoIntersection;
 352   }
 353
 354   Standard_Integer aPosHash =
 355     aPointHash[0] + aPointHash[1] + aPointHash[2] + aPointHash[3];
 356
 357   /*=========================================*/
 358   /*  1) hash code == 1:
 359
 360                     0+
 361                     /
 362            0      1/         0
 363            +======+==========+
 364
 365       2) hash code == 2:
 366
 367            0    1        1   0
 368         a) +----+========+---+
 369
 370            0       1   1     0
 371         b) +-------+===+=====+
 372
 373                                              */
 374   /*=========================================*/
 375   if ( aPosHash == 1 )
 376   {
 377     if (isConsiderPointOnSegment)
 378     {
 379       if (aPointHash[0] == 1)
 380         theIntPnt = theStartPnt1;
 381       else if (aPointHash[1] == 1)
 382         theIntPnt = theEndPnt1;
 383       else if (aPointHash[2] == 1)
 384         theIntPnt = theStartPnt2;
 385       else
 386         theIntPnt = theEndPnt2;
 387
 388       return BRepMesh_GeomTool::PointOnSegment;
 389     }
 390
 391     return BRepMesh_GeomTool::NoIntersection;
 392   }
 393   else if ( aPosHash == 2 )
 394     return BRepMesh_GeomTool::Glued;
 395
 396   Standard_Real aParam[2];
 397   IntFlag aIntFlag = IntLinLin(theStartPnt1, theEndPnt1,
 398     theStartPnt2, theEndPnt2, theIntPnt.ChangeCoord(), aParam);
 399
 400   if (aIntFlag == BRepMesh_GeomTool::NoIntersection)
 401     return BRepMesh_GeomTool::NoIntersection;
 402
 403   if (aIntFlag == BRepMesh_GeomTool::Same)
 404   {
 405     if ( aPosHash < -2 )
 406       return BRepMesh_GeomTool::Same;
 407     else if ( aPosHash == -1 )
 408       return BRepMesh_GeomTool::Glued;
 409
 410     return BRepMesh_GeomTool::NoIntersection;
 411   }
 412
 413   // Cross
 414   // Intersection is out of segments ranges
 415   const Standard_Real aPrec    = Precision::PConfusion();
 416   const Standard_Real aEndPrec = 1 - aPrec;
 417   for (Standard_Integer i = 0; i < 2; ++i)
 418   {
 419     if(aParam[i] < aPrec || aParam[i] > aEndPrec )
 420       return BRepMesh_GeomTool::NoIntersection;
 421   }
 422
 423   return BRepMesh_GeomTool::Cross;
 424 }
 425
 426 //=============================================================================
 427 //function : CellsCount
 428 //purpose  :
 429 //=============================================================================
 430 std::pair<Standard_Integer, Standard_Integer> BRepMesh_GeomTool::CellsCount (
 431   const Handle (Adaptor3d_HSurface)&   theSurface,
 432   const Standard_Integer               theVerticesNb,
 433   const Standard_Real                  theDeflection,
 434   const BRepMesh_DefaultRangeSplitter* theRangeSplitter)
 435 {
 436   if (theRangeSplitter == NULL)
 437     return std::pair<Standard_Integer, Standard_Integer>(-1, -1);
 438
 439   const GeomAbs_SurfaceType aType = theSurface->GetType ();
 440
 441   Standard_Real anErrFactorU, anErrFactorV;
 442   ComputeErrFactors(theDeflection, theSurface, anErrFactorU, anErrFactorV);
 443
 444   const std::pair<Standard_Real, Standard_Real>& aRangeU = theRangeSplitter->GetRangeU();
 445   const std::pair<Standard_Real, Standard_Real>& aRangeV = theRangeSplitter->GetRangeV();
 446   const std::pair<Standard_Real, Standard_Real>& aDelta  = theRangeSplitter->GetDelta ();
 447
 448   Standard_Integer aCellsCountU, aCellsCountV;
 449   if (aType == GeomAbs_Torus)
 450   {
 451     aCellsCountU = (Standard_Integer)Ceiling(Pow(2, Log10(
 452       (aRangeU.second - aRangeU.first) / aDelta.first)));
 453     aCellsCountV = (Standard_Integer)Ceiling(Pow(2, Log10(
 454       (aRangeV.second - aRangeV.first) / aDelta.second)));
 455   }
 456   else if (aType == GeomAbs_Cylinder)
 457   {
 458     aCellsCountU = (Standard_Integer)Ceiling(Pow(2, Log10(
 459       (aRangeU.second - aRangeU.first) / aDelta.first /
 460       (aRangeV.second - aRangeV.first))));
 461     aCellsCountV = (Standard_Integer)Ceiling(Pow(2, Log10(
 462       (aRangeV.second - aRangeV.first) / anErrFactorV)));
 463   }
 464   else
 465   {
 466     aCellsCountU = (Standard_Integer)Ceiling(Pow(2, Log10(
 467       (aRangeU.second - aRangeU.first) / aDelta.first  / anErrFactorU)));
 468     aCellsCountV = (Standard_Integer)Ceiling(Pow(2, Log10(
 469       (aRangeV.second - aRangeV.first) / aDelta.second / anErrFactorV)));
 470   }
 471
 472   AdjustCellsCounts(theSurface, theVerticesNb, aCellsCountU, aCellsCountV);
 473   return std::pair<Standard_Integer, Standard_Integer>(aCellsCountU, aCellsCountV);
 474 }
 475
 476 //=============================================================================
 477 //function : classifyPoint
 478 //purpose  :
 479 //=============================================================================
 480 Standard_Integer BRepMesh_GeomTool::classifyPoint(
 481   const gp_XY& thePoint1,
 482   const gp_XY& thePoint2,
 483   const gp_XY& thePointToCheck)
 484 {
 485   gp_XY aP1 = thePoint2       - thePoint1;
 486   gp_XY aP2 = thePointToCheck - thePoint1;
 487
 488   const Standard_Real aPrec   = Precision::PConfusion();
 489   const Standard_Real aSqPrec = aPrec * aPrec;
 490   Standard_Real aDist = Abs(aP1 ^ aP2);
 491   if (aDist > aPrec)
 492   {
 493     aDist = (aDist * aDist) / aP1.SquareModulus();
 494     if (aDist > aSqPrec)
 495       return 0; //out
 496   }
 497
 498   gp_XY aMult = aP1.Multiplied(aP2);
 499   if ( aMult.X() < 0.0 || aMult.Y() < 0.0 )
 500     return 0; //out
 501
 502   if (aP1.SquareModulus() < aP2.SquareModulus())
 503     return 0; //out
 504
 505   if (thePointToCheck.IsEqual(thePoint1, aPrec) ||
 506       thePointToCheck.IsEqual(thePoint2, aPrec))
 507   {
 508     return -1; //coinsides with an end point
 509   }
 510
 511   return 1;
 512 }