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_Surface.hxx>
  25 #include <GeomInt.hxx>
  26 #include <GeomInt_LineConstructor.hxx>
  27 #include <Geom_Circle.hxx>
  28 #include <Geom_Surface.hxx>
  29 #include <gp_Circ.hxx>
  30 #include <IntTools_Context.hxx>
  31 #include <TColStd_Array1OfListOfInteger.hxx>
  32 #include <TColStd_SequenceOfReal.hxx>
  33 #include <TColgp_SequenceOfPnt2d.hxx>
  34
  35 /////////////////////// NotUseSurfacesForApprox /////////////////////////
  36
  37 // The block is dedicated to determine whether WLine [ifprm, ilprm]
  38 // crosses the degenerated zone on each given surface or not.
  39 // If Yes -> We will not use info about surfaces during approximation
  40 // because inside degenerated zone of the surface the approx. algo.
  41 // uses wrong values of normal, etc., and resulting curve will have
  42 // oscillations that we would not like to have.
  43
  44 //=======================================================================
  45 //function : IsDegeneratedZone
  46 //purpose  : static subfunction in IsDegeneratedZone
  47 //=======================================================================
  48 static
  49 Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
  50                                    const Handle(Geom_Surface)& aS,
  51                                    const Standard_Integer iDir)
  52 {
  53   Standard_Boolean bFlag=Standard_True;
  54   Standard_Real US1, US2, VS1, VS2, dY, dX, d1, d2, dD;
  55   Standard_Real aXm, aYm, aXb, aYb, aXe, aYe;
  56   aS->Bounds(US1, US2, VS1, VS2);
  57
  58   gp_Pnt aPm, aPb, aPe;
  59
  60   aXm=aP2d.X();
  61   aYm=aP2d.Y();
  62
  63   aS->D0(aXm, aYm, aPm);
  64
  65   dX=1.e-5;
  66   dY=1.e-5;
  67   dD=1.e-12;
  68
  69   if (iDir==1) {
  70     aXb=aXm;
  71     aXe=aXm;
  72     aYb=aYm-dY;
  73     if (aYb < VS1) {
  74       aYb=VS1;
  75     }
  76     aYe=aYm+dY;
  77     if (aYe > VS2) {
  78       aYe=VS2;
  79     }
  80     aS->D0(aXb, aYb, aPb);
  81     aS->D0(aXe, aYe, aPe);
  82
  83     d1=aPm.Distance(aPb);
  84     d2=aPm.Distance(aPe);
  85     if (d1 < dD && d2 < dD) {
  86       return bFlag;
  87     }
  88     return !bFlag;
  89   }
  90   //
  91   else if (iDir==2) {
  92     aYb=aYm;
  93     aYe=aYm;
  94     aXb=aXm-dX;
  95     if (aXb < US1) {
  96       aXb=US1;
  97     }
  98     aXe=aXm+dX;
  99     if (aXe > US2) {
 100       aXe=US2;
 101     }
 102     aS->D0(aXb, aYb, aPb);
 103     aS->D0(aXe, aYe, aPe);
 104
 105     d1=aPm.Distance(aPb);
 106     d2=aPm.Distance(aPe);
 107     if (d1 < dD && d2 < dD) {
 108       return bFlag;
 109     }
 110     return !bFlag;
 111   }
 112   return !bFlag;
 113 }
 114
 115 //=======================================================================
 116 //function : IsPointInDegeneratedZone
 117 //purpose  : static subfunction in NotUseSurfacesForApprox
 118 //=======================================================================
 119 static
 120 Standard_Boolean IsPointInDegeneratedZone(const IntSurf_PntOn2S& aP2S,
 121                                           const TopoDS_Face& aF1,
 122                                           const TopoDS_Face& aF2)
 123
 124 {
 125   Standard_Boolean bFlag=Standard_True;
 126   Standard_Real US11, US12, VS11, VS12, US21, US22, VS21, VS22;
 127   Standard_Real U1, V1, U2, V2, aDelta, aD;
 128   gp_Pnt2d aP2d;
 129
 130   Handle(Geom_Surface)aS1 = BRep_Tool::Surface(aF1);
 131   aS1->Bounds(US11, US12, VS11, VS12);
 132   GeomAdaptor_Surface aGAS1(aS1);
 133
 134   Handle(Geom_Surface)aS2 = BRep_Tool::Surface(aF2);
 135   aS1->Bounds(US21, US22, VS21, VS22);
 136   GeomAdaptor_Surface aGAS2(aS2);
 137   //
 138   //const gp_Pnt& aP=aP2S.Value();
 139   aP2S.Parameters(U1, V1, U2, V2);
 140   //
 141   aDelta=1.e-7;
 142   // Check on Surf 1
 143   aD=aGAS1.UResolution(aDelta);
 144   aP2d.SetCoord(U1, V1);
 145   if (fabs(U1-US11) < aD) {
 146     bFlag=IsDegeneratedZone(aP2d, aS1, 1);
 147     if (bFlag) {
 148       return bFlag;
 149     }
 150   }
 151   if (fabs(U1-US12) < aD) {
 152     bFlag=IsDegeneratedZone(aP2d, aS1, 1);
 153     if (bFlag) {
 154       return bFlag;
 155     }
 156   }
 157   aD=aGAS1.VResolution(aDelta);
 158   if (fabs(V1-VS11) < aDelta) {
 159     bFlag=IsDegeneratedZone(aP2d, aS1, 2);
 160     if (bFlag) {
 161       return bFlag;
 162     }
 163   }
 164   if (fabs(V1-VS12) < aDelta) {
 165     bFlag=IsDegeneratedZone(aP2d, aS1, 2);
 166     if (bFlag) {
 167       return bFlag;
 168     }
 169   }
 170   // Check on Surf 2
 171   aD=aGAS2.UResolution(aDelta);
 172   aP2d.SetCoord(U2, V2);
 173   if (fabs(U2-US21) < aDelta) {
 174     bFlag=IsDegeneratedZone(aP2d, aS2, 1);
 175     if (bFlag) {
 176       return bFlag;
 177     }
 178   }
 179   if (fabs(U2-US22) < aDelta) {
 180     bFlag=IsDegeneratedZone(aP2d, aS2, 1);
 181     if (bFlag) {
 182       return bFlag;
 183     }
 184   }
 185   aD=aGAS2.VResolution(aDelta);
 186   if (fabs(V2-VS21) < aDelta) {
 187     bFlag=IsDegeneratedZone(aP2d, aS2, 2);
 188     if (bFlag) {
 189       return bFlag;
 190     }
 191   }
 192   if (fabs(V2-VS22) < aDelta) {
 193     bFlag=IsDegeneratedZone(aP2d, aS2, 2);
 194     if (bFlag) {
 195       return bFlag;
 196     }
 197   }
 198   return !bFlag;
 199 }
 200
 201 //=======================================================================
 202 //function :  NotUseSurfacesForApprox
 203 //purpose  :
 204 //=======================================================================
 205 Standard_Boolean IntTools_WLineTool::NotUseSurfacesForApprox(const TopoDS_Face& aF1,
 206                                                              const TopoDS_Face& aF2,
 207                                                              const Handle(IntPatch_WLine)& WL,
 208                                                              const Standard_Integer ifprm,
 209                                                              const Standard_Integer ilprm)
 210 {
 211   Standard_Boolean bPInDZ;
 212
 213   Handle(IntSurf_LineOn2S) aLineOn2S=WL->Curve();
 214
 215   const IntSurf_PntOn2S& aP2Sfprm=aLineOn2S->Value(ifprm);
 216   bPInDZ=IsPointInDegeneratedZone(aP2Sfprm, aF1, aF2);
 217   if (bPInDZ) {
 218     return bPInDZ;
 219   }
 220
 221   const IntSurf_PntOn2S& aP2Slprm=aLineOn2S->Value(ilprm);
 222   bPInDZ=IsPointInDegeneratedZone(aP2Slprm, aF1, aF2);
 223
 224   return bPInDZ;
 225 }
 226
 227 /////////////////////// end of NotUseSurfacesForApprox //////////////////
 228
 229 /////////////////////// DecompositionOfWLine ////////////////////////////
 230
 231 //=======================================================================
 232 //function : IsPointOnBoundary
 233 //purpose  : static subfunction in DecompositionOfWLine
 234 //=======================================================================
 235 static
 236 Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
 237                                    const Standard_Real theFirstBoundary,
 238                                    const Standard_Real theSecondBoundary,
 239                                    const Standard_Real theResolution,
 240                                    Standard_Boolean&   IsOnFirstBoundary)
 241 {
 242   Standard_Boolean bRet;
 243   Standard_Integer i;
 244   Standard_Real adist;
 245   //
 246   bRet=Standard_False;
 247   for(i = 0; i < 2; ++i) {
 248     IsOnFirstBoundary = (i == 0);
 249     if (IsOnFirstBoundary) {
 250       adist = fabs(theParameter - theFirstBoundary);
 251     }
 252     else {
 253       adist = fabs(theParameter - theSecondBoundary);
 254     }
 255     if(adist < theResolution) {
 256       return !bRet;
 257     }
 258   }
 259   return bRet;
 260 }
 261
 262 //=======================================================================
 263 //function : AdjustByNeighbour
 264 //purpose  : static subfunction in DecompositionOfWLine
 265 //=======================================================================
 266 static
 267 gp_Pnt2d AdjustByNeighbour(const gp_Pnt2d&     theaNeighbourPoint,
 268                            const gp_Pnt2d&     theOriginalPoint,
 269                            Handle(GeomAdaptor_Surface) theGASurface)
 270 {
 271   gp_Pnt2d ap1 = theaNeighbourPoint;
 272   gp_Pnt2d ap2 = theOriginalPoint;
 273
 274   if ( theGASurface->IsUPeriodic() ) {
 275     Standard_Real aPeriod     = theGASurface->UPeriod();
 276     gp_Pnt2d aPTest = ap2;
 277     Standard_Real aSqDistMin = 1.e+100;
 278
 279     for ( Standard_Integer pIt = -1; pIt <= 1; pIt++) {
 280       aPTest.SetX( theOriginalPoint.X() + aPeriod * pIt );
 281       Standard_Real dd = ap1.SquareDistance( aPTest );
 282
 283       if ( dd < aSqDistMin ) {
 284         ap2 = aPTest;
 285         aSqDistMin = dd;
 286       }
 287     }
 288   }
 289   if ( theGASurface->IsVPeriodic() ) {
 290     Standard_Real aPeriod     = theGASurface->VPeriod();
 291     gp_Pnt2d aPTest = ap2;
 292     Standard_Real aSqDistMin = 1.e+100;
 293
 294     for ( Standard_Integer pIt = -1; pIt <= 1; pIt++) {
 295       aPTest.SetY( theOriginalPoint.Y() + aPeriod * pIt );
 296       Standard_Real dd = ap1.SquareDistance( aPTest );
 297
 298       if ( dd < aSqDistMin ) {
 299         ap2 = aPTest;
 300         aSqDistMin = dd;
 301       }
 302     }
 303   }
 304   return ap2;
 305 }
 306
 307 //=======================================================================
 308 //function : RefineVector
 309 //purpose  : static subfunction in FindPoint
 310 //=======================================================================
 311 static
 312 void RefineVector(gp_Vec2d& aV2D)
 313 {
 314   Standard_Integer k,m;
 315   Standard_Real aC[2], aEps, aR1, aR2, aNum;
 316   //
 317   aEps=RealEpsilon();
 318   aR1=1.-aEps;
 319   aR2=1.+aEps;
 320   //
 321   aV2D.Coord(aC[0], aC[1]);
 322   //
 323   for (k=0; k<2; ++k) {
 324     m=(k+1)%2;
 325     aNum=fabs(aC[k]);
 326     if (aNum>aR1 && aNum<aR2) {
 327       if (aC[k]<0.) {
 328         aC[k]=-1.;
 329       }
 330       else {
 331         aC[k]=1.;
 332       }
 333       aC[m]=0.;
 334       break;
 335     }
 336   }
 337   aV2D.SetCoord(aC[0], aC[1]);
 338 }
 339
 340 //=======================================================================
 341 //function : FindPoint
 342 //purpose  : static subfunction in DecompositionOfWLine
 343 //=======================================================================
 344 static
 345 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
 346                            const gp_Pnt2d&     theLastPoint,
 347                            const Standard_Real theUmin,
 348                            const Standard_Real theUmax,
 349                            const Standard_Real theVmin,
 350                            const Standard_Real theVmax,
 351                            gp_Pnt2d&           theNewPoint)
 352 {
 353   gp_Vec2d aVec(theFirstPoint, theLastPoint);
 354   Standard_Integer i = 0, j = 0;
 355
 356   for(i = 0; i < 4; i++) {
 357     gp_Vec2d anOtherVec;
 358     gp_Vec2d anOtherVecNormal;
 359     gp_Pnt2d aprojpoint = theLastPoint;
 360
 361     if((i % 2) == 0) {
 362       anOtherVec.SetX(0.);
 363       anOtherVec.SetY(1.);
 364       anOtherVecNormal.SetX(1.);
 365       anOtherVecNormal.SetY(0.);
 366
 367       if(i < 2)
 368         aprojpoint.SetX(theUmin);
 369       else
 370         aprojpoint.SetX(theUmax);
 371     }
 372     else {
 373       anOtherVec.SetX(1.);
 374       anOtherVec.SetY(0.);
 375       anOtherVecNormal.SetX(0.);
 376       anOtherVecNormal.SetY(1.);
 377
 378       if(i < 2)
 379         aprojpoint.SetY(theVmin);
 380       else
 381         aprojpoint.SetY(theVmax);
 382     }
 383     gp_Vec2d anormvec = aVec;
 384     anormvec.Normalize();
 385     RefineVector(anormvec);
 386     Standard_Real adot1 = anormvec.Dot(anOtherVecNormal);
 387
 388     if(fabs(adot1) < Precision::Angular())
 389       continue;
 390     Standard_Real adist = 0.;
 391     Standard_Boolean bIsOut = Standard_False;
 392
 393     if((i % 2) == 0) {
 394       adist = (i < 2) ? fabs(theLastPoint.X() - theUmin) : fabs(theLastPoint.X() - theUmax);
 395       bIsOut = (i < 2) ? (theLastPoint.X() < theUmin) : (theLastPoint.X() > theUmax);
 396     }
 397     else {
 398       adist = (i < 2) ? fabs(theLastPoint.Y() - theVmin) : fabs(theLastPoint.Y() - theVmax);
 399       bIsOut = (i < 2) ? (theLastPoint.Y() < theVmin) : (theLastPoint.Y() > theVmax);
 400     }
 401     Standard_Real anoffset = adist * anOtherVec.Dot(anormvec) / adot1;
 402
 403     for(j = 0; j < 2; j++) {
 404       anoffset = (j == 0) ? anoffset : -anoffset;
 405       gp_Pnt2d acurpoint(aprojpoint.XY() + (anOtherVec.XY()*anoffset));
 406       gp_Vec2d acurvec(theLastPoint, acurpoint);
 407       if ( bIsOut )
 408         acurvec.Reverse();
 409
 410       Standard_Real aDotX, anAngleX;
 411       //
 412       aDotX = aVec.Dot(acurvec);
 413       anAngleX = aVec.Angle(acurvec);
 414       //
 415       if(aDotX > 0. && fabs(anAngleX) < Precision::PConfusion()) {
 416         if((i % 2) == 0) {
 417           if((acurpoint.Y() >= theVmin) &&
 418              (acurpoint.Y() <= theVmax)) {
 419             theNewPoint = acurpoint;
 420             return Standard_True;
 421           }
 422         }
 423         else {
 424           if((acurpoint.X() >= theUmin) &&
 425              (acurpoint.X() <= theUmax)) {
 426             theNewPoint = acurpoint;
 427             return Standard_True;
 428           }
 429         }
 430       }
 431     }
 432   }
 433   return Standard_False;
 434 }
 435
 436 //=======================================================================
 437 //function : DecompositionOfWLine
 438 //purpose  :
 439 //=======================================================================
 440 Standard_Boolean IntTools_WLineTool::
 441   DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
 442                        const Handle(GeomAdaptor_Surface)&            theSurface1,
 443                        const Handle(GeomAdaptor_Surface)&            theSurface2,
 444                        const TopoDS_Face&                             theFace1,
 445                        const TopoDS_Face&                             theFace2,
 446                        const GeomInt_LineConstructor&                 theLConstructor,
 447                        const Standard_Boolean                         theAvoidLConstructor,
 448                        const Standard_Real                            theTol,
 449                        IntPatch_SequenceOfLine&                       theNewLines,
 450                        const Handle(IntTools_Context)& aContext)
 451 {
 452   Standard_Boolean bRet, bAvoidLineConstructor;
 453   Standard_Integer aNbPnts, aNbParts;
 454   //
 455   bRet=Standard_False;
 456   aNbPnts=theWLine->NbPnts();
 457   bAvoidLineConstructor=theAvoidLConstructor;
 458   //
 459   if(!aNbPnts){
 460     return bRet;
 461   }
 462   if (!bAvoidLineConstructor) {
 463     aNbParts=theLConstructor.NbParts();
 464     if (!aNbParts) {
 465       return bRet;
 466     }
 467   }
 468   //
 469   Standard_Boolean bIsPrevPointOnBoundary, bIsPointOnBoundary, bIsCurrentPointOnBoundary;
 470   Standard_Integer nblines, pit, i, j;
 471   Standard_Real aTol;
 472   TColStd_Array1OfListOfInteger anArrayOfLines(1, aNbPnts);
 473   TColStd_Array1OfInteger       anArrayOfLineType(1, aNbPnts);
 474   TColStd_ListOfInteger aListOfPointIndex;
 475
 476   //
 477   nblines=0;
 478   aTol=Precision::Confusion();
 479   aTol=0.5*aTol;
 480   bIsPrevPointOnBoundary=Standard_False;
 481   bIsPointOnBoundary=Standard_False;
 482   //
 483   // 1. ...
 484   //
 485   // Points
 486   for(pit = 1; pit <= aNbPnts; ++pit) {
 487     Standard_Boolean bIsOnFirstBoundary, isperiodic;
 488     Standard_Real aResolution, aPeriod, alowerboundary, aupperboundary, U, V;
 489     Standard_Real aParameter, anoffset, anAdjustPar;
 490     Standard_Real umin, umax, vmin, vmax;
 491     //
 492     bIsCurrentPointOnBoundary = Standard_False;
 493     const IntSurf_PntOn2S& aPoint = theWLine->Point(pit);
 494     //
 495     // Surface
 496     for(i = 0; i < 2; ++i) {
 497       Handle(GeomAdaptor_Surface) aGASurface = (!i) ? theSurface1 : theSurface2;
 498       aGASurface->Surface()->Bounds(umin, umax, vmin, vmax);
 499       if(!i) {
 500         aPoint.ParametersOnS1(U, V);
 501       }
 502       else {
 503         aPoint.ParametersOnS2(U, V);
 504       }
 505       // U, V
 506       for(j = 0; j < 2; j++) {
 507         isperiodic = (!j) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
 508         if(!isperiodic){
 509           continue;
 510         }
 511         //
 512         if (!j) {
 513           aResolution=aGASurface->UResolution(aTol);
 514           aPeriod=aGASurface->UPeriod();
 515           alowerboundary=umin;
 516           aupperboundary=umax;
 517           aParameter=U;
 518         }
 519         else {
 520           aResolution=aGASurface->VResolution(aTol);
 521           aPeriod=aGASurface->VPeriod();
 522           alowerboundary=vmin;
 523           aupperboundary=vmax;
 524           aParameter=V;
 525         }
 526
 527         GeomInt::AdjustPeriodic(aParameter,
 528                                        alowerboundary,
 529                                        aupperboundary,
 530                                        aPeriod,
 531                                        anAdjustPar,
 532                                        anoffset);
 533         //
 534         bIsOnFirstBoundary = Standard_True;// ?
 535         bIsPointOnBoundary=
 536           IsPointOnBoundary(anAdjustPar,
 537                             alowerboundary,
 538                             aupperboundary,
 539                             aResolution,
 540                             bIsOnFirstBoundary);
 541         //
 542         if(bIsPointOnBoundary) {
 543           bIsCurrentPointOnBoundary = Standard_True;
 544           break;
 545         }
 546       }//for(j = 0; j < 2; j++) {
 547
 548       if(bIsCurrentPointOnBoundary){
 549         break;
 550       }
 551     }//for(i = 0; i < 2; ++i) {
 552     //
 553     if((bIsCurrentPointOnBoundary != bIsPrevPointOnBoundary)) {
 554       if(!aListOfPointIndex.IsEmpty()) {
 555         nblines++;
 556         anArrayOfLines.SetValue(nblines, aListOfPointIndex);
 557         anArrayOfLineType.SetValue(nblines, bIsPrevPointOnBoundary);
 558         aListOfPointIndex.Clear();
 559       }
 560       bIsPrevPointOnBoundary = bIsCurrentPointOnBoundary;
 561     }
 562     aListOfPointIndex.Append(pit);
 563   } //for(pit = 1; pit <= aNbPnts; ++pit) {
 564   //
 565   if(!aListOfPointIndex.IsEmpty()) {
 566     nblines++;
 567     anArrayOfLines.SetValue(nblines, aListOfPointIndex);
 568     anArrayOfLineType.SetValue(nblines, bIsPrevPointOnBoundary);
 569     aListOfPointIndex.Clear();
 570   }
 571   //
 572   if(nblines<=1) {
 573     return bRet; //Standard_False;
 574   }
 575   //
 576   //
 577   // 2. Correct wlines.begin
 578   TColStd_Array1OfListOfInteger anArrayOfLineEnds(1, nblines);
 579   Handle(IntSurf_LineOn2S) aSeqOfPntOn2S = new IntSurf_LineOn2S();
 580   //
 581   for(i = 1; i <= nblines; i++) {
 582     if(anArrayOfLineType.Value(i) != 0) {
 583       continue;
 584     }
 585     const TColStd_ListOfInteger& aListOfIndex = anArrayOfLines.Value(i);
 586     TColStd_ListOfInteger aListOfFLIndex;
 587
 588     for(j = 0; j < 2; j++) {
 589       Standard_Integer aneighbourindex = (j == 0) ? (i - 1) : (i + 1);
 590
 591       if((aneighbourindex < 1) || (aneighbourindex > nblines))
 592         continue;
 593
 594       if(anArrayOfLineType.Value(aneighbourindex) == 0)
 595         continue;
 596       const TColStd_ListOfInteger& aNeighbour = anArrayOfLines.Value(aneighbourindex);
 597       Standard_Integer anIndex = (j == 0) ? aNeighbour.Last() : aNeighbour.First();
 598       const IntSurf_PntOn2S& aPoint = theWLine->Point(anIndex);
 599
 600       IntSurf_PntOn2S aNewP = aPoint;
 601       if(aListOfIndex.Extent() < 2) {
 602         aSeqOfPntOn2S->Add(aNewP);
 603         aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
 604         continue;
 605       }
 606       //
 607       Standard_Integer iFirst = aListOfIndex.First();
 608       Standard_Integer iLast  = aListOfIndex.Last();
 609       //
 610       for(Standard_Integer surfit = 0; surfit < 2; surfit++) {
 611
 612         Handle(GeomAdaptor_Surface) aGASurface = (surfit == 0) ? theSurface1 : theSurface2;
 613         Standard_Real umin=0., umax=0., vmin=0., vmax=0.;
 614         aGASurface->Surface()->Bounds(umin, umax, vmin, vmax);
 615         Standard_Real U=0., V=0.;
 616
 617         if(surfit == 0)
 618           aNewP.ParametersOnS1(U, V);
 619         else
 620           aNewP.ParametersOnS2(U, V);
 621         Standard_Integer nbboundaries = 0;
 622
 623         Standard_Boolean bIsNearBoundary = Standard_False;
 624         //Standard_Integer aZoneIndex = 0;
 625         Standard_Integer bIsUBoundary = Standard_False; // use if nbboundaries == 1
 626         Standard_Integer bIsFirstBoundary = Standard_False; // use if nbboundaries == 1
 627
 628
 629         for(Standard_Integer parit = 0; parit < 2; parit++) {
 630           Standard_Boolean isperiodic = (parit == 0) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
 631
 632           Standard_Real aResolution = (parit == 0) ? aGASurface->UResolution(aTol) : aGASurface->VResolution(aTol);
 633           Standard_Real alowerboundary = (parit == 0) ? umin : vmin;
 634           Standard_Real aupperboundary = (parit == 0) ? umax : vmax;
 635
 636           Standard_Real aParameter = (parit == 0) ? U : V;
 637           Standard_Boolean bIsOnFirstBoundary = Standard_True;
 638
 639           if(!isperiodic) {
 640             bIsPointOnBoundary=
 641               IsPointOnBoundary(aParameter, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
 642             if(bIsPointOnBoundary) {
 643               bIsUBoundary = (parit == 0);
 644               bIsFirstBoundary = bIsOnFirstBoundary;
 645               nbboundaries++;
 646             }
 647           }
 648           else {
 649             Standard_Real aPeriod     = (parit == 0) ? aGASurface->UPeriod() : aGASurface->VPeriod();
 650             Standard_Real anoffset, anAdjustPar;
 651             GeomInt::AdjustPeriodic(aParameter, alowerboundary, aupperboundary,
 652                                            aPeriod, anAdjustPar, anoffset);
 653
 654             bIsPointOnBoundary=
 655               IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
 656             if(bIsPointOnBoundary) {
 657               bIsUBoundary = (parit == 0);
 658               bIsFirstBoundary = bIsOnFirstBoundary;
 659               nbboundaries++;
 660             }
 661             else {
 662             //check neighbourhood of boundary
 663             Standard_Real anEpsilon = aResolution * 100.;
 664             Standard_Real aPart = ( aupperboundary - alowerboundary ) * 0.1;
 665             anEpsilon = ( anEpsilon > aPart ) ? aPart : anEpsilon;
 666
 667               bIsNearBoundary = IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary,
 668                                                   anEpsilon, bIsOnFirstBoundary);
 669
 670             }
 671           }
 672         }
 673
 674         Standard_Boolean bComputeLineEnd = Standard_False;
 675
 676         if(nbboundaries == 2) {
 677           //xf
 678           bComputeLineEnd = Standard_True;
 679           //xt
 680         }
 681         else if(nbboundaries == 1) {
 682           Standard_Boolean isperiodic = (bIsUBoundary) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
 683
 684           if(isperiodic) {
 685             Standard_Real alowerboundary = (bIsUBoundary) ? umin : vmin;
 686             Standard_Real aupperboundary = (bIsUBoundary) ? umax : vmax;
 687             Standard_Real aPeriod     = (bIsUBoundary) ? aGASurface->UPeriod() : aGASurface->VPeriod();
 688             Standard_Real aParameter = (bIsUBoundary) ? U : V;
 689             Standard_Real anoffset, anAdjustPar;
 690             GeomInt::AdjustPeriodic(aParameter, alowerboundary, aupperboundary,
 691                                            aPeriod, anAdjustPar, anoffset);
 692
 693             Standard_Real adist = (bIsFirstBoundary) ? fabs(anAdjustPar - alowerboundary) : fabs(anAdjustPar - aupperboundary);
 694             Standard_Real anotherPar = (bIsFirstBoundary) ? (aupperboundary - adist) : (alowerboundary + adist);
 695             anotherPar += anoffset;
 696             Standard_Integer aneighbourpointindex = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
 697             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex);
 698             Standard_Real nU1, nV1;
 699
 700             if(surfit == 0)
 701               aNeighbourPoint.ParametersOnS1(nU1, nV1);
 702             else
 703               aNeighbourPoint.ParametersOnS2(nU1, nV1);
 704
 705             Standard_Real adist1 = (bIsUBoundary) ? fabs(nU1 - U) : fabs(nV1 - V);
 706             Standard_Real adist2 = (bIsUBoundary) ? fabs(nU1 - anotherPar) : fabs(nV1 - anotherPar);
 707             bComputeLineEnd = Standard_True;
 708             Standard_Boolean bCheckAngle1 = Standard_False;
 709             Standard_Boolean bCheckAngle2 = Standard_False;
 710             gp_Vec2d aNewVec;
 711             Standard_Real anewU = (bIsUBoundary) ? anotherPar : U;
 712             Standard_Real anewV = (bIsUBoundary) ? V : anotherPar;
 713
 714             if(((adist1 - adist2) > Precision::PConfusion()) &&
 715                (adist2 < (aPeriod / 4.))) {
 716               bCheckAngle1 = Standard_True;
 717               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(anewU, anewV));
 718
 719               if(aNewVec.SquareMagnitude() < gp::Resolution()) {
 720                 aNewP.SetValue((surfit == 0), anewU, anewV);
 721                 bCheckAngle1 = Standard_False;
 722               }
 723             }
 724             else if(adist1 < (aPeriod / 4.)) {
 725               bCheckAngle2 = Standard_True;
 726               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(U, V));
 727
 728               if(aNewVec.SquareMagnitude() < gp::Resolution()) {
 729                 bCheckAngle2 = Standard_False;
 730               }
 731             }
 732
 733             if(bCheckAngle1 || bCheckAngle2) {
 734               // assume there are at least two points in line (see "if" above)
 735               Standard_Integer anindexother = aneighbourpointindex;
 736
 737               while((anindexother <= iLast) && (anindexother >= iFirst)) {
 738                 anindexother = (j == 0) ? (anindexother + 1) : (anindexother - 1);
 739                 const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(anindexother);
 740                 Standard_Real nU2, nV2;
 741
 742                 if(surfit == 0)
 743                   aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
 744                 else
 745                   aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
 746                 gp_Vec2d aVecOld(gp_Pnt2d(nU2, nV2), gp_Pnt2d(nU1, nV1));
 747
 748                 if(aVecOld.SquareMagnitude() <= gp::Resolution()) {
 749                   continue;
 750                 }
 751                 else {
 752                   Standard_Real anAngle = aNewVec.Angle(aVecOld);
 753
 754                   if((fabs(anAngle) < (M_PI * 0.25)) && (aNewVec.Dot(aVecOld) > 0.)) {
 755
 756                     if(bCheckAngle1) {
 757                       Standard_Real U1, U2, V1, V2;
 758                       IntSurf_PntOn2S atmppoint = aNewP;
 759                       atmppoint.SetValue((surfit == 0), anewU, anewV);
 760                       atmppoint.Parameters(U1, V1, U2, V2);
 761                       gp_Pnt P1 = theSurface1->Value(U1, V1);
 762                       gp_Pnt P2 = theSurface2->Value(U2, V2);
 763                       gp_Pnt P0 = aPoint.Value();
 764
 765                       if(P0.IsEqual(P1, aTol) &&
 766                          P0.IsEqual(P2, aTol) &&
 767                          P1.IsEqual(P2, aTol)) {
 768                         bComputeLineEnd = Standard_False;
 769                         aNewP.SetValue((surfit == 0), anewU, anewV);
 770                       }
 771                     }
 772
 773                     if(bCheckAngle2) {
 774                       bComputeLineEnd = Standard_False;
 775                     }
 776                   }
 777                   break;
 778                 }
 779               } // end while(anindexother...)
 780             }
 781           }
 782         }
 783         else if ( bIsNearBoundary ) {
 784           bComputeLineEnd = Standard_True;
 785         }
 786
 787         if(bComputeLineEnd) {
 788
 789           gp_Pnt2d anewpoint;
 790           Standard_Boolean found = Standard_False;
 791
 792           if ( bIsNearBoundary ) {
 793             // re-compute point near natural boundary or near tangent zone
 794             Standard_Real u1, v1, u2, v2;
 795             aNewP.Parameters( u1, v1, u2, v2 );
 796             if(surfit == 0)
 797               anewpoint = gp_Pnt2d( u1, v1 );
 798             else
 799               anewpoint = gp_Pnt2d( u2, v2 );
 800
 801             Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
 802             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
 803             Standard_Real nU1, nV1;
 804
 805             if(surfit == 0)
 806               aNeighbourPoint.ParametersOnS1(nU1, nV1);
 807             else
 808               aNeighbourPoint.ParametersOnS2(nU1, nV1);
 809             gp_Pnt2d ap1(nU1, nV1);
 810             gp_Pnt2d ap2;
 811
 812             if ( aGASurface->IsUPeriodic() || aGASurface->IsVPeriodic() ) {
 813               // re-compute point near boundary if shifted on a period
 814               ap2 = AdjustByNeighbour( ap1, anewpoint, aGASurface );
 815
 816               if ( ( ap2.X() < umin ) || ( ap2.X() > umax ) ||
 817                   ( ap2.Y() < vmin ) || ( ap2.Y() > vmax ) ) {
 818                 found = FindPoint(ap1, ap2, umin, umax, vmin, vmax, anewpoint);
 819               }
 820               else {
 821                 anewpoint = ap2;
 822                 aNewP.SetValue( (surfit == 0), anewpoint.X(), anewpoint.Y() );
 823               }
 824             }
 825           }
 826           else {
 827
 828             Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
 829             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
 830             Standard_Real nU1, nV1;
 831
 832             if(surfit == 0)
 833               aNeighbourPoint.ParametersOnS1(nU1, nV1);
 834             else
 835               aNeighbourPoint.ParametersOnS2(nU1, nV1);
 836             gp_Pnt2d ap1(nU1, nV1);
 837             gp_Pnt2d ap2(nU1, nV1);
 838             Standard_Integer aneighbourpointindex2 = aneighbourpointindex1;
 839
 840             while((aneighbourpointindex2 <= iLast) && (aneighbourpointindex2 >= iFirst)) {
 841               aneighbourpointindex2 = (j == 0) ? (aneighbourpointindex2 + 1) : (aneighbourpointindex2 - 1);
 842               const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(aneighbourpointindex2);
 843               Standard_Real nU2, nV2;
 844
 845               if(surfit == 0)
 846                 aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
 847               else
 848                 aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
 849               ap2.SetX(nU2);
 850               ap2.SetY(nV2);
 851
 852               if(ap1.SquareDistance(ap2) > gp::Resolution()) {
 853                 break;
 854               }
 855             }
 856             found = FindPoint(ap2, ap1, umin, umax, vmin, vmax, anewpoint);
 857           }
 858
 859           if(found) {
 860             // check point
 861             Standard_Real aCriteria = theTol;
 862             GeomAPI_ProjectPointOnSurf& aProjector =
 863               (surfit == 0) ? aContext->ProjPS(theFace2) : aContext->ProjPS(theFace1);
 864             Handle(GeomAdaptor_Surface) aSurface = (surfit == 0) ? theSurface1 : theSurface2;
 865
 866             Handle(GeomAdaptor_Surface) aSurfaceOther = (surfit == 0) ? theSurface2 : theSurface1;
 867
 868             gp_Pnt aP3d = aSurface->Value(anewpoint.X(), anewpoint.Y());
 869             aProjector.Perform(aP3d);
 870
 871             if(aProjector.IsDone()) {
 872               if(aProjector.LowerDistance() < aCriteria) {
 873                 Standard_Real foundU = U, foundV = V;
 874                 aProjector.LowerDistanceParameters(foundU, foundV);
 875
 876                 //Correction of projected coordinates. Begin
 877                 //Note, it may be shifted on a period
 878                 Standard_Integer aneindex1 = (j == 0) ? iFirst : iLast;
 879                 const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneindex1);
 880                 Standard_Real nUn, nVn;
 881
 882                 if(surfit == 0)
 883                   aNeighbourPoint.ParametersOnS2(nUn, nVn);
 884                 else
 885                   aNeighbourPoint.ParametersOnS1(nUn, nVn);
 886                 gp_Pnt2d aNeighbour2d(nUn, nVn);
 887                 gp_Pnt2d anAdjustedPoint = AdjustByNeighbour( aNeighbour2d, gp_Pnt2d(foundU, foundV), aSurfaceOther );
 888                 foundU = anAdjustedPoint.X();
 889                 foundV = anAdjustedPoint.Y();
 890
 891                 if ( ( anAdjustedPoint.X() < umin ) && ( anAdjustedPoint.X() > umax ) &&
 892                     ( anAdjustedPoint.Y() < vmin ) && ( anAdjustedPoint.Y() > vmax ) ) {
 893                   // attempt to roughly re-compute point
 894                   foundU = ( foundU < umin ) ? umin : foundU;
 895                   foundU = ( foundU > umax ) ? umax : foundU;
 896                   foundV = ( foundV < vmin ) ? vmin : foundV;
 897                   foundV = ( foundV > vmax ) ? vmax : foundV;
 898
 899                   GeomAPI_ProjectPointOnSurf& aProjector2 =
 900                     (surfit == 0) ? aContext->ProjPS(theFace1) : aContext->ProjPS(theFace2);
 901
 902                   aP3d = aSurfaceOther->Value(foundU, foundV);
 903                   aProjector2.Perform(aP3d);
 904
 905                   if(aProjector2.IsDone()) {
 906                     if(aProjector2.LowerDistance() < aCriteria) {
 907                       Standard_Real foundU2 = anewpoint.X(), foundV2 = anewpoint.Y();
 908                       aProjector2.LowerDistanceParameters(foundU2, foundV2);
 909                       anewpoint.SetX(foundU2);
 910                       anewpoint.SetY(foundV2);
 911                     }
 912                   }
 913                 }
 914                 //Correction of projected coordinates. End
 915
 916                 if(surfit == 0)
 917                   aNewP.SetValue(aP3d, anewpoint.X(), anewpoint.Y(), foundU, foundV);
 918                 else
 919                   aNewP.SetValue(aP3d, foundU, foundV, anewpoint.X(), anewpoint.Y());
 920               }
 921             }
 922           }
 923         }
 924       }
 925       aSeqOfPntOn2S->Add(aNewP);
 926       aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
 927     }
 928     anArrayOfLineEnds.SetValue(i, aListOfFLIndex);
 929   }
 930   // Correct wlines.end
 931
 932   // Split wlines.begin
 933   Standard_Integer nbiter;
 934   //
 935   nbiter=1;
 936   if (!bAvoidLineConstructor) {
 937     nbiter=theLConstructor.NbParts();
 938   }
 939   //
 940   for(j = 1; j <= nbiter; ++j) {
 941     Standard_Real fprm, lprm;
 942     Standard_Integer ifprm, ilprm;
 943     //
 944     if(bAvoidLineConstructor) {
 945       ifprm = 1;
 946       ilprm = theWLine->NbPnts();
 947     }
 948     else {
 949       theLConstructor.Part(j, fprm, lprm);
 950       ifprm = (Standard_Integer)fprm;
 951       ilprm = (Standard_Integer)lprm;
 952     }
 953
 954     Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
 955     //
 956     for(i = 1; i <= nblines; i++) {
 957       if(anArrayOfLineType.Value(i) != 0) {
 958         continue;
 959       }
 960       const TColStd_ListOfInteger& aListOfIndex = anArrayOfLines.Value(i);
 961       const TColStd_ListOfInteger& aListOfFLIndex = anArrayOfLineEnds.Value(i);
 962       Standard_Boolean bhasfirstpoint = (aListOfFLIndex.Extent() == 2);
 963       Standard_Boolean bhaslastpoint = (aListOfFLIndex.Extent() == 2);
 964
 965       if(!bhasfirstpoint && !aListOfFLIndex.IsEmpty()) {
 966         bhasfirstpoint = (i != 1);
 967       }
 968
 969       if(!bhaslastpoint && !aListOfFLIndex.IsEmpty()) {
 970         bhaslastpoint = (i != nblines);
 971       }
 972
 973       Standard_Integer iFirst = aListOfIndex.First();
 974       Standard_Integer iLast  = aListOfIndex.Last();
 975       Standard_Boolean bIsFirstInside = ((ifprm >= iFirst) && (ifprm <= iLast));
 976       Standard_Boolean bIsLastInside =  ((ilprm >= iFirst) && (ilprm <= iLast));
 977
 978       if(!bIsFirstInside && !bIsLastInside) {
 979         if((ifprm < iFirst) && (ilprm > iLast)) {
 980           // append whole line, and boundaries if necessary
 981           if(bhasfirstpoint) {
 982             pit = aListOfFLIndex.First();
 983             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
 984             aLineOn2S->Add(aP);
 985           }
 986           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
 987
 988           for(; anIt.More(); anIt.Next()) {
 989             pit = anIt.Value();
 990             const IntSurf_PntOn2S& aP = theWLine->Point(pit);
 991             aLineOn2S->Add(aP);
 992           }
 993
 994           if(bhaslastpoint) {
 995             pit = aListOfFLIndex.Last();
 996             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
 997             aLineOn2S->Add(aP);
 998           }
 999
1000           // check end of split line (end is almost always)
1001           Standard_Integer aneighbour = i + 1;
1002           Standard_Boolean bIsEndOfLine = Standard_True;
1003
1004           if(aneighbour <= nblines) {
1005             const TColStd_ListOfInteger& aListOfNeighbourIndex = anArrayOfLines.Value(aneighbour);
1006
1007             if((anArrayOfLineType.Value(aneighbour) != 0) &&
1008                (aListOfNeighbourIndex.IsEmpty())) {
1009               bIsEndOfLine = Standard_False;
1010             }
1011           }
1012
1013           if(bIsEndOfLine) {
1014             if(aLineOn2S->NbPoints() > 1) {
1015               Handle(IntPatch_WLine) aNewWLine =
1016                 new IntPatch_WLine(aLineOn2S, Standard_False);
1017               aNewWLine->SetCreatingWayInfo(theWLine->GetCreatingWay());
1018               theNewLines.Append(aNewWLine);
1019             }
1020             aLineOn2S = new IntSurf_LineOn2S();
1021           }
1022         }
1023         continue;
1024       }
1025       // end if(!bIsFirstInside && !bIsLastInside)
1026
1027       if(bIsFirstInside && bIsLastInside) {
1028         // append inside points between ifprm and ilprm
1029         TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
1030
1031         for(; anIt.More(); anIt.Next()) {
1032           pit = anIt.Value();
1033           if((pit < ifprm) || (pit > ilprm))
1034             continue;
1035           const IntSurf_PntOn2S& aP = theWLine->Point(pit);
1036           aLineOn2S->Add(aP);
1037         }
1038       }
1039       else {
1040
1041         if(bIsFirstInside) {
1042           // append points from ifprm to last point + boundary point
1043           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
1044
1045           for(; anIt.More(); anIt.Next()) {
1046             pit = anIt.Value();
1047             if(pit < ifprm)
1048               continue;
1049             const IntSurf_PntOn2S& aP = theWLine->Point(pit);
1050             aLineOn2S->Add(aP);
1051           }
1052
1053           if(bhaslastpoint) {
1054             pit = aListOfFLIndex.Last();
1055             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
1056             aLineOn2S->Add(aP);
1057           }
1058           // check end of split line (end is almost always)
1059           Standard_Integer aneighbour = i + 1;
1060           Standard_Boolean bIsEndOfLine = Standard_True;
1061
1062           if(aneighbour <= nblines) {
1063             const TColStd_ListOfInteger& aListOfNeighbourIndex = anArrayOfLines.Value(aneighbour);
1064
1065             if((anArrayOfLineType.Value(aneighbour) != 0) &&
1066                (aListOfNeighbourIndex.IsEmpty())) {
1067               bIsEndOfLine = Standard_False;
1068             }
1069           }
1070
1071           if(bIsEndOfLine) {
1072             if(aLineOn2S->NbPoints() > 1) {
1073               Handle(IntPatch_WLine) aNewWLine =
1074                 new IntPatch_WLine(aLineOn2S, Standard_False);
1075               aNewWLine->SetCreatingWayInfo(theWLine->GetCreatingWay());
1076               theNewLines.Append(aNewWLine);
1077             }
1078             aLineOn2S = new IntSurf_LineOn2S();
1079           }
1080         }
1081         // end if(bIsFirstInside)
1082
1083         if(bIsLastInside) {
1084           // append points from first boundary point to ilprm
1085           if(bhasfirstpoint) {
1086             pit = aListOfFLIndex.First();
1087             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
1088             aLineOn2S->Add(aP);
1089           }
1090           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
1091
1092           for(; anIt.More(); anIt.Next()) {
1093             pit = anIt.Value();
1094             if(pit > ilprm)
1095               continue;
1096             const IntSurf_PntOn2S& aP = theWLine->Point(pit);
1097             aLineOn2S->Add(aP);
1098           }
1099         }
1100         //end if(bIsLastInside)
1101       }
1102     }
1103
1104     if(aLineOn2S->NbPoints() > 1) {
1105       Handle(IntPatch_WLine) aNewWLine =
1106         new IntPatch_WLine(aLineOn2S, Standard_False);
1107       aNewWLine->SetCreatingWayInfo(theWLine->GetCreatingWay());
1108       theNewLines.Append(aNewWLine);
1109     }
1110   }
1111   // Split wlines.end
1112
1113   return Standard_True;
1114 }
1115
1116 ///////////////////// end of DecompositionOfWLine ///////////////////////