src/IntTools/IntTools_WLineTool.cxx

   1 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   2 //
   3 // This file is part of Open CASCADE Technology software library.
   4 //
   5 // This library is free software; you can redistribute it and/or modify it under
   6 // the terms of the GNU Lesser General Public License version 2.1 as published
   7 // by the Free Software Foundation, with special exception defined in the file
   8 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
   9 // distribution for complete text of the license and disclaimer of any warranty.
  10 //
  11 // Alternatively, this file may be used under the terms of Open CASCADE
  12 // commercial license or contractual agreement.
  13
  14 #include <IntTools_WLineTool.hxx>
  15
  16 #include <BRep_Tool.hxx>
  17 #include <Extrema_ExtCC.hxx>
  18 #include <Geom2dAPI_InterCurveCurve.hxx>
  19 #include <Geom2d_Circle.hxx>
  20 #include <Geom2d_Line.hxx>
  21 #include <Geom2d_TrimmedCurve.hxx>
  22 #include <GeomAPI_ProjectPointOnSurf.hxx>
  23 #include <GeomAdaptor_Curve.hxx>
  24 #include <GeomAdaptor_HSurface.hxx>
  25 #include <GeomAdaptor_Surface.hxx>
  26 #include <GeomInt.hxx>
  27 #include <GeomInt_LineConstructor.hxx>
  28 #include <Geom_Circle.hxx>
  29 #include <Geom_Surface.hxx>
  30 #include <gp_Circ.hxx>
  31 #include <IntTools_Context.hxx>
  32 #include <TColStd_Array1OfListOfInteger.hxx>
  33 #include <TColStd_SequenceOfReal.hxx>
  34 #include <TColgp_SequenceOfPnt2d.hxx>
  35
  36 /////////////////////// NotUseSurfacesForApprox /////////////////////////
  37
  38 // The block is dedicated to determine whether WLine [ifprm, ilprm]
  39 // crosses the degenerated zone on each given surface or not.
  40 // If Yes -> We will not use info about surfaces during approximation
  41 // because inside degenerated zone of the surface the approx. algo.
  42 // uses wrong values of normal, etc., and resulting curve will have
  43 // oscillations that we would not like to have.
  44
  45 //=======================================================================
  46 //function : IsDegeneratedZone
  47 //purpose  : static subfunction in IsDegeneratedZone
  48 //=======================================================================
  49 static
  50 Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
  51                                    const Handle(Geom_Surface)& aS,
  52                                    const Standard_Integer iDir)
  53 {
  54   Standard_Boolean bFlag=Standard_True;
  55   Standard_Real US1, US2, VS1, VS2, dY, dX, d1, d2, dD;
  56   Standard_Real aXm, aYm, aXb, aYb, aXe, aYe;
  57   aS->Bounds(US1, US2, VS1, VS2);
  58
  59   gp_Pnt aPm, aPb, aPe;
  60
  61   aXm=aP2d.X();
  62   aYm=aP2d.Y();
  63
  64   aS->D0(aXm, aYm, aPm);
  65
  66   dX=1.e-5;
  67   dY=1.e-5;
  68   dD=1.e-12;
  69
  70   if (iDir==1) {
  71     aXb=aXm;
  72     aXe=aXm;
  73     aYb=aYm-dY;
  74     if (aYb < VS1) {
  75       aYb=VS1;
  76     }
  77     aYe=aYm+dY;
  78     if (aYe > VS2) {
  79       aYe=VS2;
  80     }
  81     aS->D0(aXb, aYb, aPb);
  82     aS->D0(aXe, aYe, aPe);
  83
  84     d1=aPm.Distance(aPb);
  85     d2=aPm.Distance(aPe);
  86     if (d1 < dD && d2 < dD) {
  87       return bFlag;
  88     }
  89     return !bFlag;
  90   }
  91   //
  92   else if (iDir==2) {
  93     aYb=aYm;
  94     aYe=aYm;
  95     aXb=aXm-dX;
  96     if (aXb < US1) {
  97       aXb=US1;
  98     }
  99     aXe=aXm+dX;
 100     if (aXe > US2) {
 101       aXe=US2;
 102     }
 103     aS->D0(aXb, aYb, aPb);
 104     aS->D0(aXe, aYe, aPe);
 105
 106     d1=aPm.Distance(aPb);
 107     d2=aPm.Distance(aPe);
 108     if (d1 < dD && d2 < dD) {
 109       return bFlag;
 110     }
 111     return !bFlag;
 112   }
 113   return !bFlag;
 114 }
 115
 116 //=======================================================================
 117 //function : IsPointInDegeneratedZone
 118 //purpose  : static subfunction in NotUseSurfacesForApprox
 119 //=======================================================================
 120 static
 121 Standard_Boolean IsPointInDegeneratedZone(const IntSurf_PntOn2S& aP2S,
 122                                           const TopoDS_Face& aF1,
 123                                           const TopoDS_Face& aF2)
 124
 125 {
 126   Standard_Boolean bFlag=Standard_True;
 127   Standard_Real US11, US12, VS11, VS12, US21, US22, VS21, VS22;
 128   Standard_Real U1, V1, U2, V2, aDelta, aD;
 129   gp_Pnt2d aP2d;
 130
 131   Handle(Geom_Surface)aS1 = BRep_Tool::Surface(aF1);
 132   aS1->Bounds(US11, US12, VS11, VS12);
 133   GeomAdaptor_Surface aGAS1(aS1);
 134
 135   Handle(Geom_Surface)aS2 = BRep_Tool::Surface(aF2);
 136   aS1->Bounds(US21, US22, VS21, VS22);
 137   GeomAdaptor_Surface aGAS2(aS2);
 138   //
 139   //const gp_Pnt& aP=aP2S.Value();
 140   aP2S.Parameters(U1, V1, U2, V2);
 141   //
 142   aDelta=1.e-7;
 143   // Check on Surf 1
 144   aD=aGAS1.UResolution(aDelta);
 145   aP2d.SetCoord(U1, V1);
 146   if (fabs(U1-US11) < aD) {
 147     bFlag=IsDegeneratedZone(aP2d, aS1, 1);
 148     if (bFlag) {
 149       return bFlag;
 150     }
 151   }
 152   if (fabs(U1-US12) < aD) {
 153     bFlag=IsDegeneratedZone(aP2d, aS1, 1);
 154     if (bFlag) {
 155       return bFlag;
 156     }
 157   }
 158   aD=aGAS1.VResolution(aDelta);
 159   if (fabs(V1-VS11) < aDelta) {
 160     bFlag=IsDegeneratedZone(aP2d, aS1, 2);
 161     if (bFlag) {
 162       return bFlag;
 163     }
 164   }
 165   if (fabs(V1-VS12) < aDelta) {
 166     bFlag=IsDegeneratedZone(aP2d, aS1, 2);
 167     if (bFlag) {
 168       return bFlag;
 169     }
 170   }
 171   // Check on Surf 2
 172   aD=aGAS2.UResolution(aDelta);
 173   aP2d.SetCoord(U2, V2);
 174   if (fabs(U2-US21) < aDelta) {
 175     bFlag=IsDegeneratedZone(aP2d, aS2, 1);
 176     if (bFlag) {
 177       return bFlag;
 178     }
 179   }
 180   if (fabs(U2-US22) < aDelta) {
 181     bFlag=IsDegeneratedZone(aP2d, aS2, 1);
 182     if (bFlag) {
 183       return bFlag;
 184     }
 185   }
 186   aD=aGAS2.VResolution(aDelta);
 187   if (fabs(V2-VS21) < aDelta) {
 188     bFlag=IsDegeneratedZone(aP2d, aS2, 2);
 189     if (bFlag) {
 190       return bFlag;
 191     }
 192   }
 193   if (fabs(V2-VS22) < aDelta) {
 194     bFlag=IsDegeneratedZone(aP2d, aS2, 2);
 195     if (bFlag) {
 196       return bFlag;
 197     }
 198   }
 199   return !bFlag;
 200 }
 201
 202 //=======================================================================
 203 //function :  NotUseSurfacesForApprox
 204 //purpose  :
 205 //=======================================================================
 206 Standard_Boolean IntTools_WLineTool::NotUseSurfacesForApprox(const TopoDS_Face& aF1,
 207                                                              const TopoDS_Face& aF2,
 208                                                              const Handle(IntPatch_WLine)& WL,
 209                                                              const Standard_Integer ifprm,
 210                                                              const Standard_Integer ilprm)
 211 {
 212   Standard_Boolean bPInDZ;
 213
 214   Handle(IntSurf_LineOn2S) aLineOn2S=WL->Curve();
 215
 216   const IntSurf_PntOn2S& aP2Sfprm=aLineOn2S->Value(ifprm);
 217   bPInDZ=IsPointInDegeneratedZone(aP2Sfprm, aF1, aF2);
 218   if (bPInDZ) {
 219     return bPInDZ;
 220   }
 221
 222   const IntSurf_PntOn2S& aP2Slprm=aLineOn2S->Value(ilprm);
 223   bPInDZ=IsPointInDegeneratedZone(aP2Slprm, aF1, aF2);
 224
 225   return bPInDZ;
 226 }
 227
 228 /////////////////////// end of NotUseSurfacesForApprox //////////////////
 229
 230 /////////////////////// DecompositionOfWLine ////////////////////////////
 231
 232 //=======================================================================
 233 //function : CheckTangentZonesExist
 234 //purpose  : static subfunction in ComputeTangentZones
 235 //=======================================================================
 236 static
 237 Standard_Boolean CheckTangentZonesExist(const Handle(GeomAdaptor_HSurface)& theSurface1,
 238                                         const Handle(GeomAdaptor_HSurface)& theSurface2)
 239 {
 240   if ( ( theSurface1->GetType() != GeomAbs_Torus ) ||
 241        ( theSurface2->GetType() != GeomAbs_Torus ) )
 242     return Standard_False;
 243
 244   gp_Torus aTor1 = theSurface1->Torus();
 245   gp_Torus aTor2 = theSurface2->Torus();
 246
 247   if ( aTor1.Location().Distance( aTor2.Location() ) > Precision::Confusion() )
 248     return Standard_False;
 249
 250   if ( ( fabs( aTor1.MajorRadius() - aTor2.MajorRadius() ) > Precision::Confusion() ) ||
 251        ( fabs( aTor1.MinorRadius() - aTor2.MinorRadius() ) > Precision::Confusion() ) )
 252     return Standard_False;
 253
 254   if ( ( aTor1.MajorRadius() < aTor1.MinorRadius() ) ||
 255        ( aTor2.MajorRadius() < aTor2.MinorRadius() ) )
 256     return Standard_False;
 257
 258   return Standard_True;
 259 }
 260
 261
 262 //=======================================================================
 263 //function : ComputeTangentZones
 264 //purpose  : static subfunction in DecompositionOfWLine
 265 //=======================================================================
 266 static
 267 Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSurface1,
 268                                      const Handle(GeomAdaptor_HSurface)&  theSurface2,
 269                                      const TopoDS_Face&                   theFace1,
 270                                      const TopoDS_Face&                   theFace2,
 271                                      Handle(TColgp_HArray1OfPnt2d)&       theResultOnS1,
 272                                      Handle(TColgp_HArray1OfPnt2d)&       theResultOnS2,
 273                                      Handle(TColStd_HArray1OfReal)&       theResultRadius,
 274                                      const Handle(IntTools_Context)& aContext)
 275 {
 276   Standard_Integer aResult = 0;
 277   if ( !CheckTangentZonesExist( theSurface1, theSurface2 ) )
 278     return aResult;
 279
 280
 281   TColgp_SequenceOfPnt2d aSeqResultS1, aSeqResultS2;
 282   TColStd_SequenceOfReal aSeqResultRad;
 283
 284   gp_Torus aTor1 = theSurface1->Torus();
 285   gp_Torus aTor2 = theSurface2->Torus();
 286
 287   gp_Ax2 anax1( aTor1.Location(), aTor1.Axis().Direction() );
 288   gp_Ax2 anax2( aTor2.Location(), aTor2.Axis().Direction() );
 289   Standard_Integer j = 0;
 290
 291   for ( j = 0; j < 2; j++ ) {
 292     Standard_Real aCoef = ( j == 0 ) ? -1 : 1;
 293     Standard_Real aRadius1 = fabs(aTor1.MajorRadius() + aCoef * aTor1.MinorRadius());
 294     Standard_Real aRadius2 = fabs(aTor2.MajorRadius() + aCoef * aTor2.MinorRadius());
 295
 296     gp_Circ aCircle1( anax1, aRadius1 );
 297     gp_Circ aCircle2( anax2, aRadius2 );
 298
 299     // roughly compute radius of tangent zone for perpendicular case
 300     Standard_Real aCriteria = Precision::Confusion() * 0.5;
 301
 302     Standard_Real aT1 = aCriteria;
 303     Standard_Real aT2 = aCriteria;
 304     if ( j == 0 ) {
 305       // internal tangency
 306       Standard_Real aR = ( aRadius1 > aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
 307       //aT1 = aCriteria * aCriteria + aR * aR - ( aR - aCriteria ) * ( aR - aCriteria );
 308       aT1 = 2. * aR * aCriteria;
 309       aT2 = aT1;
 310     }
 311     else {
 312       // external tangency
 313       Standard_Real aRb = ( aRadius1 > aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
 314       Standard_Real aRm = ( aRadius1 < aTor2.MinorRadius() ) ? aRadius1 : aTor2.MinorRadius();
 315       Standard_Real aDelta = aRb - aCriteria;
 316       aDelta *= aDelta;
 317       aDelta -= aRm * aRm;
 318       aDelta /= 2. * (aRb - aRm);
 319       aDelta -= 0.5 * (aRb - aRm);
 320
 321       aT1 = 2. * aRm * (aRm - aDelta);
 322       aT2 = aT1;
 323     }
 324     aCriteria = ( aT1 > aT2) ? aT1 : aT2;
 325     if ( aCriteria > 0 )
 326       aCriteria = sqrt( aCriteria );
 327
 328     if ( aCriteria > 0.5 * aTor1.MinorRadius() ) {
 329       // too big zone -> drop to minimum
 330       aCriteria = Precision::Confusion();
 331     }
 332
 333     GeomAdaptor_Curve aC1( new Geom_Circle(aCircle1) );
 334     GeomAdaptor_Curve aC2( new Geom_Circle(aCircle2) );
 335     Extrema_ExtCC anExtrema(aC1, aC2, 0, 2. * M_PI, 0, 2. * M_PI,
 336                             Precision::PConfusion(), Precision::PConfusion());
 337
 338     if ( anExtrema.IsDone() ) {
 339
 340       Standard_Integer i = 0;
 341       for ( i = 1; i <= anExtrema.NbExt(); i++ ) {
 342         if ( anExtrema.SquareDistance(i) > aCriteria * aCriteria )
 343           continue;
 344
 345         Extrema_POnCurv P1, P2;
 346         anExtrema.Points( i, P1, P2 );
 347
 348         Standard_Boolean bFoundResult = Standard_True;
 349         gp_Pnt2d pr1, pr2;
 350
 351         Standard_Integer surfit = 0;
 352         for ( surfit = 0; surfit < 2; surfit++ ) {
 353           GeomAPI_ProjectPointOnSurf& aProjector =
 354             (surfit == 0) ? aContext->ProjPS(theFace1) : aContext->ProjPS(theFace2);
 355
 356           gp_Pnt aP3d = (surfit == 0) ? P1.Value() : P2.Value();
 357           aProjector.Perform(aP3d);
 358
 359           if(!aProjector.IsDone())
 360             bFoundResult = Standard_False;
 361           else {
 362             if(aProjector.LowerDistance() > aCriteria) {
 363               bFoundResult = Standard_False;
 364             }
 365             else {
 366               Standard_Real foundU = 0, foundV = 0;
 367               aProjector.LowerDistanceParameters(foundU, foundV);
 368               if ( surfit == 0 )
 369                 pr1 = gp_Pnt2d( foundU, foundV );
 370               else
 371                 pr2 = gp_Pnt2d( foundU, foundV );
 372             }
 373           }
 374         }
 375         if ( bFoundResult ) {
 376           aSeqResultS1.Append( pr1 );
 377           aSeqResultS2.Append( pr2 );
 378           aSeqResultRad.Append( aCriteria );
 379
 380           // torus is u and v periodic
 381           const Standard_Real twoPI = M_PI + M_PI;
 382           Standard_Real arr1tmp[2] = {pr1.X(), pr1.Y()};
 383           Standard_Real arr2tmp[2] = {pr2.X(), pr2.Y()};
 384
 385           // iteration on period bounds
 386           for ( Standard_Integer k1 = 0; k1 < 2; k1++ ) {
 387             Standard_Real aBound = ( k1 == 0 ) ? 0 : twoPI;
 388             Standard_Real aShift = ( k1 == 0 ) ? twoPI : -twoPI;
 389
 390             // iteration on surfaces
 391             for ( Standard_Integer k2 = 0; k2 < 2; k2++ ) {
 392               Standard_Real* arr1 = ( k2 == 0 ) ? arr1tmp : arr2tmp;
 393               Standard_Real* arr2 = ( k2 != 0 ) ? arr1tmp : arr2tmp;
 394               TColgp_SequenceOfPnt2d& aSeqS1 = ( k2 == 0 ) ? aSeqResultS1 : aSeqResultS2;
 395               TColgp_SequenceOfPnt2d& aSeqS2 = ( k2 != 0 ) ? aSeqResultS1 : aSeqResultS2;
 396
 397               if (fabs(arr1[0] - aBound) < Precision::PConfusion()) {
 398                 aSeqS1.Append( gp_Pnt2d( arr1[0] + aShift, arr1[1] ) );
 399                 aSeqS2.Append( gp_Pnt2d( arr2[0], arr2[1] ) );
 400                 aSeqResultRad.Append( aCriteria );
 401               }
 402               if (fabs(arr1[1] - aBound) < Precision::PConfusion()) {
 403                 aSeqS1.Append( gp_Pnt2d( arr1[0], arr1[1] + aShift) );
 404                 aSeqS2.Append( gp_Pnt2d( arr2[0], arr2[1] ) );
 405                 aSeqResultRad.Append( aCriteria );
 406               }
 407             }
 408           } //
 409         }
 410       }
 411     }
 412   }
 413   aResult = aSeqResultRad.Length();
 414
 415   if ( aResult > 0 ) {
 416     theResultOnS1 = new TColgp_HArray1OfPnt2d( 1, aResult );
 417     theResultOnS2 = new TColgp_HArray1OfPnt2d( 1, aResult );
 418     theResultRadius = new TColStd_HArray1OfReal( 1, aResult );
 419
 420     for ( Standard_Integer i = 1 ; i <= aResult; i++ ) {
 421       theResultOnS1->SetValue( i, aSeqResultS1.Value(i) );
 422       theResultOnS2->SetValue( i, aSeqResultS2.Value(i) );
 423       theResultRadius->SetValue( i, aSeqResultRad.Value(i) );
 424     }
 425   }
 426   return aResult;
 427 }
 428
 429 //=======================================================================
 430 //function : IsPointOnBoundary
 431 //purpose  : static subfunction in DecompositionOfWLine
 432 //=======================================================================
 433 static
 434 Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
 435                                    const Standard_Real theFirstBoundary,
 436                                    const Standard_Real theSecondBoundary,
 437                                    const Standard_Real theResolution,
 438                                    Standard_Boolean&   IsOnFirstBoundary)
 439 {
 440   Standard_Boolean bRet;
 441   Standard_Integer i;
 442   Standard_Real adist;
 443   //
 444   bRet=Standard_False;
 445   for(i = 0; i < 2; ++i) {
 446     IsOnFirstBoundary = (i == 0);
 447     if (IsOnFirstBoundary) {
 448       adist = fabs(theParameter - theFirstBoundary);
 449     }
 450     else {
 451       adist = fabs(theParameter - theSecondBoundary);
 452     }
 453     if(adist < theResolution) {
 454       return !bRet;
 455     }
 456   }
 457   return bRet;
 458 }
 459
 460 //=======================================================================
 461 //function : IsInsideTanZone
 462 //purpose  : Check if point is inside a radial tangent zone.
 463 //           static subfunction in DecompositionOfWLine and FindPoint
 464 //=======================================================================
 465 static
 466 Standard_Boolean IsInsideTanZone(const gp_Pnt2d&     thePoint,
 467                                  const gp_Pnt2d&     theTanZoneCenter,
 468                                  const Standard_Real theZoneRadius,
 469                                  Handle(GeomAdaptor_HSurface) theGASurface)
 470 {
 471   Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
 472   Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
 473   Standard_Real aRadiusSQR = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
 474   aRadiusSQR *= aRadiusSQR;
 475   if ( thePoint.SquareDistance( theTanZoneCenter ) <= aRadiusSQR )
 476     return Standard_True;
 477
 478   return Standard_False;
 479 }
 480
 481 //=======================================================================
 482 //function : AdjustByNeighbour
 483 //purpose  : static subfunction in DecompositionOfWLine
 484 //=======================================================================
 485 static
 486 gp_Pnt2d AdjustByNeighbour(const gp_Pnt2d&     theaNeighbourPoint,
 487                            const gp_Pnt2d&     theOriginalPoint,
 488                            Handle(GeomAdaptor_HSurface) theGASurface)
 489 {
 490   gp_Pnt2d ap1 = theaNeighbourPoint;
 491   gp_Pnt2d ap2 = theOriginalPoint;
 492
 493   if ( theGASurface->IsUPeriodic() ) {
 494     Standard_Real aPeriod     = theGASurface->UPeriod();
 495     gp_Pnt2d aPTest = ap2;
 496     Standard_Real aSqDistMin = 1.e+100;
 497
 498     for ( Standard_Integer pIt = -1; pIt <= 1; pIt++) {
 499       aPTest.SetX( theOriginalPoint.X() + aPeriod * pIt );
 500       Standard_Real dd = ap1.SquareDistance( aPTest );
 501
 502       if ( dd < aSqDistMin ) {
 503         ap2 = aPTest;
 504         aSqDistMin = dd;
 505       }
 506     }
 507   }
 508   if ( theGASurface->IsVPeriodic() ) {
 509     Standard_Real aPeriod     = theGASurface->VPeriod();
 510     gp_Pnt2d aPTest = ap2;
 511     Standard_Real aSqDistMin = 1.e+100;
 512
 513     for ( Standard_Integer pIt = -1; pIt <= 1; pIt++) {
 514       aPTest.SetY( theOriginalPoint.Y() + aPeriod * pIt );
 515       Standard_Real dd = ap1.SquareDistance( aPTest );
 516
 517       if ( dd < aSqDistMin ) {
 518         ap2 = aPTest;
 519         aSqDistMin = dd;
 520       }
 521     }
 522   }
 523   return ap2;
 524 }
 525
 526 //=======================================================================
 527 //function : RefineVector
 528 //purpose  : static subfunction in FindPoint
 529 //=======================================================================
 530 static
 531 void RefineVector(gp_Vec2d& aV2D)
 532 {
 533   Standard_Integer k,m;
 534   Standard_Real aC[2], aEps, aR1, aR2, aNum;
 535   //
 536   aEps=RealEpsilon();
 537   aR1=1.-aEps;
 538   aR2=1.+aEps;
 539   //
 540   aV2D.Coord(aC[0], aC[1]);
 541   //
 542   for (k=0; k<2; ++k) {
 543     m=(k+1)%2;
 544     aNum=fabs(aC[k]);
 545     if (aNum>aR1 && aNum<aR2) {
 546       if (aC[k]<0.) {
 547         aC[k]=-1.;
 548       }
 549       else {
 550         aC[k]=1.;
 551       }
 552       aC[m]=0.;
 553       break;
 554     }
 555   }
 556   aV2D.SetCoord(aC[0], aC[1]);
 557 }
 558
 559 //=======================================================================
 560 //function : FindPoint
 561 //purpose  : static subfunction in DecompositionOfWLine
 562 //=======================================================================
 563 static
 564 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
 565                            const gp_Pnt2d&     theLastPoint,
 566                            const Standard_Real theUmin,
 567                            const Standard_Real theUmax,
 568                            const Standard_Real theVmin,
 569                            const Standard_Real theVmax,
 570                            gp_Pnt2d&           theNewPoint)
 571 {
 572   gp_Vec2d aVec(theFirstPoint, theLastPoint);
 573   Standard_Integer i = 0, j = 0;
 574
 575   for(i = 0; i < 4; i++) {
 576     gp_Vec2d anOtherVec;
 577     gp_Vec2d anOtherVecNormal;
 578     gp_Pnt2d aprojpoint = theLastPoint;
 579
 580     if((i % 2) == 0) {
 581       anOtherVec.SetX(0.);
 582       anOtherVec.SetY(1.);
 583       anOtherVecNormal.SetX(1.);
 584       anOtherVecNormal.SetY(0.);
 585
 586       if(i < 2)
 587         aprojpoint.SetX(theUmin);
 588       else
 589         aprojpoint.SetX(theUmax);
 590     }
 591     else {
 592       anOtherVec.SetX(1.);
 593       anOtherVec.SetY(0.);
 594       anOtherVecNormal.SetX(0.);
 595       anOtherVecNormal.SetY(1.);
 596
 597       if(i < 2)
 598         aprojpoint.SetY(theVmin);
 599       else
 600         aprojpoint.SetY(theVmax);
 601     }
 602     gp_Vec2d anormvec = aVec;
 603     anormvec.Normalize();
 604     RefineVector(anormvec);
 605     Standard_Real adot1 = anormvec.Dot(anOtherVecNormal);
 606
 607     if(fabs(adot1) < Precision::Angular())
 608       continue;
 609     Standard_Real adist = 0.;
 610     Standard_Boolean bIsOut = Standard_False;
 611
 612     if((i % 2) == 0) {
 613       adist = (i < 2) ? fabs(theLastPoint.X() - theUmin) : fabs(theLastPoint.X() - theUmax);
 614       bIsOut = (i < 2) ? (theLastPoint.X() < theUmin) : (theLastPoint.X() > theUmax);
 615     }
 616     else {
 617       adist = (i < 2) ? fabs(theLastPoint.Y() - theVmin) : fabs(theLastPoint.Y() - theVmax);
 618       bIsOut = (i < 2) ? (theLastPoint.Y() < theVmin) : (theLastPoint.Y() > theVmax);
 619     }
 620     Standard_Real anoffset = adist * anOtherVec.Dot(anormvec) / adot1;
 621
 622     for(j = 0; j < 2; j++) {
 623       anoffset = (j == 0) ? anoffset : -anoffset;
 624       gp_Pnt2d acurpoint(aprojpoint.XY() + (anOtherVec.XY()*anoffset));
 625       gp_Vec2d acurvec(theLastPoint, acurpoint);
 626       if ( bIsOut )
 627         acurvec.Reverse();
 628
 629       Standard_Real aDotX, anAngleX;
 630       //
 631       aDotX = aVec.Dot(acurvec);
 632       anAngleX = aVec.Angle(acurvec);
 633       //
 634       if(aDotX > 0. && fabs(anAngleX) < Precision::PConfusion()) {
 635         if((i % 2) == 0) {
 636           if((acurpoint.Y() >= theVmin) &&
 637              (acurpoint.Y() <= theVmax)) {
 638             theNewPoint = acurpoint;
 639             return Standard_True;
 640           }
 641         }
 642         else {
 643           if((acurpoint.X() >= theUmin) &&
 644              (acurpoint.X() <= theUmax)) {
 645             theNewPoint = acurpoint;
 646             return Standard_True;
 647           }
 648         }
 649       }
 650     }
 651   }
 652   return Standard_False;
 653 }
 654
 655 //=======================================================================
 656 //function : FindPoint
 657 //purpose  : Find point on the boundary of radial tangent zone
 658 //           static subfunction in DecompositionOfWLine
 659 //=======================================================================
 660 static
 661 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
 662                            const gp_Pnt2d&     theLastPoint,
 663                            const Standard_Real theUmin,
 664                            const Standard_Real theUmax,
 665                            const Standard_Real theVmin,
 666                            const Standard_Real theVmax,
 667                            const gp_Pnt2d&     theTanZoneCenter,
 668                            const Standard_Real theZoneRadius,
 669                            Handle(GeomAdaptor_HSurface) theGASurface,
 670                            gp_Pnt2d&           theNewPoint) {
 671   theNewPoint = theLastPoint;
 672
 673   if ( !IsInsideTanZone( theLastPoint, theTanZoneCenter, theZoneRadius, theGASurface) )
 674     return Standard_False;
 675
 676   Standard_Real aUResolution = theGASurface->UResolution( theZoneRadius );
 677   Standard_Real aVResolution = theGASurface->VResolution( theZoneRadius );
 678
 679   Standard_Real aRadius = ( aUResolution < aVResolution ) ? aUResolution : aVResolution;
 680   gp_Ax22d anAxis( theTanZoneCenter, gp_Dir2d(1, 0), gp_Dir2d(0, 1) );
 681   gp_Circ2d aCircle( anAxis, aRadius );
 682
 683   //
 684   gp_Vec2d aDir( theLastPoint.XY() - theFirstPoint.XY() );
 685   Standard_Real aLength = aDir.Magnitude();
 686   if ( aLength <= gp::Resolution() )
 687     return Standard_False;
 688   gp_Lin2d aLine( theFirstPoint, aDir );
 689
 690   //
 691   Handle(Geom2d_Line) aCLine = new Geom2d_Line( aLine );
 692   Handle(Geom2d_TrimmedCurve) aC1 = new Geom2d_TrimmedCurve( aCLine, 0, aLength );
 693   Handle(Geom2d_Circle) aC2 = new Geom2d_Circle( aCircle );
 694
 695   Standard_Real aTol = aRadius * 0.001;
 696   aTol = ( aTol < Precision::PConfusion() ) ? Precision::PConfusion() : aTol;
 697
 698   Geom2dAPI_InterCurveCurve anIntersector;
 699   anIntersector.Init( aC1, aC2, aTol );
 700
 701   if ( anIntersector.NbPoints() == 0 )
 702     return Standard_False;
 703
 704   Standard_Boolean aFound = Standard_False;
 705   Standard_Real aMinDist = aLength * aLength;
 706   Standard_Integer i = 0;
 707   for ( i = 1; i <= anIntersector.NbPoints(); i++ ) {
 708     gp_Pnt2d aPInt = anIntersector.Point( i );
 709     if ( aPInt.SquareDistance( theFirstPoint ) < aMinDist ) {
 710       if ( ( aPInt.X() >= theUmin ) && ( aPInt.X() <= theUmax ) &&
 711            ( aPInt.Y() >= theVmin ) && ( aPInt.Y() <= theVmax ) ) {
 712         theNewPoint = aPInt;
 713         aFound = Standard_True;
 714       }
 715     }
 716   }
 717
 718   return aFound;
 719 }
 720
 721 //=======================================================================
 722 //function : DecompositionOfWLine
 723 //purpose  :
 724 //=======================================================================
 725 Standard_Boolean IntTools_WLineTool::
 726   DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
 727                        const Handle(GeomAdaptor_HSurface)&            theSurface1,
 728                        const Handle(GeomAdaptor_HSurface)&            theSurface2,
 729                        const TopoDS_Face&                             theFace1,
 730                        const TopoDS_Face&                             theFace2,
 731                        const GeomInt_LineConstructor&                 theLConstructor,
 732                        const Standard_Boolean                         theAvoidLConstructor,
 733                        const Standard_Real                            theTol,
 734                        IntPatch_SequenceOfLine&                       theNewLines,
 735                        Standard_Real&                                 theReachedTol3d,
 736                        const Handle(IntTools_Context)& aContext)
 737 {
 738   Standard_Boolean bRet, bAvoidLineConstructor;
 739   Standard_Integer aNbPnts, aNbParts;
 740   //
 741   bRet=Standard_False;
 742   aNbPnts=theWLine->NbPnts();
 743   bAvoidLineConstructor=theAvoidLConstructor;
 744   //
 745   if(!aNbPnts){
 746     return bRet;
 747   }
 748   if (!bAvoidLineConstructor) {
 749     aNbParts=theLConstructor.NbParts();
 750     if (!aNbParts) {
 751       return bRet;
 752     }
 753   }
 754   //
 755   Standard_Boolean bIsPrevPointOnBoundary, bIsPointOnBoundary, bIsCurrentPointOnBoundary;
 756   Standard_Integer nblines, pit, i, j;
 757   Standard_Real aTol;
 758   TColStd_Array1OfListOfInteger anArrayOfLines(1, aNbPnts);
 759   TColStd_Array1OfInteger       anArrayOfLineType(1, aNbPnts);
 760   TColStd_ListOfInteger aListOfPointIndex;
 761
 762   Handle(TColgp_HArray1OfPnt2d) aTanZoneS1;
 763   Handle(TColgp_HArray1OfPnt2d) aTanZoneS2;
 764   Handle(TColStd_HArray1OfReal) aTanZoneRadius;
 765   Standard_Integer aNbZone = ComputeTangentZones( theSurface1, theSurface2, theFace1, theFace2,
 766                                                  aTanZoneS1, aTanZoneS2, aTanZoneRadius, aContext);
 767
 768   //
 769   nblines=0;
 770   aTol=Precision::Confusion();
 771   aTol=0.5*aTol;
 772   bIsPrevPointOnBoundary=Standard_False;
 773   bIsPointOnBoundary=Standard_False;
 774   //
 775   // 1. ...
 776   //
 777   // Points
 778   for(pit = 1; pit <= aNbPnts; ++pit) {
 779     Standard_Boolean bIsOnFirstBoundary, isperiodic;
 780     Standard_Real aResolution, aPeriod, alowerboundary, aupperboundary, U, V;
 781     Standard_Real aParameter, anoffset, anAdjustPar;
 782     Standard_Real umin, umax, vmin, vmax;
 783     //
 784     bIsCurrentPointOnBoundary = Standard_False;
 785     const IntSurf_PntOn2S& aPoint = theWLine->Point(pit);
 786     //
 787     // Surface
 788     for(i = 0; i < 2; ++i) {
 789       Handle(GeomAdaptor_HSurface) aGASurface = (!i) ? theSurface1 : theSurface2;
 790       aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
 791       if(!i) {
 792         aPoint.ParametersOnS1(U, V);
 793       }
 794       else {
 795         aPoint.ParametersOnS2(U, V);
 796       }
 797       // U, V
 798       for(j = 0; j < 2; j++) {
 799         isperiodic = (!j) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
 800         if(!isperiodic){
 801           continue;
 802         }
 803         //
 804         if (!j) {
 805           aResolution=aGASurface->UResolution(aTol);
 806           aPeriod=aGASurface->UPeriod();
 807           alowerboundary=umin;
 808           aupperboundary=umax;
 809           aParameter=U;
 810         }
 811         else {
 812           aResolution=aGASurface->VResolution(aTol);
 813           aPeriod=aGASurface->VPeriod();
 814           alowerboundary=vmin;
 815           aupperboundary=vmax;
 816           aParameter=V;
 817         }
 818
 819         GeomInt::AdjustPeriodic(aParameter,
 820                                        alowerboundary,
 821                                        aupperboundary,
 822                                        aPeriod,
 823                                        anAdjustPar,
 824                                        anoffset);
 825         //
 826         bIsOnFirstBoundary = Standard_True;// ?
 827         bIsPointOnBoundary=
 828           IsPointOnBoundary(anAdjustPar,
 829                             alowerboundary,
 830                             aupperboundary,
 831                             aResolution,
 832                             bIsOnFirstBoundary);
 833         //
 834         if(bIsPointOnBoundary) {
 835           bIsCurrentPointOnBoundary = Standard_True;
 836           break;
 837         }
 838         else {
 839           // check if a point belong to a tangent zone. Begin
 840           Standard_Integer zIt = 0;
 841           for ( zIt = 1; zIt <= aNbZone; zIt++ ) {
 842             gp_Pnt2d aPZone = (i == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt);
 843             Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt);
 844
 845             if ( IsInsideTanZone(gp_Pnt2d( U, V ), aPZone, aZoneRadius, aGASurface ) ) {
 846               // set boundary flag to split the curve by a tangent zone
 847               bIsPointOnBoundary = Standard_True;
 848               bIsCurrentPointOnBoundary = Standard_True;
 849               if ( theReachedTol3d < aZoneRadius ) {
 850                 theReachedTol3d = aZoneRadius;
 851               }
 852               break;
 853             }
 854           }
 855         }
 856       }//for(j = 0; j < 2; j++) {
 857
 858       if(bIsCurrentPointOnBoundary){
 859         break;
 860       }
 861     }//for(i = 0; i < 2; ++i) {
 862     //
 863     if((bIsCurrentPointOnBoundary != bIsPrevPointOnBoundary)) {
 864       if(!aListOfPointIndex.IsEmpty()) {
 865         nblines++;
 866         anArrayOfLines.SetValue(nblines, aListOfPointIndex);
 867         anArrayOfLineType.SetValue(nblines, bIsPrevPointOnBoundary);
 868         aListOfPointIndex.Clear();
 869       }
 870       bIsPrevPointOnBoundary = bIsCurrentPointOnBoundary;
 871     }
 872     aListOfPointIndex.Append(pit);
 873   } //for(pit = 1; pit <= aNbPnts; ++pit) {
 874   //
 875   if(!aListOfPointIndex.IsEmpty()) {
 876     nblines++;
 877     anArrayOfLines.SetValue(nblines, aListOfPointIndex);
 878     anArrayOfLineType.SetValue(nblines, bIsPrevPointOnBoundary);
 879     aListOfPointIndex.Clear();
 880   }
 881   //
 882   if(nblines<=1) {
 883     return bRet; //Standard_False;
 884   }
 885   //
 886   //
 887   // 2. Correct wlines.begin
 888   TColStd_Array1OfListOfInteger anArrayOfLineEnds(1, nblines);
 889   Handle(IntSurf_LineOn2S) aSeqOfPntOn2S = new IntSurf_LineOn2S();
 890   //
 891   for(i = 1; i <= nblines; i++) {
 892     if(anArrayOfLineType.Value(i) != 0) {
 893       continue;
 894     }
 895     const TColStd_ListOfInteger& aListOfIndex = anArrayOfLines.Value(i);
 896     TColStd_ListOfInteger aListOfFLIndex;
 897
 898     for(j = 0; j < 2; j++) {
 899       Standard_Integer aneighbourindex = (j == 0) ? (i - 1) : (i + 1);
 900
 901       if((aneighbourindex < 1) || (aneighbourindex > nblines))
 902         continue;
 903
 904       if(anArrayOfLineType.Value(aneighbourindex) == 0)
 905         continue;
 906       const TColStd_ListOfInteger& aNeighbour = anArrayOfLines.Value(aneighbourindex);
 907       Standard_Integer anIndex = (j == 0) ? aNeighbour.Last() : aNeighbour.First();
 908       const IntSurf_PntOn2S& aPoint = theWLine->Point(anIndex);
 909
 910       IntSurf_PntOn2S aNewP = aPoint;
 911       if(aListOfIndex.Extent() < 2) {
 912         aSeqOfPntOn2S->Add(aNewP);
 913         aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
 914         continue;
 915       }
 916       //
 917       Standard_Integer iFirst = aListOfIndex.First();
 918       Standard_Integer iLast  = aListOfIndex.Last();
 919       //
 920       for(Standard_Integer surfit = 0; surfit < 2; surfit++) {
 921
 922         Handle(GeomAdaptor_HSurface) aGASurface = (surfit == 0) ? theSurface1 : theSurface2;
 923         Standard_Real umin=0., umax=0., vmin=0., vmax=0.;
 924         aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
 925         Standard_Real U=0., V=0.;
 926
 927         if(surfit == 0)
 928           aNewP.ParametersOnS1(U, V);
 929         else
 930           aNewP.ParametersOnS2(U, V);
 931         Standard_Integer nbboundaries = 0;
 932
 933         Standard_Boolean bIsNearBoundary = Standard_False;
 934         Standard_Integer aZoneIndex = 0;
 935         Standard_Integer bIsUBoundary = Standard_False; // use if nbboundaries == 1
 936         Standard_Integer bIsFirstBoundary = Standard_False; // use if nbboundaries == 1
 937
 938
 939         for(Standard_Integer parit = 0; parit < 2; parit++) {
 940           Standard_Boolean isperiodic = (parit == 0) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
 941
 942           Standard_Real aResolution = (parit == 0) ? aGASurface->UResolution(aTol) : aGASurface->VResolution(aTol);
 943           Standard_Real alowerboundary = (parit == 0) ? umin : vmin;
 944           Standard_Real aupperboundary = (parit == 0) ? umax : vmax;
 945
 946           Standard_Real aParameter = (parit == 0) ? U : V;
 947           Standard_Boolean bIsOnFirstBoundary = Standard_True;
 948
 949           if(!isperiodic) {
 950             bIsPointOnBoundary=
 951               IsPointOnBoundary(aParameter, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
 952             if(bIsPointOnBoundary) {
 953               bIsUBoundary = (parit == 0);
 954               bIsFirstBoundary = bIsOnFirstBoundary;
 955               nbboundaries++;
 956             }
 957           }
 958           else {
 959             Standard_Real aPeriod     = (parit == 0) ? aGASurface->UPeriod() : aGASurface->VPeriod();
 960             Standard_Real anoffset, anAdjustPar;
 961             GeomInt::AdjustPeriodic(aParameter, alowerboundary, aupperboundary,
 962                                            aPeriod, anAdjustPar, anoffset);
 963
 964             bIsPointOnBoundary=
 965               IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
 966             if(bIsPointOnBoundary) {
 967               bIsUBoundary = (parit == 0);
 968               bIsFirstBoundary = bIsOnFirstBoundary;
 969               nbboundaries++;
 970             }
 971             else {
 972             //check neighbourhood of boundary
 973             Standard_Real anEpsilon = aResolution * 100.;
 974             Standard_Real aPart = ( aupperboundary - alowerboundary ) * 0.1;
 975             anEpsilon = ( anEpsilon > aPart ) ? aPart : anEpsilon;
 976
 977               bIsNearBoundary = IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary,
 978                                                   anEpsilon, bIsOnFirstBoundary);
 979
 980             }
 981           }
 982         }
 983
 984         // check if a point belong to a tangent zone. Begin
 985         for ( Standard_Integer zIt = 1; zIt <= aNbZone; zIt++ ) {
 986           gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt);
 987           Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt);
 988
 989           Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
 990           const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
 991           Standard_Real nU1, nV1;
 992
 993           if(surfit == 0)
 994             aNeighbourPoint.ParametersOnS1(nU1, nV1);
 995           else
 996             aNeighbourPoint.ParametersOnS2(nU1, nV1);
 997           gp_Pnt2d ap1(nU1, nV1);
 998           gp_Pnt2d ap2 = AdjustByNeighbour( ap1, gp_Pnt2d( U, V ), aGASurface );
 999
1000
1001           if ( IsInsideTanZone( ap2, aPZone, aZoneRadius, aGASurface ) ) {
1002             aZoneIndex = zIt;
1003             bIsNearBoundary = Standard_True;
1004             if ( theReachedTol3d < aZoneRadius ) {
1005               theReachedTol3d = aZoneRadius;
1006             }
1007           }
1008         }
1009         // check if a point belong to a tangent zone. End
1010         Standard_Boolean bComputeLineEnd = Standard_False;
1011
1012         if(nbboundaries == 2) {
1013           //xf
1014           bComputeLineEnd = Standard_True;
1015           //xt
1016         }
1017         else if(nbboundaries == 1) {
1018           Standard_Boolean isperiodic = (bIsUBoundary) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
1019
1020           if(isperiodic) {
1021             Standard_Real alowerboundary = (bIsUBoundary) ? umin : vmin;
1022             Standard_Real aupperboundary = (bIsUBoundary) ? umax : vmax;
1023             Standard_Real aPeriod     = (bIsUBoundary) ? aGASurface->UPeriod() : aGASurface->VPeriod();
1024             Standard_Real aParameter = (bIsUBoundary) ? U : V;
1025             Standard_Real anoffset, anAdjustPar;
1026             GeomInt::AdjustPeriodic(aParameter, alowerboundary, aupperboundary,
1027                                            aPeriod, anAdjustPar, anoffset);
1028
1029             Standard_Real adist = (bIsFirstBoundary) ? fabs(anAdjustPar - alowerboundary) : fabs(anAdjustPar - aupperboundary);
1030             Standard_Real anotherPar = (bIsFirstBoundary) ? (aupperboundary - adist) : (alowerboundary + adist);
1031             anotherPar += anoffset;
1032             Standard_Integer aneighbourpointindex = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
1033             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex);
1034             Standard_Real nU1, nV1;
1035
1036             if(surfit == 0)
1037               aNeighbourPoint.ParametersOnS1(nU1, nV1);
1038             else
1039               aNeighbourPoint.ParametersOnS2(nU1, nV1);
1040
1041             Standard_Real adist1 = (bIsUBoundary) ? fabs(nU1 - U) : fabs(nV1 - V);
1042             Standard_Real adist2 = (bIsUBoundary) ? fabs(nU1 - anotherPar) : fabs(nV1 - anotherPar);
1043             bComputeLineEnd = Standard_True;
1044             Standard_Boolean bCheckAngle1 = Standard_False;
1045             Standard_Boolean bCheckAngle2 = Standard_False;
1046             gp_Vec2d aNewVec;
1047             Standard_Real anewU = (bIsUBoundary) ? anotherPar : U;
1048             Standard_Real anewV = (bIsUBoundary) ? V : anotherPar;
1049
1050             if(((adist1 - adist2) > Precision::PConfusion()) &&
1051                (adist2 < (aPeriod / 4.))) {
1052               bCheckAngle1 = Standard_True;
1053               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(anewU, anewV));
1054
1055               if(aNewVec.SquareMagnitude() < gp::Resolution()) {
1056                 aNewP.SetValue((surfit == 0), anewU, anewV);
1057                 bCheckAngle1 = Standard_False;
1058               }
1059             }
1060             else if(adist1 < (aPeriod / 4.)) {
1061               bCheckAngle2 = Standard_True;
1062               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(U, V));
1063
1064               if(aNewVec.SquareMagnitude() < gp::Resolution()) {
1065                 bCheckAngle2 = Standard_False;
1066               }
1067             }
1068
1069             if(bCheckAngle1 || bCheckAngle2) {
1070               // assume there are at least two points in line (see "if" above)
1071               Standard_Integer anindexother = aneighbourpointindex;
1072
1073               while((anindexother <= iLast) && (anindexother >= iFirst)) {
1074                 anindexother = (j == 0) ? (anindexother + 1) : (anindexother - 1);
1075                 const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(anindexother);
1076                 Standard_Real nU2, nV2;
1077
1078                 if(surfit == 0)
1079                   aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
1080                 else
1081                   aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
1082                 gp_Vec2d aVecOld(gp_Pnt2d(nU2, nV2), gp_Pnt2d(nU1, nV1));
1083
1084                 if(aVecOld.SquareMagnitude() <= gp::Resolution()) {
1085                   continue;
1086                 }
1087                 else {
1088                   Standard_Real anAngle = aNewVec.Angle(aVecOld);
1089
1090                   if((fabs(anAngle) < (M_PI * 0.25)) && (aNewVec.Dot(aVecOld) > 0.)) {
1091
1092                     if(bCheckAngle1) {
1093                       Standard_Real U1, U2, V1, V2;
1094                       IntSurf_PntOn2S atmppoint = aNewP;
1095                       atmppoint.SetValue((surfit == 0), anewU, anewV);
1096                       atmppoint.Parameters(U1, V1, U2, V2);
1097                       gp_Pnt P1 = theSurface1->Value(U1, V1);
1098                       gp_Pnt P2 = theSurface2->Value(U2, V2);
1099                       gp_Pnt P0 = aPoint.Value();
1100
1101                       if(P0.IsEqual(P1, aTol) &&
1102                          P0.IsEqual(P2, aTol) &&
1103                          P1.IsEqual(P2, aTol)) {
1104                         bComputeLineEnd = Standard_False;
1105                         aNewP.SetValue((surfit == 0), anewU, anewV);
1106                       }
1107                     }
1108
1109                     if(bCheckAngle2) {
1110                       bComputeLineEnd = Standard_False;
1111                     }
1112                   }
1113                   break;
1114                 }
1115               } // end while(anindexother...)
1116             }
1117           }
1118         }
1119         else if ( bIsNearBoundary ) {
1120           bComputeLineEnd = Standard_True;
1121         }
1122
1123         if(bComputeLineEnd) {
1124
1125           gp_Pnt2d anewpoint;
1126           Standard_Boolean found = Standard_False;
1127
1128           if ( bIsNearBoundary ) {
1129             // re-compute point near natural boundary or near tangent zone
1130             Standard_Real u1, v1, u2, v2;
1131             aNewP.Parameters( u1, v1, u2, v2 );
1132             if(surfit == 0)
1133               anewpoint = gp_Pnt2d( u1, v1 );
1134             else
1135               anewpoint = gp_Pnt2d( u2, v2 );
1136
1137             Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
1138             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
1139             Standard_Real nU1, nV1;
1140
1141             if(surfit == 0)
1142               aNeighbourPoint.ParametersOnS1(nU1, nV1);
1143             else
1144               aNeighbourPoint.ParametersOnS2(nU1, nV1);
1145             gp_Pnt2d ap1(nU1, nV1);
1146             gp_Pnt2d ap2;
1147
1148
1149             if ( aZoneIndex ) {
1150               // exclude point from a tangent zone
1151               anewpoint = AdjustByNeighbour( ap1, anewpoint, aGASurface );
1152               gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(aZoneIndex) : aTanZoneS2->Value(aZoneIndex);
1153               Standard_Real aZoneRadius = aTanZoneRadius->Value(aZoneIndex);
1154
1155               if ( FindPoint(ap1, anewpoint, umin, umax, vmin, vmax,
1156                              aPZone, aZoneRadius, aGASurface, ap2) ) {
1157                 anewpoint = ap2;
1158                 found = Standard_True;
1159               }
1160             }
1161             else if ( aGASurface->IsUPeriodic() || aGASurface->IsVPeriodic() ) {
1162               // re-compute point near boundary if shifted on a period
1163               ap2 = AdjustByNeighbour( ap1, anewpoint, aGASurface );
1164
1165               if ( ( ap2.X() < umin ) || ( ap2.X() > umax ) ||
1166                   ( ap2.Y() < vmin ) || ( ap2.Y() > vmax ) ) {
1167                 found = FindPoint(ap1, ap2, umin, umax, vmin, vmax, anewpoint);
1168               }
1169               else {
1170                 anewpoint = ap2;
1171                 aNewP.SetValue( (surfit == 0), anewpoint.X(), anewpoint.Y() );
1172               }
1173             }
1174           }
1175           else {
1176
1177             Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
1178             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
1179             Standard_Real nU1, nV1;
1180
1181             if(surfit == 0)
1182               aNeighbourPoint.ParametersOnS1(nU1, nV1);
1183             else
1184               aNeighbourPoint.ParametersOnS2(nU1, nV1);
1185             gp_Pnt2d ap1(nU1, nV1);
1186             gp_Pnt2d ap2(nU1, nV1);
1187             Standard_Integer aneighbourpointindex2 = aneighbourpointindex1;
1188
1189             while((aneighbourpointindex2 <= iLast) && (aneighbourpointindex2 >= iFirst)) {
1190               aneighbourpointindex2 = (j == 0) ? (aneighbourpointindex2 + 1) : (aneighbourpointindex2 - 1);
1191               const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(aneighbourpointindex2);
1192               Standard_Real nU2, nV2;
1193
1194               if(surfit == 0)
1195                 aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
1196               else
1197                 aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
1198               ap2.SetX(nU2);
1199               ap2.SetY(nV2);
1200
1201               if(ap1.SquareDistance(ap2) > gp::Resolution()) {
1202                 break;
1203               }
1204             }
1205             found = FindPoint(ap2, ap1, umin, umax, vmin, vmax, anewpoint);
1206           }
1207
1208           if(found) {
1209             // check point
1210             Standard_Real aCriteria = theTol;
1211             GeomAPI_ProjectPointOnSurf& aProjector =
1212               (surfit == 0) ? aContext->ProjPS(theFace2) : aContext->ProjPS(theFace1);
1213             Handle(GeomAdaptor_HSurface) aSurface = (surfit == 0) ? theSurface1 : theSurface2;
1214
1215             Handle(GeomAdaptor_HSurface) aSurfaceOther = (surfit == 0) ? theSurface2 : theSurface1;
1216
1217             gp_Pnt aP3d = aSurface->Value(anewpoint.X(), anewpoint.Y());
1218             aProjector.Perform(aP3d);
1219
1220             if(aProjector.IsDone()) {
1221               if(aProjector.LowerDistance() < aCriteria) {
1222                 Standard_Real foundU = U, foundV = V;
1223                 aProjector.LowerDistanceParameters(foundU, foundV);
1224
1225                 //Correction of projected coordinates. Begin
1226                 //Note, it may be shifted on a period
1227                 Standard_Integer aneindex1 = (j == 0) ? iFirst : iLast;
1228                 const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneindex1);
1229                 Standard_Real nUn, nVn;
1230
1231                 if(surfit == 0)
1232                   aNeighbourPoint.ParametersOnS2(nUn, nVn);
1233                 else
1234                   aNeighbourPoint.ParametersOnS1(nUn, nVn);
1235                 gp_Pnt2d aNeighbour2d(nUn, nVn);
1236                 gp_Pnt2d anAdjustedPoint = AdjustByNeighbour( aNeighbour2d, gp_Pnt2d(foundU, foundV), aSurfaceOther );
1237                 foundU = anAdjustedPoint.X();
1238                 foundV = anAdjustedPoint.Y();
1239
1240                 if ( ( anAdjustedPoint.X() < umin ) && ( anAdjustedPoint.X() > umax ) &&
1241                     ( anAdjustedPoint.Y() < vmin ) && ( anAdjustedPoint.Y() > vmax ) ) {
1242                   // attempt to roughly re-compute point
1243                   foundU = ( foundU < umin ) ? umin : foundU;
1244                   foundU = ( foundU > umax ) ? umax : foundU;
1245                   foundV = ( foundV < vmin ) ? vmin : foundV;
1246                   foundV = ( foundV > vmax ) ? vmax : foundV;
1247
1248                   GeomAPI_ProjectPointOnSurf& aProjector2 =
1249                     (surfit == 0) ? aContext->ProjPS(theFace1) : aContext->ProjPS(theFace2);
1250
1251                   aP3d = aSurfaceOther->Value(foundU, foundV);
1252                   aProjector2.Perform(aP3d);
1253
1254                   if(aProjector2.IsDone()) {
1255                     if(aProjector2.LowerDistance() < aCriteria) {
1256                       Standard_Real foundU2 = anewpoint.X(), foundV2 = anewpoint.Y();
1257                       aProjector2.LowerDistanceParameters(foundU2, foundV2);
1258                       anewpoint.SetX(foundU2);
1259                       anewpoint.SetY(foundV2);
1260                     }
1261                   }
1262                 }
1263                 //Correction of projected coordinates. End
1264
1265                 if(surfit == 0)
1266                   aNewP.SetValue(aP3d, anewpoint.X(), anewpoint.Y(), foundU, foundV);
1267                 else
1268                   aNewP.SetValue(aP3d, foundU, foundV, anewpoint.X(), anewpoint.Y());
1269               }
1270             }
1271           }
1272         }
1273       }
1274       aSeqOfPntOn2S->Add(aNewP);
1275       aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
1276     }
1277     anArrayOfLineEnds.SetValue(i, aListOfFLIndex);
1278   }
1279   // Correct wlines.end
1280
1281   // Split wlines.begin
1282   Standard_Integer nbiter;
1283   //
1284   nbiter=1;
1285   if (!bAvoidLineConstructor) {
1286     nbiter=theLConstructor.NbParts();
1287   }
1288   //
1289   for(j = 1; j <= nbiter; ++j) {
1290     Standard_Real fprm, lprm;
1291     Standard_Integer ifprm, ilprm;
1292     //
1293     if(bAvoidLineConstructor) {
1294       ifprm = 1;
1295       ilprm = theWLine->NbPnts();
1296     }
1297     else {
1298       theLConstructor.Part(j, fprm, lprm);
1299       ifprm = (Standard_Integer)fprm;
1300       ilprm = (Standard_Integer)lprm;
1301     }
1302
1303     Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
1304     //
1305     for(i = 1; i <= nblines; i++) {
1306       if(anArrayOfLineType.Value(i) != 0) {
1307         continue;
1308       }
1309       const TColStd_ListOfInteger& aListOfIndex = anArrayOfLines.Value(i);
1310       const TColStd_ListOfInteger& aListOfFLIndex = anArrayOfLineEnds.Value(i);
1311       Standard_Boolean bhasfirstpoint = (aListOfFLIndex.Extent() == 2);
1312       Standard_Boolean bhaslastpoint = (aListOfFLIndex.Extent() == 2);
1313
1314       if(!bhasfirstpoint && !aListOfFLIndex.IsEmpty()) {
1315         bhasfirstpoint = (i != 1);
1316       }
1317
1318       if(!bhaslastpoint && !aListOfFLIndex.IsEmpty()) {
1319         bhaslastpoint = (i != nblines);
1320       }
1321
1322       Standard_Integer iFirst = aListOfIndex.First();
1323       Standard_Integer iLast  = aListOfIndex.Last();
1324       Standard_Boolean bIsFirstInside = ((ifprm >= iFirst) && (ifprm <= iLast));
1325       Standard_Boolean bIsLastInside =  ((ilprm >= iFirst) && (ilprm <= iLast));
1326
1327       if(!bIsFirstInside && !bIsLastInside) {
1328         if((ifprm < iFirst) && (ilprm > iLast)) {
1329           // append whole line, and boundaries if neccesary
1330           if(bhasfirstpoint) {
1331             pit = aListOfFLIndex.First();
1332             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
1333             aLineOn2S->Add(aP);
1334           }
1335           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
1336
1337           for(; anIt.More(); anIt.Next()) {
1338             pit = anIt.Value();
1339             const IntSurf_PntOn2S& aP = theWLine->Point(pit);
1340             aLineOn2S->Add(aP);
1341           }
1342
1343           if(bhaslastpoint) {
1344             pit = aListOfFLIndex.Last();
1345             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
1346             aLineOn2S->Add(aP);
1347           }
1348
1349           // check end of split line (end is almost always)
1350           Standard_Integer aneighbour = i + 1;
1351           Standard_Boolean bIsEndOfLine = Standard_True;
1352
1353           if(aneighbour <= nblines) {
1354             const TColStd_ListOfInteger& aListOfNeighbourIndex = anArrayOfLines.Value(aneighbour);
1355
1356             if((anArrayOfLineType.Value(aneighbour) != 0) &&
1357                (aListOfNeighbourIndex.IsEmpty())) {
1358               bIsEndOfLine = Standard_False;
1359             }
1360           }
1361
1362           if(bIsEndOfLine) {
1363             if(aLineOn2S->NbPoints() > 1) {
1364               Handle(IntPatch_WLine) aNewWLine =
1365                 new IntPatch_WLine(aLineOn2S, Standard_False);
1366               aNewWLine->SetCreatingWayInfo(theWLine->GetCreatingWay());
1367               theNewLines.Append(aNewWLine);
1368             }
1369             aLineOn2S = new IntSurf_LineOn2S();
1370           }
1371         }
1372         continue;
1373       }
1374       // end if(!bIsFirstInside && !bIsLastInside)
1375
1376       if(bIsFirstInside && bIsLastInside) {
1377         // append inside points between ifprm and ilprm
1378         TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
1379
1380         for(; anIt.More(); anIt.Next()) {
1381           pit = anIt.Value();
1382           if((pit < ifprm) || (pit > ilprm))
1383             continue;
1384           const IntSurf_PntOn2S& aP = theWLine->Point(pit);
1385           aLineOn2S->Add(aP);
1386         }
1387       }
1388       else {
1389
1390         if(bIsFirstInside) {
1391           // append points from ifprm to last point + boundary point
1392           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
1393
1394           for(; anIt.More(); anIt.Next()) {
1395             pit = anIt.Value();
1396             if(pit < ifprm)
1397               continue;
1398             const IntSurf_PntOn2S& aP = theWLine->Point(pit);
1399             aLineOn2S->Add(aP);
1400           }
1401
1402           if(bhaslastpoint) {
1403             pit = aListOfFLIndex.Last();
1404             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
1405             aLineOn2S->Add(aP);
1406           }
1407           // check end of split line (end is almost always)
1408           Standard_Integer aneighbour = i + 1;
1409           Standard_Boolean bIsEndOfLine = Standard_True;
1410
1411           if(aneighbour <= nblines) {
1412             const TColStd_ListOfInteger& aListOfNeighbourIndex = anArrayOfLines.Value(aneighbour);
1413
1414             if((anArrayOfLineType.Value(aneighbour) != 0) &&
1415                (aListOfNeighbourIndex.IsEmpty())) {
1416               bIsEndOfLine = Standard_False;
1417             }
1418           }
1419
1420           if(bIsEndOfLine) {
1421             if(aLineOn2S->NbPoints() > 1) {
1422               Handle(IntPatch_WLine) aNewWLine =
1423                 new IntPatch_WLine(aLineOn2S, Standard_False);
1424               aNewWLine->SetCreatingWayInfo(theWLine->GetCreatingWay());
1425               theNewLines.Append(aNewWLine);
1426             }
1427             aLineOn2S = new IntSurf_LineOn2S();
1428           }
1429         }
1430         // end if(bIsFirstInside)
1431
1432         if(bIsLastInside) {
1433           // append points from first boundary point to ilprm
1434           if(bhasfirstpoint) {
1435             pit = aListOfFLIndex.First();
1436             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
1437             aLineOn2S->Add(aP);
1438           }
1439           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
1440
1441           for(; anIt.More(); anIt.Next()) {
1442             pit = anIt.Value();
1443             if(pit > ilprm)
1444               continue;
1445             const IntSurf_PntOn2S& aP = theWLine->Point(pit);
1446             aLineOn2S->Add(aP);
1447           }
1448         }
1449         //end if(bIsLastInside)
1450       }
1451     }
1452
1453     if(aLineOn2S->NbPoints() > 1) {
1454       Handle(IntPatch_WLine) aNewWLine =
1455         new IntPatch_WLine(aLineOn2S, Standard_False);
1456       aNewWLine->SetCreatingWayInfo(theWLine->GetCreatingWay());
1457       theNewLines.Append(aNewWLine);
1458     }
1459   }
1460   // Split wlines.end
1461
1462   return Standard_True;
1463 }
1464
1465 ///////////////////// end of DecompositionOfWLine ///////////////////////