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