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                        IntPatch_SequenceOfLine&                       theNewLines,
 734                        Standard_Real&                                 theReachedTol3d,
 735                        const Handle(IntTools_Context)& aContext)
 736 {
 737   Standard_Boolean bRet, bAvoidLineConstructor;
 738   Standard_Integer aNbPnts, aNbParts;
 739   //
 740   bRet=Standard_False;
 741   aNbPnts=theWLine->NbPnts();
 742   bAvoidLineConstructor=theAvoidLConstructor;
 743   //
 744   if(!aNbPnts){
 745     return bRet;
 746   }
 747   if (!bAvoidLineConstructor) {
 748     aNbParts=theLConstructor.NbParts();
 749     if (!aNbParts) {
 750       return bRet;
 751     }
 752   }
 753   //
 754   Standard_Boolean bIsPrevPointOnBoundary, bIsPointOnBoundary, bIsCurrentPointOnBoundary;
 755   Standard_Integer nblines, pit, i, j;
 756   Standard_Real aTol;
 757   TColStd_Array1OfListOfInteger anArrayOfLines(1, aNbPnts);
 758   TColStd_Array1OfInteger       anArrayOfLineType(1, aNbPnts);
 759   TColStd_ListOfInteger aListOfPointIndex;
 760
 761   Handle(TColgp_HArray1OfPnt2d) aTanZoneS1;
 762   Handle(TColgp_HArray1OfPnt2d) aTanZoneS2;
 763   Handle(TColStd_HArray1OfReal) aTanZoneRadius;
 764   Standard_Integer aNbZone = ComputeTangentZones( theSurface1, theSurface2, theFace1, theFace2,
 765                                                  aTanZoneS1, aTanZoneS2, aTanZoneRadius, aContext);
 766
 767   //
 768   nblines=0;
 769   aTol=Precision::Confusion();
 770   aTol=0.5*aTol;
 771   bIsPrevPointOnBoundary=Standard_False;
 772   bIsPointOnBoundary=Standard_False;
 773   //
 774   // 1. ...
 775   //
 776   // Points
 777   for(pit = 1; pit <= aNbPnts; ++pit) {
 778     Standard_Boolean bIsOnFirstBoundary, isperiodic;
 779     Standard_Real aResolution, aPeriod, alowerboundary, aupperboundary, U, V;
 780     Standard_Real aParameter, anoffset, anAdjustPar;
 781     Standard_Real umin, umax, vmin, vmax;
 782     //
 783     bIsCurrentPointOnBoundary = Standard_False;
 784     const IntSurf_PntOn2S& aPoint = theWLine->Point(pit);
 785     //
 786     // Surface
 787     for(i = 0; i < 2; ++i) {
 788       Handle(GeomAdaptor_HSurface) aGASurface = (!i) ? theSurface1 : theSurface2;
 789       aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
 790       if(!i) {
 791         aPoint.ParametersOnS1(U, V);
 792       }
 793       else {
 794         aPoint.ParametersOnS2(U, V);
 795       }
 796       // U, V
 797       for(j = 0; j < 2; j++) {
 798         isperiodic = (!j) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
 799         if(!isperiodic){
 800           continue;
 801         }
 802         //
 803         if (!j) {
 804           aResolution=aGASurface->UResolution(aTol);
 805           aPeriod=aGASurface->UPeriod();
 806           alowerboundary=umin;
 807           aupperboundary=umax;
 808           aParameter=U;
 809         }
 810         else {
 811           aResolution=aGASurface->VResolution(aTol);
 812           aPeriod=aGASurface->VPeriod();
 813           alowerboundary=vmin;
 814           aupperboundary=vmax;
 815           aParameter=V;
 816         }
 817
 818         GeomInt::AdjustPeriodic(aParameter,
 819                                        alowerboundary,
 820                                        aupperboundary,
 821                                        aPeriod,
 822                                        anAdjustPar,
 823                                        anoffset);
 824         //
 825         bIsOnFirstBoundary = Standard_True;// ?
 826         bIsPointOnBoundary=
 827           IsPointOnBoundary(anAdjustPar,
 828                             alowerboundary,
 829                             aupperboundary,
 830                             aResolution,
 831                             bIsOnFirstBoundary);
 832         //
 833         if(bIsPointOnBoundary) {
 834           bIsCurrentPointOnBoundary = Standard_True;
 835           break;
 836         }
 837         else {
 838           // check if a point belong to a tangent zone. Begin
 839           Standard_Integer zIt = 0;
 840           for ( zIt = 1; zIt <= aNbZone; zIt++ ) {
 841             gp_Pnt2d aPZone = (i == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt);
 842             Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt);
 843
 844             if ( IsInsideTanZone(gp_Pnt2d( U, V ), aPZone, aZoneRadius, aGASurface ) ) {
 845               // set boundary flag to split the curve by a tangent zone
 846               bIsPointOnBoundary = Standard_True;
 847               bIsCurrentPointOnBoundary = Standard_True;
 848               if ( theReachedTol3d < aZoneRadius ) {
 849                 theReachedTol3d = aZoneRadius;
 850               }
 851               break;
 852             }
 853           }
 854         }
 855       }//for(j = 0; j < 2; j++) {
 856
 857       if(bIsCurrentPointOnBoundary){
 858         break;
 859       }
 860     }//for(i = 0; i < 2; ++i) {
 861     //
 862     if((bIsCurrentPointOnBoundary != bIsPrevPointOnBoundary)) {
 863       if(!aListOfPointIndex.IsEmpty()) {
 864         nblines++;
 865         anArrayOfLines.SetValue(nblines, aListOfPointIndex);
 866         anArrayOfLineType.SetValue(nblines, bIsPrevPointOnBoundary);
 867         aListOfPointIndex.Clear();
 868       }
 869       bIsPrevPointOnBoundary = bIsCurrentPointOnBoundary;
 870     }
 871     aListOfPointIndex.Append(pit);
 872   } //for(pit = 1; pit <= aNbPnts; ++pit) {
 873   //
 874   if(!aListOfPointIndex.IsEmpty()) {
 875     nblines++;
 876     anArrayOfLines.SetValue(nblines, aListOfPointIndex);
 877     anArrayOfLineType.SetValue(nblines, bIsPrevPointOnBoundary);
 878     aListOfPointIndex.Clear();
 879   }
 880   //
 881   if(nblines<=1) {
 882     return bRet; //Standard_False;
 883   }
 884   //
 885   //
 886   // 2. Correct wlines.begin
 887   TColStd_Array1OfListOfInteger anArrayOfLineEnds(1, nblines);
 888   Handle(IntSurf_LineOn2S) aSeqOfPntOn2S = new IntSurf_LineOn2S();
 889   //
 890   for(i = 1; i <= nblines; i++) {
 891     if(anArrayOfLineType.Value(i) != 0) {
 892       continue;
 893     }
 894     const TColStd_ListOfInteger& aListOfIndex = anArrayOfLines.Value(i);
 895     TColStd_ListOfInteger aListOfFLIndex;
 896
 897     for(j = 0; j < 2; j++) {
 898       Standard_Integer aneighbourindex = (j == 0) ? (i - 1) : (i + 1);
 899
 900       if((aneighbourindex < 1) || (aneighbourindex > nblines))
 901         continue;
 902
 903       if(anArrayOfLineType.Value(aneighbourindex) == 0)
 904         continue;
 905       const TColStd_ListOfInteger& aNeighbour = anArrayOfLines.Value(aneighbourindex);
 906       Standard_Integer anIndex = (j == 0) ? aNeighbour.Last() : aNeighbour.First();
 907       const IntSurf_PntOn2S& aPoint = theWLine->Point(anIndex);
 908
 909       IntSurf_PntOn2S aNewP = aPoint;
 910       if(aListOfIndex.Extent() < 2) {
 911         aSeqOfPntOn2S->Add(aNewP);
 912         aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
 913         continue;
 914       }
 915       //
 916       Standard_Integer iFirst = aListOfIndex.First();
 917       Standard_Integer iLast  = aListOfIndex.Last();
 918       //
 919       for(Standard_Integer surfit = 0; surfit < 2; surfit++) {
 920
 921         Handle(GeomAdaptor_HSurface) aGASurface = (surfit == 0) ? theSurface1 : theSurface2;
 922         Standard_Real umin=0., umax=0., vmin=0., vmax=0.;
 923         aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
 924         Standard_Real U=0., V=0.;
 925
 926         if(surfit == 0)
 927           aNewP.ParametersOnS1(U, V);
 928         else
 929           aNewP.ParametersOnS2(U, V);
 930         Standard_Integer nbboundaries = 0;
 931
 932         Standard_Boolean bIsNearBoundary = Standard_False;
 933         Standard_Integer aZoneIndex = 0;
 934         Standard_Integer bIsUBoundary = Standard_False; // use if nbboundaries == 1
 935         Standard_Integer bIsFirstBoundary = Standard_False; // use if nbboundaries == 1
 936
 937
 938         for(Standard_Integer parit = 0; parit < 2; parit++) {
 939           Standard_Boolean isperiodic = (parit == 0) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
 940
 941           Standard_Real aResolution = (parit == 0) ? aGASurface->UResolution(aTol) : aGASurface->VResolution(aTol);
 942           Standard_Real alowerboundary = (parit == 0) ? umin : vmin;
 943           Standard_Real aupperboundary = (parit == 0) ? umax : vmax;
 944
 945           Standard_Real aParameter = (parit == 0) ? U : V;
 946           Standard_Boolean bIsOnFirstBoundary = Standard_True;
 947
 948           if(!isperiodic) {
 949             bIsPointOnBoundary=
 950               IsPointOnBoundary(aParameter, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
 951             if(bIsPointOnBoundary) {
 952               bIsUBoundary = (parit == 0);
 953               bIsFirstBoundary = bIsOnFirstBoundary;
 954               nbboundaries++;
 955             }
 956           }
 957           else {
 958             Standard_Real aPeriod     = (parit == 0) ? aGASurface->UPeriod() : aGASurface->VPeriod();
 959             Standard_Real anoffset, anAdjustPar;
 960             GeomInt::AdjustPeriodic(aParameter, alowerboundary, aupperboundary,
 961                                            aPeriod, anAdjustPar, anoffset);
 962
 963             bIsPointOnBoundary=
 964               IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
 965             if(bIsPointOnBoundary) {
 966               bIsUBoundary = (parit == 0);
 967               bIsFirstBoundary = bIsOnFirstBoundary;
 968               nbboundaries++;
 969             }
 970             else {
 971             //check neighbourhood of boundary
 972             Standard_Real anEpsilon = aResolution * 100.;
 973             Standard_Real aPart = ( aupperboundary - alowerboundary ) * 0.1;
 974             anEpsilon = ( anEpsilon > aPart ) ? aPart : anEpsilon;
 975
 976               bIsNearBoundary = IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary,
 977                                                   anEpsilon, bIsOnFirstBoundary);
 978
 979             }
 980           }
 981         }
 982
 983         // check if a point belong to a tangent zone. Begin
 984         for ( Standard_Integer zIt = 1; zIt <= aNbZone; zIt++ ) {
 985           gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt);
 986           Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt);
 987
 988           Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
 989           const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
 990           Standard_Real nU1, nV1;
 991
 992           if(surfit == 0)
 993             aNeighbourPoint.ParametersOnS1(nU1, nV1);
 994           else
 995             aNeighbourPoint.ParametersOnS2(nU1, nV1);
 996           gp_Pnt2d ap1(nU1, nV1);
 997           gp_Pnt2d ap2 = AdjustByNeighbour( ap1, gp_Pnt2d( U, V ), aGASurface );
 998
 999
1000           if ( IsInsideTanZone( ap2, aPZone, aZoneRadius, aGASurface ) ) {
1001             aZoneIndex = zIt;
1002             bIsNearBoundary = Standard_True;
1003             if ( theReachedTol3d < aZoneRadius ) {
1004               theReachedTol3d = aZoneRadius;
1005             }
1006           }
1007         }
1008         // check if a point belong to a tangent zone. End
1009         Standard_Boolean bComputeLineEnd = Standard_False;
1010
1011         if(nbboundaries == 2) {
1012           //xf
1013           bComputeLineEnd = Standard_True;
1014           //xt
1015         }
1016         else if(nbboundaries == 1) {
1017           Standard_Boolean isperiodic = (bIsUBoundary) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
1018
1019           if(isperiodic) {
1020             Standard_Real alowerboundary = (bIsUBoundary) ? umin : vmin;
1021             Standard_Real aupperboundary = (bIsUBoundary) ? umax : vmax;
1022             Standard_Real aPeriod     = (bIsUBoundary) ? aGASurface->UPeriod() : aGASurface->VPeriod();
1023             Standard_Real aParameter = (bIsUBoundary) ? U : V;
1024             Standard_Real anoffset, anAdjustPar;
1025             GeomInt::AdjustPeriodic(aParameter, alowerboundary, aupperboundary,
1026                                            aPeriod, anAdjustPar, anoffset);
1027
1028             Standard_Real adist = (bIsFirstBoundary) ? fabs(anAdjustPar - alowerboundary) : fabs(anAdjustPar - aupperboundary);
1029             Standard_Real anotherPar = (bIsFirstBoundary) ? (aupperboundary - adist) : (alowerboundary + adist);
1030             anotherPar += anoffset;
1031             Standard_Integer aneighbourpointindex = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
1032             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex);
1033             Standard_Real nU1, nV1;
1034
1035             if(surfit == 0)
1036               aNeighbourPoint.ParametersOnS1(nU1, nV1);
1037             else
1038               aNeighbourPoint.ParametersOnS2(nU1, nV1);
1039
1040             Standard_Real adist1 = (bIsUBoundary) ? fabs(nU1 - U) : fabs(nV1 - V);
1041             Standard_Real adist2 = (bIsUBoundary) ? fabs(nU1 - anotherPar) : fabs(nV1 - anotherPar);
1042             bComputeLineEnd = Standard_True;
1043             Standard_Boolean bCheckAngle1 = Standard_False;
1044             Standard_Boolean bCheckAngle2 = Standard_False;
1045             gp_Vec2d aNewVec;
1046             Standard_Real anewU = (bIsUBoundary) ? anotherPar : U;
1047             Standard_Real anewV = (bIsUBoundary) ? V : anotherPar;
1048
1049             if(((adist1 - adist2) > Precision::PConfusion()) &&
1050                (adist2 < (aPeriod / 4.))) {
1051               bCheckAngle1 = Standard_True;
1052               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(anewU, anewV));
1053
1054               if(aNewVec.SquareMagnitude() < gp::Resolution()) {
1055                 aNewP.SetValue((surfit == 0), anewU, anewV);
1056                 bCheckAngle1 = Standard_False;
1057               }
1058             }
1059             else if(adist1 < (aPeriod / 4.)) {
1060               bCheckAngle2 = Standard_True;
1061               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(U, V));
1062
1063               if(aNewVec.SquareMagnitude() < gp::Resolution()) {
1064                 bCheckAngle2 = Standard_False;
1065               }
1066             }
1067
1068             if(bCheckAngle1 || bCheckAngle2) {
1069               // assume there are at least two points in line (see "if" above)
1070               Standard_Integer anindexother = aneighbourpointindex;
1071
1072               while((anindexother <= iLast) && (anindexother >= iFirst)) {
1073                 anindexother = (j == 0) ? (anindexother + 1) : (anindexother - 1);
1074                 const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(anindexother);
1075                 Standard_Real nU2, nV2;
1076
1077                 if(surfit == 0)
1078                   aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
1079                 else
1080                   aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
1081                 gp_Vec2d aVecOld(gp_Pnt2d(nU2, nV2), gp_Pnt2d(nU1, nV1));
1082
1083                 if(aVecOld.SquareMagnitude() <= gp::Resolution()) {
1084                   continue;
1085                 }
1086                 else {
1087                   Standard_Real anAngle = aNewVec.Angle(aVecOld);
1088
1089                   if((fabs(anAngle) < (M_PI * 0.25)) && (aNewVec.Dot(aVecOld) > 0.)) {
1090
1091                     if(bCheckAngle1) {
1092                       Standard_Real U1, U2, V1, V2;
1093                       IntSurf_PntOn2S atmppoint = aNewP;
1094                       atmppoint.SetValue((surfit == 0), anewU, anewV);
1095                       atmppoint.Parameters(U1, V1, U2, V2);
1096                       gp_Pnt P1 = theSurface1->Value(U1, V1);
1097                       gp_Pnt P2 = theSurface2->Value(U2, V2);
1098                       gp_Pnt P0 = aPoint.Value();
1099
1100                       if(P0.IsEqual(P1, aTol) &&
1101                          P0.IsEqual(P2, aTol) &&
1102                          P1.IsEqual(P2, aTol)) {
1103                         bComputeLineEnd = Standard_False;
1104                         aNewP.SetValue((surfit == 0), anewU, anewV);
1105                       }
1106                     }
1107
1108                     if(bCheckAngle2) {
1109                       bComputeLineEnd = Standard_False;
1110                     }
1111                   }
1112                   break;
1113                 }
1114               } // end while(anindexother...)
1115             }
1116           }
1117         }
1118         else if ( bIsNearBoundary ) {
1119           bComputeLineEnd = Standard_True;
1120         }
1121
1122         if(bComputeLineEnd) {
1123
1124           gp_Pnt2d anewpoint;
1125           Standard_Boolean found = Standard_False;
1126
1127           if ( bIsNearBoundary ) {
1128             // re-compute point near natural boundary or near tangent zone
1129             Standard_Real u1, v1, u2, v2;
1130             aNewP.Parameters( u1, v1, u2, v2 );
1131             if(surfit == 0)
1132               anewpoint = gp_Pnt2d( u1, v1 );
1133             else
1134               anewpoint = gp_Pnt2d( u2, v2 );
1135
1136             Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
1137             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
1138             Standard_Real nU1, nV1;
1139
1140             if(surfit == 0)
1141               aNeighbourPoint.ParametersOnS1(nU1, nV1);
1142             else
1143               aNeighbourPoint.ParametersOnS2(nU1, nV1);
1144             gp_Pnt2d ap1(nU1, nV1);
1145             gp_Pnt2d ap2;
1146
1147
1148             if ( aZoneIndex ) {
1149               // exclude point from a tangent zone
1150               anewpoint = AdjustByNeighbour( ap1, anewpoint, aGASurface );
1151               gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(aZoneIndex) : aTanZoneS2->Value(aZoneIndex);
1152               Standard_Real aZoneRadius = aTanZoneRadius->Value(aZoneIndex);
1153
1154               if ( FindPoint(ap1, anewpoint, umin, umax, vmin, vmax,
1155                              aPZone, aZoneRadius, aGASurface, ap2) ) {
1156                 anewpoint = ap2;
1157                 found = Standard_True;
1158               }
1159             }
1160             else if ( aGASurface->IsUPeriodic() || aGASurface->IsVPeriodic() ) {
1161               // re-compute point near boundary if shifted on a period
1162               ap2 = AdjustByNeighbour( ap1, anewpoint, aGASurface );
1163
1164               if ( ( ap2.X() < umin ) || ( ap2.X() > umax ) ||
1165                   ( ap2.Y() < vmin ) || ( ap2.Y() > vmax ) ) {
1166                 found = FindPoint(ap1, ap2, umin, umax, vmin, vmax, anewpoint);
1167               }
1168               else {
1169                 anewpoint = ap2;
1170                 aNewP.SetValue( (surfit == 0), anewpoint.X(), anewpoint.Y() );
1171               }
1172             }
1173           }
1174           else {
1175
1176             Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
1177             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
1178             Standard_Real nU1, nV1;
1179
1180             if(surfit == 0)
1181               aNeighbourPoint.ParametersOnS1(nU1, nV1);
1182             else
1183               aNeighbourPoint.ParametersOnS2(nU1, nV1);
1184             gp_Pnt2d ap1(nU1, nV1);
1185             gp_Pnt2d ap2(nU1, nV1);
1186             Standard_Integer aneighbourpointindex2 = aneighbourpointindex1;
1187
1188             while((aneighbourpointindex2 <= iLast) && (aneighbourpointindex2 >= iFirst)) {
1189               aneighbourpointindex2 = (j == 0) ? (aneighbourpointindex2 + 1) : (aneighbourpointindex2 - 1);
1190               const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(aneighbourpointindex2);
1191               Standard_Real nU2, nV2;
1192
1193               if(surfit == 0)
1194                 aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
1195               else
1196                 aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
1197               ap2.SetX(nU2);
1198               ap2.SetY(nV2);
1199
1200               if(ap1.SquareDistance(ap2) > gp::Resolution()) {
1201                 break;
1202               }
1203             }
1204             found = FindPoint(ap2, ap1, umin, umax, vmin, vmax, anewpoint);
1205           }
1206
1207           if(found) {
1208             // check point
1209             Standard_Real aCriteria = BRep_Tool::Tolerance(theFace1) + BRep_Tool::Tolerance(theFace2);
1210             GeomAPI_ProjectPointOnSurf& aProjector =
1211               (surfit == 0) ? aContext->ProjPS(theFace2) : aContext->ProjPS(theFace1);
1212             Handle(GeomAdaptor_HSurface) aSurface = (surfit == 0) ? theSurface1 : theSurface2;
1213
1214             Handle(GeomAdaptor_HSurface) aSurfaceOther = (surfit == 0) ? theSurface2 : theSurface1;
1215
1216             gp_Pnt aP3d = aSurface->Value(anewpoint.X(), anewpoint.Y());
1217             aProjector.Perform(aP3d);
1218
1219             if(aProjector.IsDone()) {
1220               if(aProjector.LowerDistance() < aCriteria) {
1221                 Standard_Real foundU = U, foundV = V;
1222                 aProjector.LowerDistanceParameters(foundU, foundV);
1223
1224                 //Correction of projected coordinates. Begin
1225                 //Note, it may be shifted on a period
1226                 Standard_Integer aneindex1 = (j == 0) ? iFirst : iLast;
1227                 const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneindex1);
1228                 Standard_Real nUn, nVn;
1229
1230                 if(surfit == 0)
1231                   aNeighbourPoint.ParametersOnS2(nUn, nVn);
1232                 else
1233                   aNeighbourPoint.ParametersOnS1(nUn, nVn);
1234                 gp_Pnt2d aNeighbour2d(nUn, nVn);
1235                 gp_Pnt2d anAdjustedPoint = AdjustByNeighbour( aNeighbour2d, gp_Pnt2d(foundU, foundV), aSurfaceOther );
1236                 foundU = anAdjustedPoint.X();
1237                 foundV = anAdjustedPoint.Y();
1238
1239                 if ( ( anAdjustedPoint.X() < umin ) && ( anAdjustedPoint.X() > umax ) &&
1240                     ( anAdjustedPoint.Y() < vmin ) && ( anAdjustedPoint.Y() > vmax ) ) {
1241                   // attempt to roughly re-compute point
1242                   foundU = ( foundU < umin ) ? umin : foundU;
1243                   foundU = ( foundU > umax ) ? umax : foundU;
1244                   foundV = ( foundV < vmin ) ? vmin : foundV;
1245                   foundV = ( foundV > vmax ) ? vmax : foundV;
1246
1247                   GeomAPI_ProjectPointOnSurf& aProjector2 =
1248                     (surfit == 0) ? aContext->ProjPS(theFace1) : aContext->ProjPS(theFace2);
1249
1250                   aP3d = aSurfaceOther->Value(foundU, foundV);
1251                   aProjector2.Perform(aP3d);
1252
1253                   if(aProjector2.IsDone()) {
1254                     if(aProjector2.LowerDistance() < aCriteria) {
1255                       Standard_Real foundU2 = anewpoint.X(), foundV2 = anewpoint.Y();
1256                       aProjector2.LowerDistanceParameters(foundU2, foundV2);
1257                       anewpoint.SetX(foundU2);
1258                       anewpoint.SetY(foundV2);
1259                     }
1260                   }
1261                 }
1262                 //Correction of projected coordinates. End
1263
1264                 if(surfit == 0)
1265                   aNewP.SetValue(aP3d, anewpoint.X(), anewpoint.Y(), foundU, foundV);
1266                 else
1267                   aNewP.SetValue(aP3d, foundU, foundV, anewpoint.X(), anewpoint.Y());
1268               }
1269             }
1270           }
1271         }
1272       }
1273       aSeqOfPntOn2S->Add(aNewP);
1274       aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
1275     }
1276     anArrayOfLineEnds.SetValue(i, aListOfFLIndex);
1277   }
1278   // Correct wlines.end
1279
1280   // Split wlines.begin
1281   Standard_Integer nbiter;
1282   //
1283   nbiter=1;
1284   if (!bAvoidLineConstructor) {
1285     nbiter=theLConstructor.NbParts();
1286   }
1287   //
1288   for(j = 1; j <= nbiter; ++j) {
1289     Standard_Real fprm, lprm;
1290     Standard_Integer ifprm, ilprm;
1291     //
1292     if(bAvoidLineConstructor) {
1293       ifprm = 1;
1294       ilprm = theWLine->NbPnts();
1295     }
1296     else {
1297       theLConstructor.Part(j, fprm, lprm);
1298       ifprm = (Standard_Integer)fprm;
1299       ilprm = (Standard_Integer)lprm;
1300     }
1301
1302     Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
1303     //
1304     for(i = 1; i <= nblines; i++) {
1305       if(anArrayOfLineType.Value(i) != 0) {
1306         continue;
1307       }
1308       const TColStd_ListOfInteger& aListOfIndex = anArrayOfLines.Value(i);
1309       const TColStd_ListOfInteger& aListOfFLIndex = anArrayOfLineEnds.Value(i);
1310       Standard_Boolean bhasfirstpoint = (aListOfFLIndex.Extent() == 2);
1311       Standard_Boolean bhaslastpoint = (aListOfFLIndex.Extent() == 2);
1312
1313       if(!bhasfirstpoint && !aListOfFLIndex.IsEmpty()) {
1314         bhasfirstpoint = (i != 1);
1315       }
1316
1317       if(!bhaslastpoint && !aListOfFLIndex.IsEmpty()) {
1318         bhaslastpoint = (i != nblines);
1319       }
1320
1321       Standard_Integer iFirst = aListOfIndex.First();
1322       Standard_Integer iLast  = aListOfIndex.Last();
1323       Standard_Boolean bIsFirstInside = ((ifprm >= iFirst) && (ifprm <= iLast));
1324       Standard_Boolean bIsLastInside =  ((ilprm >= iFirst) && (ilprm <= iLast));
1325
1326       if(!bIsFirstInside && !bIsLastInside) {
1327         if((ifprm < iFirst) && (ilprm > iLast)) {
1328           // append whole line, and boundaries if neccesary
1329           if(bhasfirstpoint) {
1330             pit = aListOfFLIndex.First();
1331             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
1332             aLineOn2S->Add(aP);
1333           }
1334           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
1335
1336           for(; anIt.More(); anIt.Next()) {
1337             pit = anIt.Value();
1338             const IntSurf_PntOn2S& aP = theWLine->Point(pit);
1339             aLineOn2S->Add(aP);
1340           }
1341
1342           if(bhaslastpoint) {
1343             pit = aListOfFLIndex.Last();
1344             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
1345             aLineOn2S->Add(aP);
1346           }
1347
1348           // check end of split line (end is almost always)
1349           Standard_Integer aneighbour = i + 1;
1350           Standard_Boolean bIsEndOfLine = Standard_True;
1351
1352           if(aneighbour <= nblines) {
1353             const TColStd_ListOfInteger& aListOfNeighbourIndex = anArrayOfLines.Value(aneighbour);
1354
1355             if((anArrayOfLineType.Value(aneighbour) != 0) &&
1356                (aListOfNeighbourIndex.IsEmpty())) {
1357               bIsEndOfLine = Standard_False;
1358             }
1359           }
1360
1361           if(bIsEndOfLine) {
1362             if(aLineOn2S->NbPoints() > 1) {
1363               Handle(IntPatch_WLine) aNewWLine =
1364                 new IntPatch_WLine(aLineOn2S, Standard_False);
1365               theNewLines.Append(aNewWLine);
1366             }
1367             aLineOn2S = new IntSurf_LineOn2S();
1368           }
1369         }
1370         continue;
1371       }
1372       // end if(!bIsFirstInside && !bIsLastInside)
1373
1374       if(bIsFirstInside && bIsLastInside) {
1375         // append inside points between ifprm and ilprm
1376         TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
1377
1378         for(; anIt.More(); anIt.Next()) {
1379           pit = anIt.Value();
1380           if((pit < ifprm) || (pit > ilprm))
1381             continue;
1382           const IntSurf_PntOn2S& aP = theWLine->Point(pit);
1383           aLineOn2S->Add(aP);
1384         }
1385       }
1386       else {
1387
1388         if(bIsFirstInside) {
1389           // append points from ifprm to last point + boundary point
1390           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
1391
1392           for(; anIt.More(); anIt.Next()) {
1393             pit = anIt.Value();
1394             if(pit < ifprm)
1395               continue;
1396             const IntSurf_PntOn2S& aP = theWLine->Point(pit);
1397             aLineOn2S->Add(aP);
1398           }
1399
1400           if(bhaslastpoint) {
1401             pit = aListOfFLIndex.Last();
1402             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
1403             aLineOn2S->Add(aP);
1404           }
1405           // check end of split line (end is almost always)
1406           Standard_Integer aneighbour = i + 1;
1407           Standard_Boolean bIsEndOfLine = Standard_True;
1408
1409           if(aneighbour <= nblines) {
1410             const TColStd_ListOfInteger& aListOfNeighbourIndex = anArrayOfLines.Value(aneighbour);
1411
1412             if((anArrayOfLineType.Value(aneighbour) != 0) &&
1413                (aListOfNeighbourIndex.IsEmpty())) {
1414               bIsEndOfLine = Standard_False;
1415             }
1416           }
1417
1418           if(bIsEndOfLine) {
1419             if(aLineOn2S->NbPoints() > 1) {
1420               Handle(IntPatch_WLine) aNewWLine =
1421                 new IntPatch_WLine(aLineOn2S, Standard_False);
1422               theNewLines.Append(aNewWLine);
1423             }
1424             aLineOn2S = new IntSurf_LineOn2S();
1425           }
1426         }
1427         // end if(bIsFirstInside)
1428
1429         if(bIsLastInside) {
1430           // append points from first boundary point to ilprm
1431           if(bhasfirstpoint) {
1432             pit = aListOfFLIndex.First();
1433             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
1434             aLineOn2S->Add(aP);
1435           }
1436           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
1437
1438           for(; anIt.More(); anIt.Next()) {
1439             pit = anIt.Value();
1440             if(pit > ilprm)
1441               continue;
1442             const IntSurf_PntOn2S& aP = theWLine->Point(pit);
1443             aLineOn2S->Add(aP);
1444           }
1445         }
1446         //end if(bIsLastInside)
1447       }
1448     }
1449
1450     if(aLineOn2S->NbPoints() > 1) {
1451       Handle(IntPatch_WLine) aNewWLine =
1452         new IntPatch_WLine(aLineOn2S, Standard_False);
1453       theNewLines.Append(aNewWLine);
1454     }
1455   }
1456   // Split wlines.end
1457
1458   return Standard_True;
1459 }
1460
1461 ///////////////////// end of DecompositionOfWLine ///////////////////////