src/GeomInt/GeomInt_IntSS_1.cxx

   1 // Created on: 1995-01-27
   2 // Created by: Jacques GOUSSARD
   3 // Copyright (c) 1995-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 #include <algorithm>
  18
  19 #include <GeomInt_IntSS.ixx>
  20
  21 #include <GeomInt.hxx>
  22
  23 #include <Precision.hxx>
  24 #include <gp_Pnt.hxx>
  25 #include <gp_Pnt2d.hxx>
  26 #include <gp_Pln.hxx>
  27
  28 #include <TColStd_Array1OfReal.hxx>
  29 #include <TColStd_Array1OfInteger.hxx>
  30 #include <TColStd_Array1OfListOfInteger.hxx>
  31 #include <TColStd_SequenceOfReal.hxx>
  32 #include <TColStd_ListOfInteger.hxx>
  33 #include <TColStd_ListIteratorOfListOfInteger.hxx>
  34
  35 #include <TColgp_Array1OfPnt.hxx>
  36 #include <TColgp_Array1OfPnt2d.hxx>
  37
  38 #include <Adaptor3d_TopolTool.hxx>
  39 #include <IntPatch_Line.hxx>
  40 #include <IntPatch_WLine.hxx>
  41 #include <IntPatch_GLine.hxx>
  42 #include <IntPatch_RLine.hxx>
  43 #include <IntPatch_ALineToWLine.hxx>
  44 #include <IntPatch_IType.hxx>
  45 #include <NCollection_IncAllocator.hxx>
  46 #include <NCollection_List.hxx>
  47 #include <NCollection_LocalArray.hxx>
  48 #include <NCollection_StdAllocator.hxx>
  49 #include <vector>
  50
  51 #include <AppParCurves_MultiBSpCurve.hxx>
  52
  53 #include <GeomAbs_SurfaceType.hxx>
  54 #include <GeomAbs_CurveType.hxx>
  55
  56 #include <GeomAdaptor.hxx>
  57 #include <Geom_Curve.hxx>
  58 #include <Geom_Line.hxx>
  59 #include <Geom_Parabola.hxx>
  60 #include <Geom_Hyperbola.hxx>
  61 #include <Geom_TrimmedCurve.hxx>
  62 #include <Geom_Circle.hxx>
  63 #include <Geom_Ellipse.hxx>
  64 #include <Geom_BSplineCurve.hxx>
  65
  66 #include <Geom2dAdaptor.hxx>
  67 #include <Adaptor2d_HCurve2d.hxx>
  68 #include <Geom2d_Curve.hxx>
  69 #include <Geom2d_BSplineCurve.hxx>
  70 #include <Geom2d_TrimmedCurve.hxx>
  71
  72 #include <GeomLib_CheckBSplineCurve.hxx>
  73 #include <GeomLib_Check2dBSplineCurve.hxx>
  74 #include <GeomProjLib.hxx>
  75 #include <Approx_CurveOnSurface.hxx>
  76
  77 #include <ElSLib.hxx>
  78
  79 #include <GeomInt_WLApprox.hxx>
  80 #include <Extrema_ExtPS.hxx>
  81
  82 #include <GeomAdaptor_HSurface.hxx>
  83 #include <gp_Lin2d.hxx>
  84 #include <Geom2d_Line.hxx>
  85 #include <IntPatch_RLine.hxx>
  86 #include <Geom2dAdaptor.hxx>
  87 #include <Adaptor2d_HCurve2d.hxx>
  88 #include <Approx_CurveOnSurface.hxx>
  89 #include <GeomAdaptor.hxx>
  90 #include <GeomProjLib.hxx>
  91 #include <TColStd_Array1OfReal.hxx>
  92 #include <TColgp_Array1OfPnt2d.hxx>
  93 #include <TColStd_Array1OfInteger.hxx>
  94 #include <Geom2d_BSplineCurve.hxx>
  95 #include <Geom2dAdaptor_Curve.hxx>
  96 #include <IntRes2d_IntersectionPoint.hxx>
  97 #include <IntRes2d_IntersectionSegment.hxx>
  98 #include <Geom2dInt_GInter.hxx>
  99
 100 static
 101   void Parameters(const Handle(GeomAdaptor_HSurface)& HS1,
 102                   const Handle(GeomAdaptor_HSurface)& HS2,
 103                   const gp_Pnt& Ptref,
 104                   Standard_Real& U1,
 105                   Standard_Real& V1,
 106                   Standard_Real& U2,
 107                   Standard_Real& V2);
 108
 109 static
 110   Handle(Geom_Curve) MakeBSpline  (const Handle(IntPatch_WLine)& WL,
 111                                    const Standard_Integer ideb,
 112                                    const Standard_Integer ifin);
 113
 114 static
 115   Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
 116                                             const Standard_Integer                         ideb,
 117                                             const Standard_Integer                         ifin,
 118                                             const Standard_Boolean                         onFirst);
 119
 120 static
 121   Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine);
 122
 123
 124 static
 125   Standard_Boolean DecompositionOfWLine  (const Handle(IntPatch_WLine)& theWLine,
 126                                           const Handle(GeomAdaptor_HSurface)&            theSurface1,
 127                                           const Handle(GeomAdaptor_HSurface)&            theSurface2,
 128                                           const Standard_Real aTolSum,
 129                                           const GeomInt_LineConstructor&                 theLConstructor,
 130                                           IntPatch_SequenceOfLine&                       theNewLines);
 131
 132 static
 133   Standard_Real AdjustPeriodic(const Standard_Real theParameter,
 134                                const Standard_Real parmin,
 135                                const Standard_Real parmax,
 136                                const Standard_Real thePeriod,
 137                                Standard_Real&      theOffset) ;
 138 static
 139   Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
 140                                      const Standard_Real theFirstBoundary,
 141                                      const Standard_Real theSecondBoundary,
 142                                      const Standard_Real theResolution,
 143                                      Standard_Boolean&   IsOnFirstBoundary);
 144
 145 static
 146   Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
 147                              const gp_Pnt2d&     theLastPoint,
 148                              const Standard_Real theUmin,
 149                              const Standard_Real theUmax,
 150                              const Standard_Real theVmin,
 151                              const Standard_Real theVmax,
 152                              gp_Pnt2d&           theNewPoint);
 153
 154
 155 static
 156   void AdjustUPeriodic (const Handle(Geom_Surface)& aS,
 157                         Handle(Geom2d_Curve)& aC2D);
 158 static
 159   void GetQuadric(const Handle(GeomAdaptor_HSurface)& HS1,
 160                   IntSurf_Quadric& quad1);
 161
 162 class ProjectPointOnSurf
 163 {
 164  public:
 165   ProjectPointOnSurf() : myIsDone (Standard_False),myIndex(0) {}
 166   void Init(const Handle(Geom_Surface)& Surface,
 167             const Standard_Real Umin,
 168             const Standard_Real Usup,
 169             const Standard_Real Vmin,
 170             const Standard_Real Vsup);
 171   void Init ();
 172   void Perform(const gp_Pnt& P);
 173   Standard_Boolean IsDone () const { return myIsDone; }
 174   void LowerDistanceParameters(Standard_Real& U, Standard_Real& V ) const;
 175   Standard_Real LowerDistance() const;
 176  protected:
 177   Standard_Boolean myIsDone;
 178   Standard_Integer myIndex;
 179   Extrema_ExtPS myExtPS;
 180   GeomAdaptor_Surface myGeomAdaptor;
 181 };
 182
 183 //=======================================================================
 184 //function : Init
 185 //purpose  :
 186 //=======================================================================
 187   void ProjectPointOnSurf::Init (const Handle(Geom_Surface)& Surface,
 188                                  const Standard_Real       Umin,
 189                                  const Standard_Real       Usup,
 190                                  const Standard_Real       Vmin,
 191                                  const Standard_Real       Vsup )
 192 {
 193   const Standard_Real Tolerance = Precision::PConfusion();
 194   //
 195   myGeomAdaptor.Load(Surface, Umin,Usup,Vmin,Vsup);
 196   myExtPS.Initialize(myGeomAdaptor, Umin, Usup, Vmin, Vsup, Tolerance, Tolerance);
 197   myIsDone = Standard_False;
 198 }
 199 //=======================================================================
 200 //function : Init
 201 //purpose  :
 202 //=======================================================================
 203   void ProjectPointOnSurf::Init ()
 204 {
 205   myIsDone = myExtPS.IsDone() && (myExtPS.NbExt() > 0);
 206   if (myIsDone) {
 207     // evaluate the lower distance and its index;
 208     Standard_Real Dist2Min = myExtPS.SquareDistance(1);
 209     myIndex = 1;
 210     for (Standard_Integer i = 2; i <= myExtPS.NbExt(); i++)
 211     {
 212       const Standard_Real Dist2 = myExtPS.SquareDistance(i);
 213       if (Dist2 < Dist2Min) {
 214         Dist2Min = Dist2;
 215         myIndex = i;
 216       }
 217     }
 218   }
 219 }
 220 //=======================================================================
 221 //function : Perform
 222 //purpose  :
 223 //=======================================================================
 224 void ProjectPointOnSurf::Perform(const gp_Pnt& P)
 225 {
 226   myExtPS.Perform(P);
 227   Init ();
 228 }
 229 //=======================================================================
 230 //function : LowerDistanceParameters
 231 //purpose  :
 232 //=======================================================================
 233 void ProjectPointOnSurf::LowerDistanceParameters (Standard_Real& U,
 234                                                   Standard_Real& V ) const
 235 {
 236   StdFail_NotDone_Raise_if(!myIsDone, "GeomInt_IntSS::ProjectPointOnSurf::LowerDistanceParameters");
 237   (myExtPS.Point(myIndex)).Parameter(U,V);
 238 }
 239 //=======================================================================
 240 //function : LowerDistance
 241 //purpose  :
 242 //=======================================================================
 243 Standard_Real ProjectPointOnSurf::LowerDistance() const
 244 {
 245   StdFail_NotDone_Raise_if(!myIsDone, "GeomInt_IntSS::ProjectPointOnSurf::LowerDistance");
 246   return sqrt(myExtPS.SquareDistance(myIndex));
 247 }
 248 //
 249
 250 //=======================================================================
 251 //function : MakeCurve
 252 //purpose  :
 253 //=======================================================================
 254 void GeomInt_IntSS::MakeCurve(const Standard_Integer Index,
 255                               const Handle(Adaptor3d_TopolTool) & dom1,
 256                               const Handle(Adaptor3d_TopolTool) & dom2,
 257                               const Standard_Real Tol,
 258                               const Standard_Boolean Approx,
 259                               const Standard_Boolean ApproxS1,
 260                               const Standard_Boolean ApproxS2)
 261
 262 {
 263   Standard_Boolean myApprox1, myApprox2, myApprox;
 264   Standard_Real Tolpc, myTolApprox;
 265   IntPatch_IType typl;
 266   Handle(Geom2d_BSplineCurve) H1;
 267   Handle(Geom_Surface) aS1, aS2;
 268   //
 269   Tolpc = Tol;
 270   myApprox=Approx;
 271   myApprox1=ApproxS1;
 272   myApprox2=ApproxS2;
 273   myTolApprox=0.0000001;
 274   //
 275   aS1=myHS1->ChangeSurface().Surface();
 276   aS2=myHS2->ChangeSurface().Surface();
 277   //
 278   Handle(IntPatch_Line) L = myIntersector.Line(Index);
 279   typl = L->ArcType();
 280   //
 281   if(typl==IntPatch_Walking) {
 282     Handle(IntPatch_Line) anewL =
 283       ComputePurgedWLine(Handle(IntPatch_WLine)::DownCast(L));
 284     //
 285     if(anewL.IsNull()) {
 286       return;
 287     }
 288     L = anewL;
 289   }
 290   //
 291   // Line Constructor
 292   myLConstruct.Perform(L);
 293   if (!myLConstruct.IsDone() || myLConstruct.NbParts() <= 0) {
 294     return;
 295   }
 296   // Do the Curve
 297   Standard_Boolean ok;
 298   Standard_Integer i, j,  aNbParts;
 299   Standard_Real fprm, lprm;
 300   Handle(Geom_Curve) newc;
 301
 302   switch (typl) {
 303     //########################################
 304     // Line, Parabola, Hyperbola
 305     //########################################
 306   case IntPatch_Lin:
 307   case IntPatch_Parabola:
 308   case IntPatch_Hyperbola: {
 309     if (typl == IntPatch_Lin) {
 310       newc=new Geom_Line (Handle(IntPatch_GLine)::DownCast(L)->Line());
 311     }
 312     else if (typl == IntPatch_Parabola) {
 313       newc=new Geom_Parabola(Handle(IntPatch_GLine)::DownCast(L)->Parabola());
 314     }
 315     else if (typl == IntPatch_Hyperbola) {
 316       newc=new Geom_Hyperbola (Handle(IntPatch_GLine)::DownCast(L)->Hyperbola());
 317     }
 318     //
 319     aNbParts=myLConstruct.NbParts();
 320     for (i=1; i<=aNbParts; i++) {
 321       myLConstruct.Part(i, fprm, lprm);
 322
 323       if (!Precision::IsNegativeInfinite(fprm) &&
 324         !Precision::IsPositiveInfinite(lprm)) {
 325           Handle(Geom_TrimmedCurve) aCT3D=new Geom_TrimmedCurve(newc, fprm, lprm);
 326           sline.Append(aCT3D);
 327           //
 328           if(myApprox1) {
 329             Handle (Geom2d_Curve) C2d;
 330             BuildPCurves(fprm, lprm, Tolpc, myHS1->ChangeSurface().Surface(), newc, C2d);
 331             if(Tolpc>myTolReached2d || myTolReached2d==0.) {
 332               myTolReached2d=Tolpc;
 333             }
 334             slineS1.Append(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
 335           }
 336           else {
 337             slineS1.Append(H1);
 338           }
 339           //
 340           if(myApprox2) {
 341             Handle (Geom2d_Curve) C2d;
 342             BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
 343             if(Tolpc>myTolReached2d || myTolReached2d==0.) {
 344               myTolReached2d=Tolpc;
 345             }
 346             //
 347             slineS2.Append(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
 348           }
 349           else {
 350             slineS2.Append(H1);
 351           }
 352       } // if (!Precision::IsNegativeInfinite(fprm) &&  !Precision::IsPositiveInfinite(lprm))
 353
 354       else {
 355         GeomAbs_SurfaceType typS1 = myHS1->Surface().GetType();
 356         GeomAbs_SurfaceType typS2 = myHS2->Surface().GetType();
 357         if( typS1 == GeomAbs_SurfaceOfExtrusion ||
 358           typS1 == GeomAbs_OffsetSurface ||
 359           typS1 == GeomAbs_SurfaceOfRevolution ||
 360           typS2 == GeomAbs_SurfaceOfExtrusion ||
 361           typS2 == GeomAbs_OffsetSurface ||
 362           typS2 == GeomAbs_SurfaceOfRevolution) {
 363             sline.Append(newc);
 364             slineS1.Append(H1);
 365             slineS2.Append(H1);
 366             continue;
 367         }
 368         Standard_Boolean bFNIt, bLPIt;
 369         Standard_Real aTestPrm, dT=100.;
 370         Standard_Real u1, v1, u2, v2, TolX;
 371         //
 372         bFNIt=Precision::IsNegativeInfinite(fprm);
 373         bLPIt=Precision::IsPositiveInfinite(lprm);
 374
 375         aTestPrm=0.;
 376
 377         if (bFNIt && !bLPIt) {
 378           aTestPrm=lprm-dT;
 379         }
 380         else if (!bFNIt && bLPIt) {
 381           aTestPrm=fprm+dT;
 382         }
 383         //
 384         gp_Pnt ptref(newc->Value(aTestPrm));
 385         //
 386         TolX = Precision::Confusion();
 387         Parameters(myHS1, myHS2, ptref,  u1, v1, u2, v2);
 388         ok = (dom1->Classify(gp_Pnt2d(u1, v1), TolX) != TopAbs_OUT);
 389         if(ok) {
 390           ok = (dom2->Classify(gp_Pnt2d(u2,v2),TolX) != TopAbs_OUT);
 391         }
 392         if (ok) {
 393           sline.Append(newc);
 394           slineS1.Append(H1);
 395           slineS2.Append(H1);
 396         }
 397       }
 398     }// end of for (i=1; i<=myLConstruct.NbParts(); i++)
 399   }// case IntPatch_Lin:  case IntPatch_Parabola:  case IntPatch_Hyperbola:
 400   break;
 401
 402                            //########################################
 403                            // Circle and Ellipse
 404                            //########################################
 405   case IntPatch_Circle:
 406   case IntPatch_Ellipse: {
 407
 408     if (typl == IntPatch_Circle) {
 409       newc = new Geom_Circle
 410         (Handle(IntPatch_GLine)::DownCast(L)->Circle());
 411     }
 412     else {
 413       newc = new Geom_Ellipse
 414         (Handle(IntPatch_GLine)::DownCast(L)->Ellipse());
 415     }
 416     //
 417     Standard_Real aPeriod, aRealEpsilon;
 418     //
 419     aRealEpsilon=RealEpsilon();
 420     aPeriod=M_PI+M_PI;
 421     //
 422     aNbParts=myLConstruct.NbParts();
 423     //
 424     for (i=1; i<=aNbParts; i++) {
 425       myLConstruct.Part(i, fprm, lprm);
 426       //
 427       if (Abs(fprm) > aRealEpsilon || Abs(lprm-aPeriod) > aRealEpsilon) {
 428         //==============================================
 429         Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
 430         //
 431         sline.Append(aTC3D);
 432         //
 433         fprm=aTC3D->FirstParameter();
 434         lprm=aTC3D->LastParameter ();
 435         ////
 436         if(myApprox1) {
 437           Handle (Geom2d_Curve) C2d;
 438           BuildPCurves(fprm,lprm,Tolpc,myHS1->ChangeSurface().Surface(),newc,C2d);
 439           if(Tolpc>myTolReached2d || myTolReached2d==0.) {
 440             myTolReached2d=Tolpc;
 441           }
 442           slineS1.Append(C2d);
 443         }
 444         else { ////
 445           slineS1.Append(H1);
 446         }
 447         //
 448         if(myApprox2) {
 449           Handle (Geom2d_Curve) C2d;
 450           BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
 451           if(Tolpc>myTolReached2d || myTolReached2d==0) {
 452             myTolReached2d=Tolpc;
 453           }
 454           slineS2.Append(C2d);
 455         }
 456         else {
 457           slineS2.Append(H1);
 458         }
 459         //==============================================
 460       } //if (Abs(fprm) > RealEpsilon() || Abs(lprm-2.*M_PI) > RealEpsilon())
 461       //
 462       else {//  on regarde si on garde
 463         //
 464         if (aNbParts==1) {
 465           if (Abs(fprm) < RealEpsilon() &&  Abs(lprm-2.*M_PI) < RealEpsilon()) {
 466             Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
 467             //
 468             sline.Append(aTC3D);
 469             fprm=aTC3D->FirstParameter();
 470             lprm=aTC3D->LastParameter ();
 471
 472             if(myApprox1) {
 473               Handle (Geom2d_Curve) C2d;
 474               BuildPCurves(fprm,lprm,Tolpc,myHS1->ChangeSurface().Surface(),newc,C2d);
 475               if(Tolpc>myTolReached2d || myTolReached2d==0) {
 476                 myTolReached2d=Tolpc;
 477               }
 478               slineS1.Append(C2d);
 479             }
 480             else { ////
 481               slineS1.Append(H1);
 482             }
 483
 484             if(myApprox2) {
 485               Handle (Geom2d_Curve) C2d;
 486               BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
 487               if(Tolpc>myTolReached2d || myTolReached2d==0) {
 488                 myTolReached2d=Tolpc;
 489               }
 490               slineS2.Append(C2d);
 491             }
 492             else {
 493               slineS2.Append(H1);
 494             }
 495             break;
 496           }
 497         }
 498         //
 499         Standard_Real aTwoPIdiv17, u1, v1, u2, v2, TolX;
 500         //
 501         aTwoPIdiv17=2.*M_PI/17.;
 502         //
 503         for (j=0; j<=17; j++) {
 504           gp_Pnt ptref (newc->Value (j*aTwoPIdiv17));
 505           TolX = Precision::Confusion();
 506
 507           Parameters(myHS1, myHS2, ptref, u1, v1, u2, v2);
 508           ok = (dom1->Classify(gp_Pnt2d(u1,v1),TolX) != TopAbs_OUT);
 509           if(ok) {
 510             ok = (dom2->Classify(gp_Pnt2d(u2,v2),TolX) != TopAbs_OUT);
 511           }
 512           if (ok) {
 513             sline.Append(newc);
 514             //==============================================
 515             if(myApprox1) {
 516               Handle (Geom2d_Curve) C2d;
 517               BuildPCurves(fprm, lprm, Tolpc, myHS1->ChangeSurface().Surface(), newc, C2d);
 518               if(Tolpc>myTolReached2d || myTolReached2d==0) {
 519                 myTolReached2d=Tolpc;
 520               }
 521               slineS1.Append(C2d);
 522             }
 523             else {
 524               slineS1.Append(H1);
 525             }
 526
 527             if(myApprox2) {
 528               Handle (Geom2d_Curve) C2d;
 529               BuildPCurves(fprm, lprm, Tolpc,myHS2->ChangeSurface().Surface(), newc, C2d);
 530               if(Tolpc>myTolReached2d || myTolReached2d==0) {
 531                 myTolReached2d=Tolpc;
 532               }
 533               slineS2.Append(C2d);
 534             }
 535             else {
 536               slineS2.Append(H1);
 537             }
 538             break;
 539           }//  end of if (ok) {
 540         }//  end of for (Standard_Integer j=0; j<=17; j++)
 541       }//  end of else { on regarde si on garde
 542     }// for (i=1; i<=myLConstruct.NbParts(); i++)
 543   }// IntPatch_Circle: IntPatch_Ellipse
 544   break;
 545
 546                          //########################################
 547                          // Analytic
 548                          //########################################
 549   case IntPatch_Analytic: {
 550     IntSurf_Quadric quad1,quad2;
 551     //
 552     GetQuadric(myHS1, quad1);
 553     GetQuadric(myHS2, quad2);
 554     //=========
 555     IntPatch_ALineToWLine convert (quad1, quad2);
 556
 557     if (!myApprox) {
 558       Handle(Geom2d_BSplineCurve) aH1, aH2;
 559       //
 560       aNbParts=myLConstruct.NbParts();
 561       for (i=1; i<=aNbParts; i++) {
 562         myLConstruct.Part(i, fprm, lprm);
 563         Handle(IntPatch_WLine) WL =
 564           convert.MakeWLine(Handle(IntPatch_ALine)::DownCast(L), fprm, lprm);
 565         //
 566         if(myApprox1) {
 567           aH1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
 568         }
 569
 570         if(myApprox2) {
 571           aH2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
 572         }
 573         sline.Append(MakeBSpline(WL,1,WL->NbPnts()));
 574         slineS1.Append(aH1);
 575         slineS2.Append(aH2);
 576       }
 577     } // if (!myApprox)
 578
 579     else { // myApprox=TRUE
 580       GeomInt_WLApprox theapp3d;
 581       Standard_Real tol2d = myTolApprox;
 582       //
 583       theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
 584
 585       aNbParts=myLConstruct.NbParts();
 586       for (i=1; i<=aNbParts; i++) {
 587         myLConstruct.Part(i, fprm, lprm);
 588         Handle(IntPatch_WLine) WL =
 589           convert.MakeWLine(Handle(IntPatch_ALine):: DownCast(L),fprm,lprm);
 590
 591         theapp3d.Perform(myHS1,myHS2,WL,Standard_True,myApprox1,myApprox2, 1, WL->NbPnts());
 592         if (!theapp3d.IsDone()) {
 593           //
 594           Handle(Geom2d_BSplineCurve) aH1, aH2;
 595
 596           if(myApprox1) {
 597             aH1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
 598           }
 599
 600           if(myApprox2) {
 601             aH2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
 602           }
 603           sline.Append(MakeBSpline(WL,1,WL->NbPnts()));
 604           slineS1.Append(aH1);
 605           slineS2.Append(aH1);
 606         }
 607         //
 608         else {
 609           if(myApprox1 || myApprox2) {
 610             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0) {
 611               myTolReached2d = theapp3d.TolReached2d();
 612             }
 613           }
 614
 615           if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0) {
 616             myTolReached3d = theapp3d.TolReached3d();
 617           }
 618
 619           Standard_Integer aNbMultiCurves, nbpoles;
 620           aNbMultiCurves=theapp3d.NbMultiCurves();
 621           for (j=1; j<=aNbMultiCurves; j++) {
 622             const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
 623
 624             nbpoles = mbspc.NbPoles();
 625             TColgp_Array1OfPnt tpoles(1, nbpoles);
 626             mbspc.Curve(1, tpoles);
 627             Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
 628               mbspc.Knots(),
 629               mbspc.Multiplicities(),
 630               mbspc.Degree());
 631
 632             GeomLib_CheckBSplineCurve Check(BS, myTolCheck, myTolAngCheck);
 633             Check.FixTangent(Standard_True,Standard_True);
 634             //
 635             sline.Append(BS);
 636             //
 637             if(myApprox1) {
 638               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
 639               mbspc.Curve(2,tpoles2d);
 640               Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
 641                 mbspc.Knots(),
 642                 mbspc.Multiplicities(),
 643                 mbspc.Degree());
 644
 645               GeomLib_Check2dBSplineCurve newCheck(BS2,myTolCheck,myTolAngCheck);
 646               newCheck.FixTangent(Standard_True,Standard_True);
 647               slineS1.Append(BS2);
 648             }
 649             else {
 650               slineS1.Append(H1);
 651             }
 652
 653             if(myApprox2) {
 654               TColgp_Array1OfPnt2d tpoles2d(1, nbpoles);
 655               Standard_Integer TwoOrThree;
 656               TwoOrThree=myApprox1 ? 3 : 2;
 657               mbspc.Curve(TwoOrThree, tpoles2d);
 658               Handle(Geom2d_BSplineCurve) BS2 =new Geom2d_BSplineCurve(tpoles2d,
 659                 mbspc.Knots(),
 660                 mbspc.Multiplicities(),
 661                 mbspc.Degree());
 662
 663               GeomLib_Check2dBSplineCurve newCheck(BS2,myTolCheck,myTolAngCheck);
 664               newCheck.FixTangent(Standard_True,Standard_True);
 665               //
 666               slineS2.Append(BS2);
 667             }
 668             else {
 669               slineS2.Append(H1);
 670             }
 671             //
 672           }// for (j=1; j<=aNbMultiCurves; j++) {
 673         }// else from if (!theapp3d.IsDone())
 674       }// for (i=1; i<=aNbParts; i++) {
 675     }// else { // myApprox=TRUE
 676   }// case IntPatch_Analytic:
 677   break;
 678
 679                           //########################################
 680                           // Walking
 681                           //########################################
 682   case IntPatch_Walking:{
 683     Handle(IntPatch_WLine) WL =
 684       Handle(IntPatch_WLine)::DownCast(L);
 685     //
 686     Standard_Integer ifprm, ilprm;
 687     //
 688     if (!myApprox) {
 689       aNbParts=myLConstruct.NbParts();
 690       for (i=1; i<=aNbParts; i++) {
 691         myLConstruct.Part(i, fprm, lprm);
 692         ifprm=(Standard_Integer)fprm;
 693         ilprm=(Standard_Integer)lprm;
 694         //
 695         Handle(Geom2d_BSplineCurve) aH1, aH2;
 696
 697         if(myApprox1) {
 698           aH1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
 699         }
 700         if(myApprox2) {
 701           aH2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
 702         }
 703         //
 704         Handle(Geom_Curve) aBSp=MakeBSpline(WL, ifprm, ilprm);
 705         //
 706         sline.Append(aBSp);
 707         slineS1.Append(aH1);
 708         slineS2.Append(aH2);
 709       }
 710     }
 711     //
 712     else {
 713       Standard_Boolean bIsDecomposited;
 714       Standard_Integer nbiter, aNbSeqOfL;
 715       GeomInt_WLApprox theapp3d;
 716       IntPatch_SequenceOfLine aSeqOfL;
 717       Standard_Real tol2d, aTolSS;
 718       //
 719       tol2d = myTolApprox;
 720       aTolSS=2.e-7;
 721       if(myHS1 == myHS2) {
 722         theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False);
 723       }
 724       else {
 725         theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
 726       }
 727       //
 728       bIsDecomposited =
 729         DecompositionOfWLine(WL, myHS1, myHS2, aTolSS, myLConstruct, aSeqOfL);
 730       //
 731       aNbParts=myLConstruct.NbParts();
 732       aNbSeqOfL=aSeqOfL.Length();
 733       //
 734       nbiter = (bIsDecomposited) ? aNbSeqOfL : aNbParts;
 735       //
 736       for(i = 1; i <= nbiter; i++) {
 737         if(bIsDecomposited) {
 738           WL = Handle(IntPatch_WLine)::DownCast(aSeqOfL.Value(i));
 739           ifprm = 1;
 740           ilprm = WL->NbPnts();
 741         }
 742         else {
 743           myLConstruct.Part(i, fprm, lprm);
 744           ifprm = (Standard_Integer)fprm;
 745           ilprm = (Standard_Integer)lprm;
 746         }
 747         //-- lbr :
 748         //-- Si une des surfaces est un plan , on approxime en 2d
 749         //-- sur cette surface et on remonte les points 2d en 3d.
 750         GeomAbs_SurfaceType typs1, typs2;
 751         typs1 = myHS1->Surface().GetType();
 752         typs2 = myHS2->Surface().GetType();
 753         //
 754         if(typs1 == GeomAbs_Plane) {
 755           theapp3d.Perform(myHS1, myHS2, WL, Standard_False,
 756             Standard_True, myApprox2,
 757             ifprm,  ilprm);
 758         }
 759         else if(typs2 == GeomAbs_Plane) {
 760           theapp3d.Perform(myHS1,myHS2,WL,Standard_False,
 761             myApprox1,Standard_True,
 762             ifprm,  ilprm);
 763         }
 764         else {
 765           //
 766           if (myHS1 != myHS2){
 767             if ((typs1==GeomAbs_BezierSurface || typs1==GeomAbs_BSplineSurface) &&
 768               (typs2==GeomAbs_BezierSurface || typs2==GeomAbs_BSplineSurface)) {
 769
 770                 theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
 771                 //Standard_Boolean bUseSurfaces;
 772                 //bUseSurfaces=NotUseSurfacesForApprox(myFace1, myFace2, WL, ifprm,  ilprm);
 773                 //if (bUseSurfaces) {
 774                 //theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False);
 775                 //}
 776             }
 777           }
 778           //
 779           theapp3d.Perform(myHS1,myHS2,WL,Standard_True,
 780             myApprox1,myApprox2,
 781             ifprm,  ilprm);
 782         }
 783
 784         if (!theapp3d.IsDone()) {
 785           //
 786           Handle(Geom2d_BSplineCurve) aH1, aH2;
 787           //
 788           Handle(Geom_Curve) aBSp=MakeBSpline(WL, ifprm, ilprm);
 789           if(myApprox1) {
 790             aH1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
 791           }
 792           if(myApprox2) {
 793             aH2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
 794           }
 795           //
 796           sline.Append(aBSp);
 797           slineS1.Append(aH1);
 798           slineS2.Append(aH2);
 799         }//if (!theapp3d.IsDone())
 800
 801         else {
 802           if(myApprox1 || myApprox2 || (typs1==GeomAbs_Plane || typs2==GeomAbs_Plane)) {
 803             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0.) {
 804               myTolReached2d = theapp3d.TolReached2d();
 805             }
 806           }
 807           if(typs1==GeomAbs_Plane || typs2==GeomAbs_Plane) {
 808             myTolReached3d = myTolReached2d;
 809           }
 810           else  if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0.) {
 811             myTolReached3d = theapp3d.TolReached3d();
 812           }
 813
 814           Standard_Integer aNbMultiCurves, nbpoles;
 815           //
 816           aNbMultiCurves=theapp3d.NbMultiCurves();
 817           for (j=1; j<=aNbMultiCurves; j++) {
 818             if(typs1 == GeomAbs_Plane) {
 819               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
 820               nbpoles = mbspc.NbPoles();
 821
 822               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
 823               TColgp_Array1OfPnt   tpoles(1,nbpoles);
 824
 825               mbspc.Curve(1,tpoles2d);
 826               const gp_Pln&  Pln = myHS1->Surface().Plane();
 827               //
 828               Standard_Integer ik;
 829               for(ik = 1; ik<= nbpoles; ik++) {
 830                 tpoles.SetValue(ik,
 831                   ElSLib::Value(tpoles2d.Value(ik).X(),
 832                   tpoles2d.Value(ik).Y(),
 833                   Pln));
 834               }
 835               //
 836               Handle(Geom_BSplineCurve) BS =
 837                 new Geom_BSplineCurve(tpoles,
 838                 mbspc.Knots(),
 839                 mbspc.Multiplicities(),
 840                 mbspc.Degree());
 841               GeomLib_CheckBSplineCurve Check(BS,myTolCheck,myTolAngCheck);
 842               Check.FixTangent(Standard_True, Standard_True);
 843               //
 844               sline.Append(BS);
 845               //
 846               if(myApprox1) {
 847                 Handle(Geom2d_BSplineCurve) BS1 =
 848                   new Geom2d_BSplineCurve(tpoles2d,
 849                   mbspc.Knots(),
 850                   mbspc.Multiplicities(),
 851                   mbspc.Degree());
 852                 GeomLib_Check2dBSplineCurve Check1(BS1,myTolCheck,myTolAngCheck);
 853                 Check1.FixTangent(Standard_True,Standard_True);
 854                 //
 855                 AdjustUPeriodic (aS1, BS1);
 856                 //
 857                 slineS1.Append(BS1);
 858               }
 859               else {
 860                 slineS1.Append(H1);
 861               }
 862
 863               if(myApprox2) {
 864                 mbspc.Curve(2, tpoles2d);
 865
 866                 Handle(Geom2d_BSplineCurve) BS2 = new Geom2d_BSplineCurve(tpoles2d,
 867                   mbspc.Knots(),
 868                   mbspc.Multiplicities(),
 869                   mbspc.Degree());
 870                 GeomLib_Check2dBSplineCurve newCheck(BS2,myTolCheck,myTolAngCheck);
 871                 newCheck.FixTangent(Standard_True,Standard_True);
 872                 //
 873                 AdjustUPeriodic (aS2, BS2);
 874                 //
 875                 slineS2.Append(BS2);
 876               }
 877               else {
 878                 slineS2.Append(H1);
 879               }
 880             }//if(typs1 == GeomAbs_Plane)
 881             //
 882             else if(typs2 == GeomAbs_Plane) {
 883               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
 884               nbpoles = mbspc.NbPoles();
 885
 886               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
 887               TColgp_Array1OfPnt   tpoles(1,nbpoles);
 888               mbspc.Curve((myApprox1==Standard_True)? 2 : 1,tpoles2d);
 889               const gp_Pln&  Pln = myHS2->Surface().Plane();
 890               //
 891               Standard_Integer ik;
 892               for(ik = 1; ik<= nbpoles; ik++) {
 893                 tpoles.SetValue(ik,
 894                   ElSLib::Value(tpoles2d.Value(ik).X(),
 895                   tpoles2d.Value(ik).Y(),
 896                   Pln));
 897
 898               }
 899               //
 900               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
 901                 mbspc.Knots(),
 902                 mbspc.Multiplicities(),
 903                 mbspc.Degree());
 904               GeomLib_CheckBSplineCurve Check(BS,myTolCheck,myTolAngCheck);
 905               Check.FixTangent(Standard_True,Standard_True);
 906               //
 907               sline.Append(BS);
 908               //
 909               if(myApprox2) {
 910                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
 911                   mbspc.Knots(),
 912                   mbspc.Multiplicities(),
 913                   mbspc.Degree());
 914                 GeomLib_Check2dBSplineCurve Check1(BS1,myTolCheck,myTolAngCheck);
 915                 Check1.FixTangent(Standard_True,Standard_True);
 916                 //
 917                 //
 918                 AdjustUPeriodic (aS2, BS1);
 919                 //
 920                 slineS2.Append(BS1);
 921               }
 922               else {
 923                 slineS2.Append(H1);
 924               }
 925
 926               if(myApprox1) {
 927                 mbspc.Curve(1,tpoles2d);
 928                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
 929                   mbspc.Knots(),
 930                   mbspc.Multiplicities(),
 931                   mbspc.Degree());
 932                 GeomLib_Check2dBSplineCurve Check2(BS2,myTolCheck,myTolAngCheck);
 933                 Check2.FixTangent(Standard_True,Standard_True);
 934                 //
 935                 //
 936                 AdjustUPeriodic (aS1, BS2);
 937                 //
 938                 slineS1.Append(BS2);
 939               }
 940               else {
 941                 slineS1.Append(H1);
 942               }
 943             } // else if(typs2 == GeomAbs_Plane)
 944             //
 945             else { // typs1!=GeomAbs_Plane && typs2!=GeomAbs_Plane
 946               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
 947               nbpoles = mbspc.NbPoles();
 948               TColgp_Array1OfPnt tpoles(1,nbpoles);
 949               mbspc.Curve(1,tpoles);
 950               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
 951                 mbspc.Knots(),
 952                 mbspc.Multiplicities(),
 953                 mbspc.Degree());
 954               GeomLib_CheckBSplineCurve Check(BS,myTolCheck,myTolAngCheck);
 955               Check.FixTangent(Standard_True,Standard_True);
 956               //
 957               //Check IsClosed()
 958               Standard_Real aDist = Max(BS->StartPoint().XYZ().SquareModulus(),
 959                 BS->EndPoint().XYZ().SquareModulus());
 960               Standard_Real eps = Epsilon(aDist);
 961               if(BS->StartPoint().SquareDistance(BS->EndPoint()) < 2.*eps &&
 962                 !BS->IsClosed() && !BS->IsPeriodic())
 963               {
 964                 //force Closed()
 965                 gp_Pnt aPm((BS->Pole(1).XYZ() + BS->Pole(BS->NbPoles()).XYZ()) / 2.);
 966                 BS->SetPole(1, aPm);
 967                 BS->SetPole(BS->NbPoles(), aPm);
 968               }
 969               sline.Append(BS);
 970
 971               if(myApprox1) {
 972                 TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
 973                 mbspc.Curve(2,tpoles2d);
 974                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
 975                   mbspc.Knots(),
 976                   mbspc.Multiplicities(),
 977                   mbspc.Degree());
 978                 GeomLib_Check2dBSplineCurve newCheck(BS1,myTolCheck,myTolAngCheck);
 979                 newCheck.FixTangent(Standard_True,Standard_True);
 980                 //
 981                 AdjustUPeriodic (aS1, BS1);
 982                 //
 983                 slineS1.Append(BS1);
 984               }
 985               else {
 986                 slineS1.Append(H1);
 987               }
 988               if(myApprox2) {
 989                 TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
 990                 mbspc.Curve((myApprox1==Standard_True)? 3 : 2,tpoles2d);
 991                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
 992                   mbspc.Knots(),
 993                   mbspc.Multiplicities(),
 994                   mbspc.Degree());
 995                 GeomLib_Check2dBSplineCurve newCheck(BS2,myTolCheck,myTolAngCheck);
 996                 newCheck.FixTangent(Standard_True,Standard_True);
 997                 //
 998                 AdjustUPeriodic (aS2, BS2);
 999                 //
1000                 slineS2.Append(BS2);
1001               }
1002               else {
1003                 slineS2.Append(H1);
1004               }
1005             }// else { // typs1!=GeomAbs_Plane && typs2!=GeomAbs_Plane
1006           }// for (j=1; j<=aNbMultiCurves; j++
1007         }
1008       }
1009     }
1010   }
1011   break;
1012
1013   case IntPatch_Restriction:
1014     {
1015       GeomAbs_SurfaceType typS1 = myHS1->Surface().GetType();
1016       GeomAbs_SurfaceType typS2 = myHS2->Surface().GetType();
1017       Standard_Boolean isAnalS1 = Standard_False;
1018       switch (typS1)
1019       {
1020       case GeomAbs_Plane:
1021       case GeomAbs_Cylinder:
1022       case GeomAbs_Sphere:
1023       case GeomAbs_Cone:
1024       case GeomAbs_Torus: isAnalS1 = Standard_True; break;
1025       default: break;
1026       }
1027
1028       Standard_Integer isAnalS2 = Standard_False;
1029       switch (typS2)
1030       {
1031       case GeomAbs_Plane:
1032       case GeomAbs_Cylinder:
1033       case GeomAbs_Sphere:
1034       case GeomAbs_Cone:
1035       case GeomAbs_Torus: isAnalS2 = Standard_True; break;
1036       default: break;
1037       }
1038
1039       Handle(IntPatch_RLine) RL =
1040         Handle(IntPatch_RLine)::DownCast(L);
1041       Handle(Geom_Curve) aC3d;
1042       Handle(Geom2d_Curve) aC2d1, aC2d2;
1043       Standard_Real aTolReached;
1044       TreatRLine(RL, myHS1, myHS2, aC3d,
1045                   aC2d1, aC2d2, aTolReached);
1046
1047       if(aC3d.IsNull())
1048         break;
1049
1050       Bnd_Box2d aBox1, aBox2;
1051
1052       const Standard_Real aU1f = myHS1->FirstUParameter(),
1053                           aV1f = myHS1->FirstVParameter(),
1054                           aU1l = myHS1->LastUParameter(),
1055                           aV1l = myHS1->LastVParameter();
1056       const Standard_Real aU2f = myHS2->FirstUParameter(),
1057                           aV2f = myHS2->FirstVParameter(),
1058                           aU2l = myHS2->LastUParameter(),
1059                           aV2l = myHS2->LastVParameter();
1060
1061       aBox1.Add(gp_Pnt2d(aU1f, aV1f));
1062       aBox1.Add(gp_Pnt2d(aU1l, aV1l));
1063       aBox2.Add(gp_Pnt2d(aU2f, aV2f));
1064       aBox2.Add(gp_Pnt2d(aU2l, aV2l));
1065
1066       GeomInt_VectorOfReal anArrayOfParameters;
1067
1068       //We consider here that the intersection line is same-parameter-line
1069       anArrayOfParameters.Append(aC3d->FirstParameter());
1070       anArrayOfParameters.Append(aC3d->LastParameter());
1071
1072       TrimILineOnSurfBoundaries(aC2d1, aC2d2, aBox1, aBox2, anArrayOfParameters);
1073
1074       const Standard_Integer aNbIntersSolutionsm1 = anArrayOfParameters.Length() - 1;
1075
1076       //Trim RLine found.
1077       for(Standard_Integer anInd = 0; anInd < aNbIntersSolutionsm1; anInd++)
1078       {
1079         const Standard_Real aParF = anArrayOfParameters(anInd),
1080                             aParL = anArrayOfParameters(anInd+1);
1081
1082         if((aParL - aParF) <= Precision::PConfusion())
1083           continue;
1084
1085         const Standard_Real aPar = 0.5*(aParF + aParL);
1086         gp_Pnt2d aPt;
1087
1088         Handle(Geom2d_Curve) aCurv2d1, aCurv2d2;
1089         if(!aC2d1.IsNull())
1090         {
1091           aC2d1->D0(aPar, aPt);
1092
1093           if(aBox1.IsOut(aPt))
1094             continue;
1095
1096           if(myApprox1)
1097             aCurv2d1 = new Geom2d_TrimmedCurve(aC2d1, aParF, aParL);
1098         }
1099
1100         if(!aC2d2.IsNull())
1101         {
1102           aC2d2->D0(aPar, aPt);
1103
1104           if(aBox2.IsOut(aPt))
1105             continue;
1106
1107           if(myApprox2)
1108             aCurv2d2 = new Geom2d_TrimmedCurve(aC2d2, aParF, aParL);
1109         }
1110
1111         Handle(Geom_Curve) aCurv3d = new Geom_TrimmedCurve(aC3d, aParF, aParL);
1112
1113         sline.Append(aCurv3d);
1114         slineS1.Append(aCurv2d1);
1115         slineS2.Append(aCurv2d2);
1116       }
1117     }
1118     break;
1119   }
1120 }
1121 //
1122 //=======================================================================
1123 //function : AdjustUPeriodic
1124 //purpose  :
1125 //=======================================================================
1126  void AdjustUPeriodic (const Handle(Geom_Surface)& aS, Handle(Geom2d_Curve)& aC2D)
1127 {
1128   if (aC2D.IsNull() || !aS->IsUPeriodic())
1129     return;
1130   //
1131   const Standard_Real aEps=Precision::PConfusion();//1.e-9
1132   const Standard_Real aEpsilon=Epsilon(10.);//1.77e-15
1133   //
1134   Standard_Real umin,umax,vmin,vmax;
1135   aS->Bounds(umin,umax,vmin,vmax);
1136   const Standard_Real aPeriod = aS->UPeriod();
1137
1138   const Standard_Real aT1=aC2D->FirstParameter();
1139   const Standard_Real aT2=aC2D->LastParameter();
1140   const Standard_Real aTx=aT1+0.467*(aT2-aT1);
1141   const gp_Pnt2d aPx=aC2D->Value(aTx);
1142   //
1143   Standard_Real aUx=aPx.X();
1144   if (fabs(aUx)<aEpsilon)
1145     aUx=0.;
1146   if (fabs(aUx-aPeriod)<aEpsilon)
1147     aUx=aPeriod;
1148   //
1149   Standard_Real dU=0.;
1150   while(aUx <(umin-aEps)) {
1151     aUx+=aPeriod;
1152     dU+=aPeriod;
1153   }
1154   while(aUx>(umax+aEps)) {
1155     aUx-=aPeriod;
1156     dU-=aPeriod;
1157   }
1158   //
1159   if (dU!=0.) {
1160     gp_Vec2d aV2D(dU, 0.);
1161     aC2D->Translate(aV2D);
1162   }
1163 }
1164 //
1165 //
1166 //=======================================================================
1167 //function : GetQuadric
1168 //purpose  :
1169 //=======================================================================
1170  void GetQuadric(const Handle(GeomAdaptor_HSurface)& HS1, IntSurf_Quadric& quad1)
1171 {
1172   switch (HS1->Surface().GetType())
1173   {
1174     case GeomAbs_Plane:    quad1.SetValue(HS1->Surface().Plane()); break;
1175     case GeomAbs_Cylinder: quad1.SetValue(HS1->Surface().Cylinder()); break;
1176     case GeomAbs_Cone:     quad1.SetValue(HS1->Surface().Cone()); break;
1177     case GeomAbs_Sphere:   quad1.SetValue(HS1->Surface().Sphere()); break;
1178     case GeomAbs_Torus:    quad1.SetValue(HS1->Surface().Torus()); break;
1179     default: Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
1180   }
1181 }
1182
1183 //=======================================================================
1184 //function : Parameters
1185 //purpose  :
1186 //=======================================================================
1187  void Parameters(const Handle(GeomAdaptor_HSurface)& HS1,
1188                  const Handle(GeomAdaptor_HSurface)& HS2,
1189                  const gp_Pnt& Ptref,
1190                  Standard_Real& U1,
1191                  Standard_Real& V1,
1192                  Standard_Real& U2,
1193                  Standard_Real& V2)
1194 {
1195   IntSurf_Quadric quad1,quad2;
1196   //
1197   GetQuadric(HS1, quad1);
1198   GetQuadric(HS2, quad2);
1199   //
1200   quad1.Parameters(Ptref,U1,V1);
1201   quad2.Parameters(Ptref,U2,V2);
1202 }
1203
1204 //=======================================================================
1205 //function : MakeBSpline
1206 //purpose  :
1207 //=======================================================================
1208 Handle(Geom_Curve) MakeBSpline  (const Handle(IntPatch_WLine)& WL,
1209                                  const Standard_Integer ideb,
1210                                  const Standard_Integer ifin)
1211 {
1212   const Standard_Integer nbpnt = ifin-ideb+1;
1213   TColgp_Array1OfPnt poles(1,nbpnt);
1214   TColStd_Array1OfReal knots(1,nbpnt);
1215   TColStd_Array1OfInteger mults(1,nbpnt);
1216   Standard_Integer i = 1, ipidebm1 = ideb;
1217   for(; i<=nbpnt; ipidebm1++, i++)
1218   {
1219     poles(i) = WL->Point(ipidebm1).Value();
1220     mults(i) = 1;
1221     knots(i) = i-1;
1222   }
1223   mults(1) = mults(nbpnt) = 2;
1224   return new Geom_BSplineCurve(poles,knots,mults,1);
1225 }
1226
1227 //=======================================================================
1228 //function : MakeBSpline2d
1229 //purpose  :
1230 //=======================================================================
1231 Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
1232                                           const Standard_Integer ideb,
1233                                           const Standard_Integer ifin,
1234                                           const Standard_Boolean onFirst)
1235 {
1236   const Standard_Integer nbpnt = ifin-ideb+1;
1237   TColgp_Array1OfPnt2d poles(1,nbpnt);
1238   TColStd_Array1OfReal knots(1,nbpnt);
1239   TColStd_Array1OfInteger mults(1,nbpnt);
1240   Standard_Integer i = 1, ipidebm1 = ideb;
1241   for(; i <= nbpnt; ipidebm1++, i++)
1242   {
1243     Standard_Real U, V;
1244     if(onFirst)
1245           theWLine->Point(ipidebm1).ParametersOnS1(U, V);
1246     else
1247           theWLine->Point(ipidebm1).ParametersOnS2(U, V);
1248     poles(i).SetCoord(U, V);
1249     mults(i) = 1;
1250     knots(i) = i-1;
1251   }
1252   mults(1) = mults(nbpnt) = 2;
1253   return new Geom2d_BSplineCurve(poles,knots,mults,1);
1254 }
1255
1256 //=========================================================================
1257 // static function : ComputePurgedWLine
1258 // purpose : Removes equal points (leave one of equal points) from theWLine
1259 //           and recompute vertex parameters.
1260 //           Returns new WLine or null WLine if the number
1261 //           of the points is less than 2.
1262 //=========================================================================
1263 Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine)
1264 {
1265   Handle(IntPatch_WLine) aResult;
1266   Handle(IntPatch_WLine) aLocalWLine;
1267   Handle(IntPatch_WLine) aTmpWLine = theWLine;
1268
1269   Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
1270   aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
1271   Standard_Integer i, k, v, nb, nbvtx;
1272   nbvtx = theWLine->NbVertex();
1273   nb = theWLine->NbPnts();
1274
1275   for(i = 1; i <= nb; i++) {
1276     aLineOn2S->Add(theWLine->Point(i));
1277   }
1278
1279   for(v = 1; v <= nbvtx; v++) {
1280     aLocalWLine->AddVertex(theWLine->Vertex(v));
1281   }
1282
1283   for(i = 1; i <= aLineOn2S->NbPoints(); i++) {
1284     Standard_Integer aStartIndex = i + 1;
1285     Standard_Integer anEndIndex = i + 5;
1286     nb = aLineOn2S->NbPoints();
1287     anEndIndex = (anEndIndex > nb) ? nb : anEndIndex;
1288
1289     if((aStartIndex >= nb) || (anEndIndex <= 1)) {
1290       continue;
1291     }
1292     k = aStartIndex;
1293
1294     while(k <= anEndIndex) {
1295
1296       if(i != k) {
1297         IntSurf_PntOn2S p1 = aLineOn2S->Value(i);
1298         IntSurf_PntOn2S p2 = aLineOn2S->Value(k);
1299
1300         if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) {
1301           aTmpWLine = aLocalWLine;
1302           aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
1303
1304           for(v = 1; v <= aTmpWLine->NbVertex(); v++) {
1305             IntPatch_Point aVertex = aTmpWLine->Vertex(v);
1306             Standard_Integer avertexindex = (Standard_Integer)aVertex.ParameterOnLine();
1307
1308             if(avertexindex >= k) {
1309               aVertex.SetParameter(aVertex.ParameterOnLine() - 1.);
1310             }
1311             aLocalWLine->AddVertex(aVertex);
1312           }
1313           aLineOn2S->RemovePoint(k);
1314           anEndIndex--;
1315           continue;
1316         }
1317       }
1318       k++;
1319     }
1320   }
1321
1322   if(aLineOn2S->NbPoints() > 1) {
1323     aResult = aLocalWLine;
1324   }
1325   return aResult;
1326 }
1327 //=======================================================================
1328 //function : DecompositionOfWLine
1329 //purpose  :
1330 //=======================================================================
1331 Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
1332                                       const Handle(GeomAdaptor_HSurface)&            theSurface1,
1333                                       const Handle(GeomAdaptor_HSurface)&            theSurface2,
1334                                       const Standard_Real aTolSum,
1335                                       const GeomInt_LineConstructor&                 theLConstructor,
1336                                       IntPatch_SequenceOfLine&                       theNewLines)
1337 {
1338   typedef NCollection_List<Standard_Integer> ListOfInteger;
1339   //have to use std::vector, not NCollection_Vector in order to use copy constructor of
1340   //ListOfInteger which will be created with specific allocator instance
1341   typedef std::vector<ListOfInteger, NCollection_StdAllocator<
1342       ListOfInteger> > ArrayOfListOfInteger;
1343
1344   Standard_Boolean bIsPrevPointOnBoundary, bIsCurrentPointOnBoundary;
1345   Standard_Integer nblines, aNbPnts, aNbParts, pit, i, j, aNbListOfPointIndex;
1346   Standard_Real aTol, umin, umax, vmin, vmax;
1347
1348   //an inc allocator, it will contain wasted space (upon list's Clear()) but it should
1349   //still be faster than the standard allocator, and wasted memory should not be
1350   //significant and will be limited by time span of this function;
1351   //this is a separate allocator from the anIncAlloc below what provides better data
1352   //locality in the latter (by avoiding wastes which will only be in anIdxAlloc)
1353   Handle(NCollection_IncAllocator) anIdxAlloc = new NCollection_IncAllocator();
1354   ListOfInteger aListOfPointIndex (anIdxAlloc);
1355
1356   //GeomAPI_ProjectPointOnSurf aPrj1, aPrj2;
1357   ProjectPointOnSurf aPrj1, aPrj2;
1358   Handle(Geom_Surface) aSurf1,  aSurf2;
1359   //
1360   aNbParts=theLConstructor.NbParts();
1361   aNbPnts=theWLine->NbPnts();
1362   //
1363   if((!aNbPnts) || (!aNbParts)){
1364     return Standard_False;
1365   }
1366   //
1367   Handle(NCollection_IncAllocator) anIncAlloc = new NCollection_IncAllocator();
1368   NCollection_StdAllocator<ListOfInteger> anAlloc (anIncAlloc);
1369   const ListOfInteger aDummy (anIncAlloc); //empty list to be copy constructed from
1370   ArrayOfListOfInteger anArrayOfLines (aNbPnts + 1, aDummy, anAlloc);
1371
1372   NCollection_LocalArray<Standard_Integer> anArrayOfLineTypeArr (aNbPnts + 1);
1373   Standard_Integer* anArrayOfLineType = anArrayOfLineTypeArr;
1374   //
1375   nblines = 0;
1376   aTol = Precision::Confusion();
1377   //
1378   aSurf1 = theSurface1->ChangeSurface().Surface();
1379   aSurf1->Bounds(umin, umax, vmin, vmax);
1380   aPrj1.Init(aSurf1, umin, umax, vmin, vmax);
1381   //
1382   aSurf2 = theSurface2->ChangeSurface().Surface();
1383   aSurf2->Bounds(umin, umax, vmin, vmax);
1384   aPrj2.Init(aSurf2, umin, umax, vmin, vmax);
1385   //
1386   //
1387   bIsPrevPointOnBoundary=Standard_False;
1388   for(pit=1; pit<=aNbPnts; pit++) {
1389     const IntSurf_PntOn2S& aPoint = theWLine->Point(pit);
1390     bIsCurrentPointOnBoundary=Standard_False;
1391     //
1392     // whether aPoint is on boundary or not
1393     //
1394     for(i=0; i<2; i++) {// exploration Surface 1,2
1395       Handle(GeomAdaptor_HSurface) aGASurface = (!i) ? theSurface1 : theSurface2;
1396       aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
1397       //
1398       for(j=0; j<2; j++) {// exploration of coordinate U,V
1399         Standard_Boolean isperiodic;
1400         //
1401         isperiodic = (!j) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
1402         if(!isperiodic) {
1403           continue;
1404         }
1405         //
1406         Standard_Real aResolution, aPeriod, alowerboundary, aupperboundary, U, V;
1407         Standard_Real aParameter, anoffset, anAdjustPar;
1408         Standard_Boolean bIsOnFirstBoundary, bIsPointOnBoundary;
1409         //
1410         aResolution = (!j) ? aGASurface->UResolution(aTol) : aGASurface->VResolution(aTol);
1411         aPeriod     = (!j) ? aGASurface->UPeriod() : aGASurface->VPeriod();
1412         alowerboundary = (!j) ? umin : vmin;
1413         aupperboundary = (!j) ? umax : vmax;
1414         U=0.;V=0.;//?
1415         //
1416         if(!i){
1417           aPoint.ParametersOnS1(U, V);
1418         }
1419         else{
1420           aPoint.ParametersOnS2(U, V);
1421         }
1422         //
1423         aParameter = (!j) ? U : V;
1424         anoffset=0.;
1425         anAdjustPar=AdjustPeriodic(aParameter, alowerboundary, aupperboundary, aPeriod, anoffset);
1426         //
1427         bIsOnFirstBoundary=Standard_True;
1428         //
1429         bIsPointOnBoundary=
1430           IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
1431
1432         if(bIsPointOnBoundary) {
1433           bIsCurrentPointOnBoundary = Standard_True;
1434           break;
1435         }
1436       }// for(j=0; j<2; j++)
1437
1438       if(bIsCurrentPointOnBoundary){
1439         break;
1440       }
1441     }// for(i=0; i<2; i++)
1442     //
1443     if((bIsCurrentPointOnBoundary != bIsPrevPointOnBoundary)) {
1444
1445       if(!aListOfPointIndex.IsEmpty()) {
1446         nblines++;
1447         anArrayOfLines[nblines] = aListOfPointIndex;
1448         anArrayOfLineType[nblines] = bIsPrevPointOnBoundary;
1449         aListOfPointIndex.Clear();
1450       }
1451       bIsPrevPointOnBoundary = bIsCurrentPointOnBoundary;
1452     }
1453     aListOfPointIndex.Append(pit);
1454   } // for(pit=1; pit<=aNbPnts; pit++)
1455   //
1456   aNbListOfPointIndex=aListOfPointIndex.Extent();
1457   if(aNbListOfPointIndex) {
1458     nblines++;
1459     anArrayOfLines[nblines].Assign (aListOfPointIndex);
1460     anArrayOfLineType[nblines] = bIsPrevPointOnBoundary;
1461     aListOfPointIndex.Clear();
1462   }
1463   //
1464   if(nblines <= 1){
1465     return Standard_False;
1466   }
1467   //
1468   // Correct wlines.begin
1469   Standard_Integer aLineType;
1470   TColStd_Array1OfListOfInteger anArrayOfLineEnds(1, nblines);
1471   Handle(IntSurf_LineOn2S) aSeqOfPntOn2S = new IntSurf_LineOn2S (new NCollection_IncAllocator());
1472   //
1473   for(i = 1; i <= nblines; i++) {
1474     aLineType=anArrayOfLineType[i];
1475     if(aLineType) {
1476       continue;
1477     }
1478     //
1479     const ListOfInteger& aListOfIndex = anArrayOfLines[i];
1480     if(aListOfIndex.Extent() < 2) {
1481       continue;
1482     }
1483     //
1484     TColStd_ListOfInteger aListOfFLIndex;
1485     Standard_Integer aneighbourindex, aLineTypeNeib;
1486     //
1487     for(j = 0; j < 2; j++) {// neighbour line choice
1488       aneighbourindex = (!j) ? (i-1) : (i+1);
1489       if((aneighbourindex < 1) || (aneighbourindex > nblines)){
1490         continue;
1491       }
1492       //
1493       aLineTypeNeib=anArrayOfLineType[aneighbourindex];
1494       if(!aLineTypeNeib){
1495         continue;
1496       }
1497       //
1498       const ListOfInteger& aNeighbour = anArrayOfLines[aneighbourindex];
1499       Standard_Integer anIndex = (!j) ? aNeighbour.Last() : aNeighbour.First();
1500       const IntSurf_PntOn2S& aPoint = theWLine->Point(anIndex);
1501       // check if need use derivative.begin .end [absence]
1502       //
1503       IntSurf_PntOn2S aNewP = aPoint;
1504       Standard_Integer surfit, parit;
1505       //
1506       for(surfit = 0; surfit < 2; ++surfit) {
1507
1508         Handle(GeomAdaptor_HSurface) aGASurface = (!surfit) ? theSurface1 : theSurface2;
1509
1510         umin = aGASurface->FirstUParameter();
1511         umax = aGASurface->LastUParameter();
1512         vmin = aGASurface->FirstVParameter();
1513         vmax = aGASurface->LastVParameter();
1514         Standard_Real U=0., V=0.;
1515
1516         if(!surfit) {
1517           aNewP.ParametersOnS1(U, V);
1518         }
1519         else {
1520           aNewP.ParametersOnS2(U, V);
1521         }
1522         //
1523         Standard_Integer nbboundaries = 0;
1524         Standard_Integer bIsUBoundary = Standard_False; // use if nbboundaries == 1
1525         Standard_Integer bIsFirstBoundary = Standard_False; // use if nbboundaries == 1
1526         //
1527         for(parit = 0; parit < 2; parit++) {
1528           Standard_Boolean isperiodic = (!parit) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
1529
1530           Standard_Real aResolution = (!parit) ? aGASurface->UResolution(aTol) : aGASurface->VResolution(aTol);
1531           Standard_Real alowerboundary = (!parit) ? umin : vmin;
1532           Standard_Real aupperboundary = (!parit) ? umax : vmax;
1533
1534           Standard_Real aParameter = (!parit) ? U : V;
1535           Standard_Boolean bIsOnFirstBoundary = Standard_True;
1536
1537           if(!isperiodic) {
1538             if(IsPointOnBoundary(aParameter, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary)) {
1539               bIsUBoundary = (!parit);
1540               bIsFirstBoundary = bIsOnFirstBoundary;
1541               nbboundaries++;
1542             }
1543           }
1544           else {
1545             Standard_Real aPeriod     = (!parit) ? aGASurface->UPeriod() : aGASurface->VPeriod();
1546             Standard_Real anoffset = 0.;
1547             Standard_Real anAdjustPar = AdjustPeriodic(aParameter, alowerboundary, aupperboundary, aPeriod, anoffset);
1548
1549             if(IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary)) {
1550               bIsUBoundary = (parit == 0);
1551               bIsFirstBoundary = bIsOnFirstBoundary;
1552               nbboundaries++;
1553             }
1554           }
1555         }
1556         //
1557         Standard_Boolean bComputeLineEnd = Standard_False;
1558
1559         if(nbboundaries == 2) {
1560           bComputeLineEnd = Standard_True;
1561         }
1562         else if(nbboundaries == 1) {
1563           Standard_Boolean isperiodic = (bIsUBoundary) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
1564
1565           if(isperiodic) {
1566             Standard_Real alowerboundary = (bIsUBoundary) ? umin : vmin;
1567             Standard_Real aupperboundary = (bIsUBoundary) ? umax : vmax;
1568             Standard_Real aPeriod     = (bIsUBoundary) ? aGASurface->UPeriod() : aGASurface->VPeriod();
1569             Standard_Real aParameter = (bIsUBoundary) ? U : V;
1570             Standard_Real anoffset = 0.;
1571             Standard_Real anAdjustPar = AdjustPeriodic(aParameter, alowerboundary, aupperboundary, aPeriod, anoffset);
1572
1573             Standard_Real adist = (bIsFirstBoundary) ? fabs(anAdjustPar - alowerboundary) : fabs(anAdjustPar - aupperboundary);
1574             Standard_Real anotherPar = (bIsFirstBoundary) ? (aupperboundary - adist) : (alowerboundary + adist);
1575             anotherPar += anoffset;
1576             Standard_Integer aneighbourpointindex = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
1577             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex);
1578             Standard_Real nU1, nV1;
1579
1580             if(surfit == 0)
1581               aNeighbourPoint.ParametersOnS1(nU1, nV1);
1582             else
1583               aNeighbourPoint.ParametersOnS2(nU1, nV1);
1584
1585             Standard_Real adist1 = (bIsUBoundary) ? fabs(nU1 - U) : fabs(nV1 - V);
1586             Standard_Real adist2 = (bIsUBoundary) ? fabs(nU1 - anotherPar) : fabs(nV1 - anotherPar);
1587             bComputeLineEnd = Standard_True;
1588             Standard_Boolean bCheckAngle1 = Standard_False;
1589             Standard_Boolean bCheckAngle2 = Standard_False;
1590             gp_Vec2d aNewVec;
1591             Standard_Real anewU = (bIsUBoundary) ? anotherPar : U;
1592             Standard_Real anewV = (bIsUBoundary) ? V : anotherPar;
1593             //
1594             if(((adist1 - adist2) > Precision::PConfusion()) &&
1595                (adist2 < (aPeriod / 4.))) {
1596               bCheckAngle1 = Standard_True;
1597               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(anewU, anewV));
1598
1599               if(aNewVec.SquareMagnitude() < (gp::Resolution() * gp::Resolution())) {
1600                 aNewP.SetValue((surfit == 0), anewU, anewV);
1601                 bCheckAngle1 = Standard_False;
1602               }
1603             }
1604             else if(adist1 < (aPeriod / 4.)) {
1605               bCheckAngle2 = Standard_True;
1606               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(U, V));
1607
1608               if(aNewVec.SquareMagnitude() < (gp::Resolution() * gp::Resolution())) {
1609                 bCheckAngle2 = Standard_False;
1610               }
1611             }
1612             //
1613             if(bCheckAngle1 || bCheckAngle2) {
1614               // assume there are at least two points in line (see "if" above)
1615               Standard_Integer anindexother = aneighbourpointindex;
1616
1617               while((anindexother <= aListOfIndex.Last()) && (anindexother >= aListOfIndex.First())) {
1618                 anindexother = (j == 0) ? (anindexother + 1) : (anindexother - 1);
1619                 const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(anindexother);
1620                 Standard_Real nU2, nV2;
1621
1622                 if(surfit == 0)
1623                   aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
1624                 else
1625                   aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
1626                 gp_Vec2d aVecOld(gp_Pnt2d(nU2, nV2), gp_Pnt2d(nU1, nV1));
1627
1628                 if(aVecOld.SquareMagnitude() <= (gp::Resolution() * gp::Resolution())) {
1629                   continue;
1630                 }
1631                 else {
1632                   Standard_Real anAngle = aNewVec.Angle(aVecOld);
1633
1634                   if((fabs(anAngle) < (M_PI * 0.25)) && (aNewVec.Dot(aVecOld) > 0.)) {
1635
1636                     if(bCheckAngle1) {
1637                       Standard_Real U1, U2, V1, V2;
1638                       IntSurf_PntOn2S atmppoint = aNewP;
1639                       atmppoint.SetValue((surfit == 0), anewU, anewV);
1640                       atmppoint.Parameters(U1, V1, U2, V2);
1641                       gp_Pnt P1 = theSurface1->Value(U1, V1);
1642                       gp_Pnt P2 = theSurface2->Value(U2, V2);
1643                       gp_Pnt P0 = aPoint.Value();
1644
1645                       if(P0.IsEqual(P1, aTol) &&
1646                          P0.IsEqual(P2, aTol) &&
1647                          P1.IsEqual(P2, aTol)) {
1648                         bComputeLineEnd = Standard_False;
1649                         aNewP.SetValue((surfit == 0), anewU, anewV);
1650                       }
1651                     }
1652
1653                     if(bCheckAngle2) {
1654                       bComputeLineEnd = Standard_False;
1655                     }
1656                   }
1657                   break;
1658                 }
1659               } // end while(anindexother...)
1660             }
1661           }
1662         }
1663         //
1664         if(bComputeLineEnd) {
1665           Standard_Integer aneighbourpointindex1 = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
1666           const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
1667           Standard_Real nU1, nV1;
1668
1669           if(surfit == 0)
1670             aNeighbourPoint.ParametersOnS1(nU1, nV1);
1671           else
1672             aNeighbourPoint.ParametersOnS2(nU1, nV1);
1673           gp_Pnt2d ap1(nU1, nV1);
1674           gp_Pnt2d ap2(nU1, nV1);
1675           Standard_Integer aneighbourpointindex2 = aneighbourpointindex1;
1676
1677           while((aneighbourpointindex2 <= aListOfIndex.Last()) && (aneighbourpointindex2 >= aListOfIndex.First())) {
1678             aneighbourpointindex2 = (j == 0) ? (aneighbourpointindex2 + 1) : (aneighbourpointindex2 - 1);
1679             const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(aneighbourpointindex2);
1680             Standard_Real nU2, nV2;
1681
1682             if(surfit == 0)
1683               aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
1684             else
1685               aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
1686             ap2.SetX(nU2);
1687             ap2.SetY(nV2);
1688
1689             if(ap1.SquareDistance(ap2) > (gp::Resolution() * gp::Resolution())) {
1690               break;
1691             }
1692           }
1693           gp_Pnt2d anewpoint;
1694           Standard_Boolean found = FindPoint(ap2, ap1, umin, umax, vmin, vmax, anewpoint);
1695
1696           if(found) {
1697             // check point
1698             Standard_Real aCriteria =aTolSum;// BRep_Tool::Tolerance(theFace1) + BRep_Tool::Tolerance(theFace2);
1699             //GeomAPI_ProjectPointOnSurf& aProjector = (surfit == 0) ? aPrj2 : aPrj1;
1700             ProjectPointOnSurf& aProjector = (surfit == 0) ? aPrj2 : aPrj1;
1701             Handle(GeomAdaptor_HSurface) aSurface = (surfit == 0) ? theSurface1 : theSurface2;
1702
1703             gp_Pnt aP3d = aSurface->Value(anewpoint.X(), anewpoint.Y());
1704             aProjector.Perform(aP3d);
1705
1706             if(aProjector.IsDone()) {
1707               if(aProjector.LowerDistance() < aCriteria) {
1708                 Standard_Real foundU = U, foundV = V;
1709                 aProjector.LowerDistanceParameters(foundU, foundV);
1710
1711                 if(surfit == 0)
1712                   aNewP.SetValue(aP3d, anewpoint.X(), anewpoint.Y(), foundU, foundV);
1713                 else
1714                   aNewP.SetValue(aP3d, foundU, foundV, anewpoint.X(), anewpoint.Y());
1715               }
1716             }
1717           }
1718         }
1719       }
1720       aSeqOfPntOn2S->Add(aNewP);
1721       aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
1722     }
1723     anArrayOfLineEnds.SetValue(i, aListOfFLIndex);
1724   }
1725   // Correct wlines.end
1726
1727   // Split wlines.begin
1728   for(j = 1; j <= theLConstructor.NbParts(); j++) {
1729     Standard_Real fprm = 0., lprm = 0.;
1730     theLConstructor.Part(j, fprm, lprm);
1731     Standard_Integer ifprm = (Standard_Integer)fprm;
1732     Standard_Integer ilprm = (Standard_Integer)lprm;
1733     //
1734     Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
1735     //
1736     for(i = 1; i <= nblines; i++) {
1737       if(anArrayOfLineType[i] != 0) {
1738         continue;
1739       }
1740       const ListOfInteger& aListOfIndex = anArrayOfLines[i];
1741
1742       if(aListOfIndex.Extent() < 2) {
1743         continue;
1744       }
1745       const TColStd_ListOfInteger& aListOfFLIndex = anArrayOfLineEnds.Value(i);
1746       Standard_Boolean bhasfirstpoint = (aListOfFLIndex.Extent() == 2);
1747       Standard_Boolean bhaslastpoint = (aListOfFLIndex.Extent() == 2);
1748
1749       if(!bhasfirstpoint && !aListOfFLIndex.IsEmpty()) {
1750         bhasfirstpoint = (i != 1);
1751       }
1752
1753       if(!bhaslastpoint && !aListOfFLIndex.IsEmpty()) {
1754         bhaslastpoint = (i != nblines);
1755       }
1756       Standard_Boolean bIsFirstInside = ((ifprm >= aListOfIndex.First()) && (ifprm <= aListOfIndex.Last()));
1757       Standard_Boolean bIsLastInside =  ((ilprm >= aListOfIndex.First()) && (ilprm <= aListOfIndex.Last()));
1758
1759       if(!bIsFirstInside && !bIsLastInside) {
1760         if((ifprm < aListOfIndex.First()) && (ilprm > aListOfIndex.Last())) {
1761           // append whole line, and boundaries if neccesary
1762           if(bhasfirstpoint) {
1763             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.First());
1764             aLineOn2S->Add(aP);
1765           }
1766           ListOfInteger::Iterator anIt(aListOfIndex);
1767
1768           for(; anIt.More(); anIt.Next()) {
1769             const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
1770             aLineOn2S->Add(aP);
1771           }
1772
1773           if(bhaslastpoint) {
1774             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.Last());
1775             aLineOn2S->Add(aP);
1776           }
1777
1778           // check end of split line (end is almost always)
1779           Standard_Integer aneighbour = i + 1;
1780           Standard_Boolean bIsEndOfLine = Standard_True;
1781
1782           if(aneighbour <= nblines) {
1783             const ListOfInteger& aListOfNeighbourIndex = anArrayOfLines[aneighbour];
1784
1785             if((anArrayOfLineType[aneighbour] != 0) &&
1786                (aListOfNeighbourIndex.IsEmpty())) {
1787               bIsEndOfLine = Standard_False;
1788             }
1789           }
1790
1791           if(bIsEndOfLine) {
1792             if(aLineOn2S->NbPoints() > 1) {
1793               Handle(IntPatch_WLine) aNewWLine =
1794                 new IntPatch_WLine(aLineOn2S, Standard_False);
1795               theNewLines.Append(aNewWLine);
1796             }
1797             aLineOn2S = new IntSurf_LineOn2S();
1798           }
1799         }
1800         continue;
1801       }
1802       // end if(!bIsFirstInside && !bIsLastInside)
1803
1804       if(bIsFirstInside && bIsLastInside) {
1805         // append inside points between ifprm and ilprm
1806         ListOfInteger::Iterator anIt(aListOfIndex);
1807
1808         for(; anIt.More(); anIt.Next()) {
1809           if((anIt.Value() < ifprm) || (anIt.Value() > ilprm))
1810             continue;
1811           const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
1812           aLineOn2S->Add(aP);
1813         }
1814       }
1815       else {
1816
1817         if(bIsFirstInside) {
1818           // append points from ifprm to last point + boundary point
1819           ListOfInteger::Iterator anIt(aListOfIndex);
1820
1821           for(; anIt.More(); anIt.Next()) {
1822             if(anIt.Value() < ifprm)
1823               continue;
1824             const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
1825             aLineOn2S->Add(aP);
1826           }
1827
1828           if(bhaslastpoint) {
1829             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.Last());
1830             aLineOn2S->Add(aP);
1831           }
1832           // check end of split line (end is almost always)
1833           Standard_Integer aneighbour = i + 1;
1834           Standard_Boolean bIsEndOfLine = Standard_True;
1835
1836           if(aneighbour <= nblines) {
1837             const ListOfInteger& aListOfNeighbourIndex = anArrayOfLines[aneighbour];
1838
1839             if((anArrayOfLineType[aneighbour] != 0) &&
1840                (aListOfNeighbourIndex.IsEmpty())) {
1841               bIsEndOfLine = Standard_False;
1842             }
1843           }
1844
1845           if(bIsEndOfLine) {
1846             if(aLineOn2S->NbPoints() > 1) {
1847               Handle(IntPatch_WLine) aNewWLine =
1848                 new IntPatch_WLine(aLineOn2S, Standard_False);
1849               theNewLines.Append(aNewWLine);
1850             }
1851             aLineOn2S = new IntSurf_LineOn2S();
1852           }
1853         }
1854         // end if(bIsFirstInside)
1855
1856         if(bIsLastInside) {
1857           // append points from first boundary point to ilprm
1858           if(bhasfirstpoint) {
1859             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.First());
1860             aLineOn2S->Add(aP);
1861           }
1862           ListOfInteger::Iterator anIt(aListOfIndex);
1863
1864           for(; anIt.More(); anIt.Next()) {
1865             if(anIt.Value() > ilprm)
1866               continue;
1867             const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
1868             aLineOn2S->Add(aP);
1869           }
1870         }
1871         //end if(bIsLastInside)
1872       }
1873     }
1874
1875     if(aLineOn2S->NbPoints() > 1) {
1876       Handle(IntPatch_WLine) aNewWLine =
1877         new IntPatch_WLine(aLineOn2S, Standard_False);
1878       theNewLines.Append(aNewWLine);
1879     }
1880   }
1881   // Split wlines.end
1882   //
1883   // cda002/I3
1884   Standard_Real fprm, lprm;
1885   Standard_Integer ifprm, ilprm, aNbPoints, aIndex;
1886   //
1887   aNbParts=theLConstructor.NbParts();
1888   //
1889   for(j = 1; j <= aNbParts; j++) {
1890     theLConstructor.Part(j, fprm, lprm);
1891     ifprm=(Standard_Integer)fprm;
1892     ilprm=(Standard_Integer)lprm;
1893     //
1894     if ((ilprm-ifprm)==1) {
1895       for(i = 1; i <= nblines; i++) {
1896         aLineType=anArrayOfLineType[i];
1897         if(aLineType) {
1898           continue;
1899         }
1900         //
1901         const ListOfInteger& aListOfIndex = anArrayOfLines[i];
1902         aNbPoints=aListOfIndex.Extent();
1903         if(aNbPoints==1) {
1904           aIndex=aListOfIndex.First();
1905           if (aIndex==ifprm || aIndex==ilprm) {
1906             Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
1907             const IntSurf_PntOn2S& aP1 = theWLine->Point(ifprm);
1908             const IntSurf_PntOn2S& aP2 = theWLine->Point(ilprm);
1909             aLineOn2S->Add(aP1);
1910             aLineOn2S->Add(aP2);
1911             Handle(IntPatch_WLine) aNewWLine =
1912               new IntPatch_WLine(aLineOn2S, Standard_False);
1913             theNewLines.Append(aNewWLine);
1914           }
1915         }
1916       }
1917     }
1918   }
1919   //
1920   return Standard_True;
1921 }
1922
1923 //=======================================================================
1924 //function : AdjustPeriodic
1925 //purpose  :
1926 //=======================================================================
1927 Standard_Real AdjustPeriodic(const Standard_Real theParameter,
1928                              const Standard_Real parmin,
1929                              const Standard_Real parmax,
1930                              const Standard_Real thePeriod,
1931                              Standard_Real&      theOffset)
1932 {
1933   Standard_Real aresult = theParameter;
1934   theOffset = 0.;
1935   while(aresult < parmin) {
1936     aresult += thePeriod;
1937     theOffset += thePeriod;
1938   }
1939   while(aresult > parmax) {
1940     aresult -= thePeriod;
1941     theOffset -= thePeriod;
1942   }
1943   return aresult;
1944 }
1945
1946 //=======================================================================
1947 //function : IsPointOnBoundary
1948 //purpose  :
1949 //=======================================================================
1950 Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
1951                                    const Standard_Real theFirstBoundary,
1952                                    const Standard_Real theSecondBoundary,
1953                                    const Standard_Real theResolution,
1954                                    Standard_Boolean&   IsOnFirstBoundary)
1955 {
1956   IsOnFirstBoundary = Standard_True;
1957   if(fabs(theParameter - theFirstBoundary) < theResolution)
1958     return Standard_True;
1959   if(fabs(theParameter - theSecondBoundary) < theResolution)
1960   {
1961     IsOnFirstBoundary = Standard_False;
1962     return Standard_True;
1963   }
1964   return Standard_False;
1965 }
1966 //=======================================================================
1967 //function : FindPoint
1968 //purpose  :
1969 //=======================================================================
1970 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
1971                            const gp_Pnt2d&     theLastPoint,
1972                            const Standard_Real theUmin,
1973                            const Standard_Real theUmax,
1974                            const Standard_Real theVmin,
1975                            const Standard_Real theVmax,
1976                            gp_Pnt2d&           theNewPoint)
1977 {
1978   gp_Vec2d aVec(theFirstPoint, theLastPoint);
1979   Standard_Integer i = 0, j = 0;
1980
1981   for(i = 0; i < 4; i++) {
1982     gp_Vec2d anOtherVec;
1983     gp_Vec2d anOtherVecNormal;
1984     gp_Pnt2d aprojpoint = theLastPoint;
1985
1986     if((i % 2) == 0) {
1987       anOtherVec.SetX(0.);
1988       anOtherVec.SetY(1.);
1989       anOtherVecNormal.SetX(1.);
1990       anOtherVecNormal.SetY(0.);
1991
1992       if(i < 2)
1993         aprojpoint.SetX(theUmin);
1994       else
1995         aprojpoint.SetX(theUmax);
1996     }
1997     else {
1998       anOtherVec.SetX(1.);
1999       anOtherVec.SetY(0.);
2000       anOtherVecNormal.SetX(0.);
2001       anOtherVecNormal.SetY(1.);
2002
2003       if(i < 2)
2004         aprojpoint.SetY(theVmin);
2005       else
2006         aprojpoint.SetY(theVmax);
2007     }
2008     gp_Vec2d anormvec = aVec;
2009     anormvec.Normalize();
2010     Standard_Real adot1 = anormvec.Dot(anOtherVecNormal);
2011
2012     if(fabs(adot1) < Precision::Angular())
2013       continue;
2014     Standard_Real adist = 0.;
2015
2016     if((i % 2) == 0) {
2017       adist = (i < 2) ? fabs(theLastPoint.X() - theUmin) : fabs(theLastPoint.X() - theUmax);
2018     }
2019     else {
2020       adist = (i < 2) ? fabs(theLastPoint.Y() - theVmin) : fabs(theLastPoint.Y() - theVmax);
2021     }
2022     Standard_Real anoffset = adist * anOtherVec.Dot(anormvec) / adot1;
2023
2024     for(j = 0; j < 2; j++) {
2025       anoffset = (j == 0) ? anoffset : -anoffset;
2026       gp_Pnt2d acurpoint(aprojpoint.XY() + (anOtherVec.XY()*anoffset));
2027       gp_Vec2d acurvec(theLastPoint, acurpoint);
2028
2029       //
2030       Standard_Real aDotX, anAngleX, aPC;
2031       //
2032       aDotX=aVec.Dot(acurvec);
2033       anAngleX=aVec.Angle(acurvec);
2034       aPC=Precision::PConfusion();
2035       //
2036       if(aDotX > 0. && fabs(anAngleX) < aPC) {
2037       //
2038         if((i % 2) == 0) {
2039           if((acurpoint.Y() >= theVmin) &&
2040              (acurpoint.Y() <= theVmax)) {
2041             theNewPoint = acurpoint;
2042             return Standard_True;
2043           }
2044         }
2045         else {
2046           if((acurpoint.X() >= theUmin) &&
2047              (acurpoint.X() <= theUmax)) {
2048             theNewPoint = acurpoint;
2049             return Standard_True;
2050           }
2051         }
2052       }
2053     }
2054   }
2055   return Standard_False;
2056 }
2057
2058 //=======================================================================
2059 //function : ParametersOfNearestPointOnSurface
2060 //purpose  :
2061 //=======================================================================
2062 static Standard_Boolean ParametersOfNearestPointOnSurface(const Extrema_ExtPS theExtr,
2063                                                           Standard_Real& theU,
2064                                                           Standard_Real& theV)
2065 {
2066   if(!theExtr.IsDone() || !theExtr.NbExt())
2067     return Standard_False;
2068
2069   Standard_Integer anIndex = 1;
2070   Standard_Real aMinSQDist = theExtr.SquareDistance(anIndex);
2071   for(Standard_Integer i = 2; i <= theExtr.NbExt(); i++)
2072   {
2073     Standard_Real aSQD = theExtr.SquareDistance(i);
2074     if (aSQD < aMinSQDist)
2075     {
2076       aMinSQDist = aSQD;
2077       anIndex = i;
2078     }
2079   }
2080
2081   theExtr.Point(anIndex).Parameter(theU, theV);
2082
2083   return Standard_True;
2084 }
2085
2086 //=======================================================================
2087 //function : GetSegmentBoundary
2088 //purpose  :
2089 //=======================================================================
2090 static void GetSegmentBoundary( const IntRes2d_IntersectionSegment& theSegm,
2091                                 const Handle(Geom2d_Curve)& theCurve,
2092                                 GeomInt_VectorOfReal& theArrayOfParameters)
2093 {
2094   Standard_Real aU1 = theCurve->FirstParameter(), aU2 = theCurve->LastParameter();
2095
2096   if(theSegm.HasFirstPoint())
2097   {
2098     const IntRes2d_IntersectionPoint& anIPF = theSegm.FirstPoint();
2099     aU1 = anIPF.ParamOnFirst();
2100   }
2101
2102   if(theSegm.HasLastPoint())
2103   {
2104     const IntRes2d_IntersectionPoint& anIPL = theSegm.LastPoint();
2105     aU2 = anIPL.ParamOnFirst();
2106   }
2107
2108   theArrayOfParameters.Append(aU1);
2109   theArrayOfParameters.Append(aU2);
2110 }
2111
2112 //=======================================================================
2113 //function : IntersectCurveAndBoundary
2114 //purpose  :
2115 //=======================================================================
2116 static void IntersectCurveAndBoundary(const Handle(Geom2d_Curve)& theC2d,
2117                                       const Handle(Geom2d_Curve)* const theArrBounds,
2118                                       const Standard_Integer theNumberOfCurves,
2119                                       const Standard_Real theTol,
2120                                       GeomInt_VectorOfReal& theArrayOfParameters)
2121 {
2122   if(theC2d.IsNull())
2123     return;
2124
2125   Geom2dAdaptor_Curve anAC1(theC2d);
2126   for(Standard_Integer aCurID = 0; aCurID < theNumberOfCurves; aCurID++)
2127   {
2128     if(theArrBounds[aCurID].IsNull())
2129       continue;
2130
2131     Geom2dAdaptor_Curve anAC2(theArrBounds[aCurID]);
2132     Geom2dInt_GInter anIntCC2d(anAC1, anAC2, theTol, theTol);
2133
2134     if(!anIntCC2d.IsDone() || anIntCC2d.IsEmpty())
2135       continue;
2136
2137     for (Standard_Integer aPntID = 1; aPntID <= anIntCC2d.NbPoints(); aPntID++)
2138     {
2139       const Standard_Real aParam = anIntCC2d.Point(aPntID).ParamOnFirst();
2140       theArrayOfParameters.Append(aParam);
2141     }
2142
2143     for (Standard_Integer aSegmID = 1; aSegmID <= anIntCC2d.NbSegments(); aSegmID++)
2144     {
2145       GetSegmentBoundary(anIntCC2d.Segment(aSegmID), theC2d, theArrayOfParameters);
2146     }
2147   }
2148 }
2149
2150 //=======================================================================
2151 //function : TreatRLine
2152 //purpose  : Approx of Restriction line
2153 //=======================================================================
2154 void GeomInt_IntSS::TreatRLine(const Handle(IntPatch_RLine)& theRL,
2155                 const Handle(GeomAdaptor_HSurface)& theHS1,
2156                 const Handle(GeomAdaptor_HSurface)& theHS2,
2157                 Handle(Geom_Curve)& theC3d,
2158                 Handle(Geom2d_Curve)& theC2d1,
2159                 Handle(Geom2d_Curve)& theC2d2,
2160                 Standard_Real& theTolReached)
2161 {
2162   Handle(GeomAdaptor_HSurface) aGAHS;
2163   Handle(Adaptor2d_HCurve2d) anAHC2d;
2164   Standard_Real tf, tl;
2165   gp_Lin2d aL;
2166   // It is assumed that 2d curve is 2d line (rectangular surface domain)
2167   if(theRL->IsArcOnS1())
2168   {
2169     aGAHS = theHS1;
2170     anAHC2d = theRL->ArcOnS1();
2171     theRL->ParamOnS1(tf, tl);
2172     theC2d1 = Geom2dAdaptor::MakeCurve(anAHC2d->Curve2d());
2173   }
2174   else if (theRL->IsArcOnS2())
2175   {
2176     aGAHS = theHS2;
2177     anAHC2d = theRL->ArcOnS2();
2178     theRL->ParamOnS2(tf, tl);
2179     theC2d2 = Geom2dAdaptor::MakeCurve(anAHC2d->Curve2d());
2180   }
2181   else
2182   {
2183     return;
2184   }
2185   //
2186   //To provide sameparameter it is necessary to get 3d curve as
2187   //approximation of curve on surface.
2188   Standard_Integer aMaxDeg = 8;
2189   Standard_Integer aMaxSeg = 1000;
2190   Approx_CurveOnSurface anApp(anAHC2d, aGAHS, tf, tl, Precision::Confusion(),
2191                               GeomAbs_C1, aMaxDeg, aMaxSeg,
2192                               Standard_True, Standard_False);
2193   if(!anApp.HasResult())
2194     return;
2195
2196   theC3d = anApp.Curve3d();
2197   theTolReached = anApp.MaxError3d();
2198   Standard_Real aTol = Precision::Confusion();
2199   if(theRL->IsArcOnS1())
2200   {
2201     Handle(Geom_Surface) aS = GeomAdaptor::MakeSurface(theHS2->Surface());
2202     BuildPCurves (tf, tl, aTol,
2203                   aS, theC3d, theC2d2);
2204   }
2205   if(theRL->IsArcOnS2())
2206   {
2207     Handle(Geom_Surface) aS = GeomAdaptor::MakeSurface(theHS1->Surface());
2208     BuildPCurves (tf, tl, aTol,
2209                   aS, theC3d, theC2d1);
2210   }
2211   theTolReached = Max(theTolReached, aTol);
2212 }
2213
2214 //=======================================================================
2215 //function : BuildPCurves
2216 //purpose  :
2217 //=======================================================================
2218 void GeomInt_IntSS::BuildPCurves (Standard_Real f,
2219                                   Standard_Real l,
2220                                   Standard_Real& Tol,
2221                                   const Handle (Geom_Surface)& S,
2222                                   const Handle (Geom_Curve)&   C,
2223                                   Handle (Geom2d_Curve)& C2d)
2224 {
2225   if (!C2d.IsNull()) {
2226     return;
2227   }
2228   //
2229   Standard_Real umin,umax,vmin,vmax;
2230   //
2231   S->Bounds(umin, umax, vmin, vmax);
2232   // in class ProjLib_Function the range of parameters is shrank by 1.e-09
2233   if((l - f) > 2.e-09) {
2234     C2d = GeomProjLib::Curve2d(C,f,l,S,umin,umax,vmin,vmax,Tol);
2235     //
2236     if (C2d.IsNull()) {
2237       // proj. a circle that goes through the pole on a sphere to the sphere
2238       Tol += Precision::Confusion();
2239       C2d = GeomProjLib::Curve2d(C,f,l,S,Tol);
2240     }
2241   }
2242   else {
2243     if((l - f) > Epsilon(Abs(f)))
2244     {
2245       //The domain of C2d is [Epsilon(Abs(f)), 2.e-09]
2246       //On this small range C2d can be considered as segment
2247       //of line.
2248
2249       Standard_Real aU=0., aV=0.;
2250       GeomAdaptor_Surface anAS;
2251       anAS.Load(S);
2252       Extrema_ExtPS anExtr;
2253       const gp_Pnt aP3d1 = C->Value(f);
2254       const gp_Pnt aP3d2 = C->Value(l);
2255
2256       anExtr.SetAlgo(Extrema_ExtAlgo_Grad);
2257       anExtr.Initialize(anAS, umin, umax, vmin, vmax,
2258                                 Precision::Confusion(), Precision::Confusion());
2259       anExtr.Perform(aP3d1);
2260
2261       if(ParametersOfNearestPointOnSurface(anExtr, aU, aV))
2262       {
2263         const gp_Pnt2d aP2d1(aU, aV);
2264
2265         anExtr.Perform(aP3d2);
2266
2267         if(ParametersOfNearestPointOnSurface(anExtr, aU, aV))
2268         {
2269           const gp_Pnt2d aP2d2(aU, aV);
2270
2271           if(aP2d1.Distance(aP2d2) > gp::Resolution())
2272           {
2273             TColgp_Array1OfPnt2d poles(1,2);
2274             TColStd_Array1OfReal knots(1,2);
2275             TColStd_Array1OfInteger mults(1,2);
2276             poles(1) = aP2d1;
2277             poles(2) = aP2d2;
2278             knots(1) = f;
2279             knots(2) = l;
2280             mults(1) = mults(2) = 2;
2281
2282             C2d = new Geom2d_BSplineCurve(poles,knots,mults,1);
2283
2284             //Check same parameter in middle point .begin
2285             const gp_Pnt PMid(C->Value(0.5*(f+l)));
2286             const gp_Pnt2d pmidcurve2d(0.5*(aP2d1.XY() + aP2d2.XY()));
2287             const gp_Pnt aPC(anAS.Value(pmidcurve2d.X(), pmidcurve2d.Y()));
2288             const Standard_Real aDist = PMid.Distance(aPC);
2289             Tol = Max(aDist, Tol);
2290             //Check same parameter in middle point .end
2291           }
2292         }
2293       }
2294     }
2295   }
2296   //
2297   if (S->IsUPeriodic() && !C2d.IsNull()) {
2298     // Recadre dans le domaine UV de la face
2299     Standard_Real aTm, U0, aEps, period, du, U0x;
2300     Standard_Boolean bAdjust;
2301     //
2302     aEps = Precision::PConfusion();
2303     period = S->UPeriod();
2304     //
2305     aTm = .5*(f + l);
2306     gp_Pnt2d pm = C2d->Value(aTm);
2307     U0 = pm.X();
2308     //
2309     bAdjust =
2310       GeomInt::AdjustPeriodic(U0, umin, umax, period, U0x, du, aEps);
2311     if (bAdjust) {
2312       gp_Vec2d T1(du, 0.);
2313       C2d->Translate(T1);
2314     }
2315   }
2316 }
2317
2318 //=======================================================================
2319 //function : TrimILineOnSurfBoundaries
2320 //purpose  : This function finds intersection points of given curves with
2321 //            surface boundaries and fills theArrayOfParameters by parameters
2322 //            along the given curves corresponding of these points.
2323 //=======================================================================
2324 void GeomInt_IntSS::TrimILineOnSurfBoundaries(const Handle(Geom2d_Curve)& theC2d1,
2325                                               const Handle(Geom2d_Curve)& theC2d2,
2326                                               const Bnd_Box2d& theBound1,
2327                                               const Bnd_Box2d& theBound2,
2328                                               GeomInt_VectorOfReal& theArrayOfParameters)
2329 {
2330   //Rectangular boundaries of two surfaces: [0]:U=Ufirst, [1]:U=Ulast,
2331   //                                        [2]:V=Vfirst, [3]:V=Vlast
2332   const Standard_Integer aNumberOfCurves = 4;
2333   Handle(Geom2d_Curve) aCurS1Bounds[aNumberOfCurves];
2334   Handle(Geom2d_Curve) aCurS2Bounds[aNumberOfCurves];
2335
2336   Standard_Real aU1f=0.0, aV1f=0.0, aU1l=0.0, aV1l=0.0;
2337   Standard_Real aU2f=0.0, aV2f=0.0, aU2l=0.0, aV2l=0.0;
2338
2339   theBound1.Get(aU1f, aV1f, aU1l, aV1l);
2340   theBound2.Get(aU2f, aV2f, aU2l, aV2l);
2341
2342   Standard_Real aDelta = aV1l-aV1f;
2343   if(Abs(aDelta) > RealSmall())
2344   {
2345     if(!Precision::IsInfinite(aU1f))
2346     {
2347       aCurS1Bounds[0] = new Geom2d_Line(gp_Pnt2d(aU1f, aV1f), gp_Dir2d(0.0, 1.0));
2348
2349       if(!Precision::IsInfinite(aDelta))
2350         aCurS1Bounds[0] = new Geom2d_TrimmedCurve(aCurS1Bounds[0], 0, aDelta);
2351     }
2352
2353     if(!Precision::IsInfinite(aU1l))
2354     {
2355       aCurS1Bounds[1] = new Geom2d_Line(gp_Pnt2d(aU1l, aV1f), gp_Dir2d(0.0, 1.0));
2356       if(!Precision::IsInfinite(aDelta))
2357         aCurS1Bounds[1] = new Geom2d_TrimmedCurve(aCurS1Bounds[1], 0, aDelta);
2358     }
2359   }
2360
2361   aDelta = aU1l-aU1f;
2362   if(Abs(aDelta) > RealSmall())
2363   {
2364     if(!Precision::IsInfinite(aV1f))
2365     {
2366       aCurS1Bounds[2] = new Geom2d_Line(gp_Pnt2d(aU1f, aV1f), gp_Dir2d(1.0, 0.0));
2367       if(!Precision::IsInfinite(aDelta))
2368         aCurS1Bounds[2] = new Geom2d_TrimmedCurve(aCurS1Bounds[2], 0, aDelta);
2369     }
2370
2371     if(!Precision::IsInfinite(aV1l))
2372     {
2373       aCurS1Bounds[3] = new Geom2d_Line(gp_Pnt2d(aU1l, aV1l), gp_Dir2d(1.0, 0.0));
2374       if(!Precision::IsInfinite(aDelta))
2375         aCurS1Bounds[3] = new Geom2d_TrimmedCurve(aCurS1Bounds[3], 0, aDelta);
2376     }
2377   }
2378
2379   aDelta = aV2l-aV2f;
2380   if(Abs(aDelta) > RealSmall())
2381   {
2382     if(!Precision::IsInfinite(aU2f))
2383     {
2384       aCurS2Bounds[0] = new Geom2d_Line(gp_Pnt2d(aU2f, aV2f), gp_Dir2d(0.0, 1.0));
2385       if(!Precision::IsInfinite(aDelta))
2386         aCurS2Bounds[0] = new Geom2d_TrimmedCurve(aCurS2Bounds[0], 0, aDelta);
2387     }
2388
2389     if(!Precision::IsInfinite(aU2l))
2390     {
2391       aCurS2Bounds[1] = new Geom2d_Line(gp_Pnt2d(aU2l, aV2f), gp_Dir2d(0.0, 1.0));
2392       if(!Precision::IsInfinite(aDelta))
2393         aCurS2Bounds[1] = new Geom2d_TrimmedCurve(aCurS2Bounds[1], 0, aDelta);
2394     }
2395   }
2396
2397   aDelta = aU2l-aU2f;
2398   if(Abs(aDelta) > RealSmall())
2399   {
2400     if(!Precision::IsInfinite(aV2f))
2401     {
2402       aCurS2Bounds[2] = new Geom2d_Line(gp_Pnt2d(aU2f, aV2f), gp_Dir2d(1.0, 0.0));
2403       if(!Precision::IsInfinite(aDelta))
2404         aCurS2Bounds[2] = new Geom2d_TrimmedCurve(aCurS2Bounds[2], 0, aDelta);
2405     }
2406
2407     if(!Precision::IsInfinite(aV2l))
2408     {
2409       aCurS2Bounds[3] = new Geom2d_Line(gp_Pnt2d(aU2l, aV2l), gp_Dir2d(1.0, 0.0));
2410       if(!Precision::IsInfinite(aDelta))
2411         aCurS2Bounds[3] = new Geom2d_TrimmedCurve(aCurS2Bounds[3], 0, aDelta);
2412     }
2413   }
2414
2415   const Standard_Real anIntTol = 10.0*Precision::Confusion();
2416
2417   IntersectCurveAndBoundary(theC2d1, aCurS1Bounds,
2418                         aNumberOfCurves, anIntTol, theArrayOfParameters);
2419
2420   IntersectCurveAndBoundary(theC2d2, aCurS2Bounds,
2421                         aNumberOfCurves, anIntTol, theArrayOfParameters);
2422
2423   std::sort(theArrayOfParameters.begin(), theArrayOfParameters.end());
2424 }