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
   9 // under the terms of the GNU Lesser General Public version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 #include <GeomInt_IntSS.ixx>
  18
  19 #include <Precision.hxx>
  20 #include <gp_Pnt.hxx>
  21 #include <gp_Pnt2d.hxx>
  22 #include <gp_Pln.hxx>
  23
  24 #include <TColStd_Array1OfReal.hxx>
  25 #include <TColStd_Array1OfInteger.hxx>
  26 #include <TColStd_Array1OfListOfInteger.hxx>
  27 #include <TColStd_SequenceOfReal.hxx>
  28 #include <TColStd_ListOfInteger.hxx>
  29 #include <TColStd_ListIteratorOfListOfInteger.hxx>
  30
  31 #include <TColgp_Array1OfPnt.hxx>
  32 #include <TColgp_Array1OfPnt2d.hxx>
  33
  34 #include <Adaptor3d_TopolTool.hxx>
  35 #include <IntPatch_Line.hxx>
  36 #include <IntPatch_WLine.hxx>
  37 #include <IntPatch_GLine.hxx>
  38 #include <IntPatch_ALineToWLine.hxx>
  39 #include <IntPatch_IType.hxx>
  40 #include <NCollection_IncAllocator.hxx>
  41 #include <NCollection_List.hxx>
  42 #include <NCollection_LocalArray.hxx>
  43 #include <NCollection_StdAllocator.hxx>
  44 #include <vector>
  45
  46 #include <AppParCurves_MultiBSpCurve.hxx>
  47
  48 #include <GeomAbs_SurfaceType.hxx>
  49 #include <GeomAbs_CurveType.hxx>
  50
  51 #include <Geom_Curve.hxx>
  52 #include <Geom_Line.hxx>
  53 #include <Geom_Parabola.hxx>
  54 #include <Geom_Hyperbola.hxx>
  55 #include <Geom_TrimmedCurve.hxx>
  56 #include <Geom_Circle.hxx>
  57 #include <Geom_Ellipse.hxx>
  58 #include <Geom_BSplineCurve.hxx>
  59
  60 #include <Geom2d_Curve.hxx>
  61 #include <Geom2d_BSplineCurve.hxx>
  62 #include <Geom2d_TrimmedCurve.hxx>
  63
  64 #include <GeomLib_CheckBSplineCurve.hxx>
  65 #include <GeomLib_Check2dBSplineCurve.hxx>
  66 #include <GeomProjLib.hxx>
  67
  68
  69 #include <ElSLib.hxx>
  70
  71 #include <GeomInt_WLApprox.hxx>
  72 #include <Extrema_ExtPS.hxx>
  73
  74 static
  75   void BuildPCurves (Standard_Real f,
  76                      Standard_Real l,
  77                      Standard_Real& Tol,
  78                      const Handle (Geom_Surface)& S,
  79                      const Handle (Geom_Curve)&   C,
  80                      Handle (Geom2d_Curve)& C2d);
  81
  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
 145 //
 146
 147
 148 class ProjectPointOnSurf
 149 {
 150  public:
 151   ProjectPointOnSurf() : myIsDone (Standard_False),myIndex(0) {}
 152   void Init(const Handle(Geom_Surface)& Surface,
 153             const Standard_Real Umin,
 154             const Standard_Real Usup,
 155             const Standard_Real Vmin,
 156             const Standard_Real Vsup);
 157   void Init ();
 158   void Perform(const gp_Pnt& P);
 159   Standard_Boolean IsDone () const { return myIsDone; }
 160   void LowerDistanceParameters(Standard_Real& U, Standard_Real& V ) const;
 161   Standard_Real LowerDistance() const;
 162  protected:
 163   Standard_Boolean myIsDone;
 164   Standard_Integer myIndex;
 165   Extrema_ExtPS myExtPS;
 166   GeomAdaptor_Surface myGeomAdaptor;
 167 };
 168
 169 //=======================================================================
 170 //function : Init
 171 //purpose  :
 172 //=======================================================================
 173   void ProjectPointOnSurf::Init (const Handle(Geom_Surface)& Surface,
 174                                  const Standard_Real       Umin,
 175                                  const Standard_Real       Usup,
 176                                  const Standard_Real       Vmin,
 177                                  const Standard_Real       Vsup )
 178 {
 179   const Standard_Real Tolerance = Precision::PConfusion();
 180   //
 181   myGeomAdaptor.Load(Surface, Umin,Usup,Vmin,Vsup);
 182   myExtPS.Initialize(myGeomAdaptor, Umin, Usup, Vmin, Vsup, Tolerance, Tolerance);
 183   myIsDone = Standard_False;
 184 }
 185 //=======================================================================
 186 //function : Init
 187 //purpose  :
 188 //=======================================================================
 189   void ProjectPointOnSurf::Init ()
 190 {
 191   myIsDone = myExtPS.IsDone() && (myExtPS.NbExt() > 0);
 192   if (myIsDone) {
 193     // evaluate the lower distance and its index;
 194     Standard_Real Dist2Min = myExtPS.SquareDistance(1);
 195     myIndex = 1;
 196     for (Standard_Integer i = 2; i <= myExtPS.NbExt(); i++)
 197     {
 198       const Standard_Real Dist2 = myExtPS.SquareDistance(i);
 199       if (Dist2 < Dist2Min) {
 200         Dist2Min = Dist2;
 201         myIndex = i;
 202       }
 203     }
 204   }
 205 }
 206 //=======================================================================
 207 //function : Perform
 208 //purpose  :
 209 //=======================================================================
 210 void ProjectPointOnSurf::Perform(const gp_Pnt& P)
 211 {
 212   myExtPS.Perform(P);
 213   Init ();
 214 }
 215 //=======================================================================
 216 //function : LowerDistanceParameters
 217 //purpose  :
 218 //=======================================================================
 219 void ProjectPointOnSurf::LowerDistanceParameters (Standard_Real& U,
 220                                                   Standard_Real& V ) const
 221 {
 222   StdFail_NotDone_Raise_if(!myIsDone, "GeomInt_IntSS::ProjectPointOnSurf::LowerDistanceParameters");
 223   (myExtPS.Point(myIndex)).Parameter(U,V);
 224 }
 225 //=======================================================================
 226 //function : LowerDistance
 227 //purpose  :
 228 //=======================================================================
 229 Standard_Real ProjectPointOnSurf::LowerDistance() const
 230 {
 231   StdFail_NotDone_Raise_if(!myIsDone, "GeomInt_IntSS::ProjectPointOnSurf::LowerDistance");
 232   return sqrt(myExtPS.SquareDistance(myIndex));
 233 }
 234 //
 235
 236 //=======================================================================
 237 //function : MakeCurve
 238 //purpose  :
 239 //=======================================================================
 240   void GeomInt_IntSS::MakeCurve(const Standard_Integer Index,
 241                                 const Handle(Adaptor3d_TopolTool) & dom1,
 242                                 const Handle(Adaptor3d_TopolTool) & dom2,
 243                                 const Standard_Real Tol,
 244                                 const Standard_Boolean Approx,
 245                                 const Standard_Boolean ApproxS1,
 246                                 const Standard_Boolean ApproxS2)
 247
 248 {
 249   Standard_Boolean myApprox1, myApprox2, myApprox;
 250   Standard_Real Tolpc, myTolApprox;
 251   IntPatch_IType typl;
 252   Handle(Geom2d_BSplineCurve) H1;
 253   Handle(Geom_Surface) aS1, aS2;
 254   //
 255   Tolpc = Tol;
 256   myApprox=Approx;
 257   myApprox1=ApproxS1;
 258   myApprox2=ApproxS2;
 259   myTolApprox=0.0000001;
 260   //
 261   aS1=myHS1->ChangeSurface().Surface();
 262   aS2=myHS2->ChangeSurface().Surface();
 263   //
 264   Handle(IntPatch_Line) L = myIntersector.Line(Index);
 265   typl = L->ArcType();
 266   //
 267   if(typl==IntPatch_Walking) {
 268     Handle(IntPatch_Line) anewL =
 269       ComputePurgedWLine(Handle(IntPatch_WLine)::DownCast(L));
 270     //
 271     if(anewL.IsNull()) {
 272       return;
 273     }
 274     L = anewL;
 275   }
 276   //
 277   // Line Constructor
 278   myLConstruct.Perform(L);
 279   if (!myLConstruct.IsDone() || myLConstruct.NbParts() <= 0) {
 280     return;
 281   }
 282   // Do the Curve
 283   Standard_Boolean ok;
 284   Standard_Integer i, j,  aNbParts;
 285   Standard_Real fprm, lprm;
 286   Handle(Geom_Curve) newc;
 287
 288   switch (typl) {
 289   //########################################
 290   // Line, Parabola, Hyperbola
 291   //########################################
 292   case IntPatch_Lin:
 293   case IntPatch_Parabola:
 294   case IntPatch_Hyperbola: {
 295     if (typl == IntPatch_Lin) {
 296       newc=new Geom_Line (Handle(IntPatch_GLine)::DownCast(L)->Line());
 297     }
 298     else if (typl == IntPatch_Parabola) {
 299       newc=new Geom_Parabola(Handle(IntPatch_GLine)::DownCast(L)->Parabola());
 300     }
 301     else if (typl == IntPatch_Hyperbola) {
 302       newc=new Geom_Hyperbola (Handle(IntPatch_GLine)::DownCast(L)->Hyperbola());
 303     }
 304     //
 305     aNbParts=myLConstruct.NbParts();
 306     for (i=1; i<=aNbParts; i++) {
 307       myLConstruct.Part(i, fprm, lprm);
 308
 309       if (!Precision::IsNegativeInfinite(fprm) &&
 310           !Precision::IsPositiveInfinite(lprm)) {
 311         // /cto/900/F1
 312         //if (typl == IntPatch_Lin) {
 313           //Standard_Real dPrm=1.3e-6;
 314           //fprm=fprm-dPrm;
 315           //lprm=lprm+dPrm;
 316         //}
 317         //
 318         Handle(Geom_TrimmedCurve) aCT3D=new Geom_TrimmedCurve(newc, fprm, lprm);
 319         sline.Append(aCT3D);
 320         //
 321         if(myApprox1) {
 322           Handle (Geom2d_Curve) C2d;
 323           BuildPCurves(fprm, lprm, Tolpc, myHS1->ChangeSurface().Surface(), newc, C2d);
 324           if(Tolpc>myTolReached2d || myTolReached2d==0.) {
 325             myTolReached2d=Tolpc;
 326           }
 327           slineS1.Append(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
 328         }
 329         else {
 330           slineS1.Append(H1);
 331         }
 332         //
 333         if(myApprox2) {
 334     Handle (Geom2d_Curve) C2d;
 335     BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
 336     if(Tolpc>myTolReached2d || myTolReached2d==0.) {
 337       myTolReached2d=Tolpc;
 338     }
 339     //
 340     slineS2.Append(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
 341   }
 342   else {
 343     slineS2.Append(H1);
 344   }
 345       } // if (!Precision::IsNegativeInfinite(fprm) &&  !Precision::IsPositiveInfinite(lprm))
 346
 347       else {
 348         GeomAbs_SurfaceType typS1 = myHS1->Surface().GetType();
 349         GeomAbs_SurfaceType typS2 = myHS2->Surface().GetType();
 350         if( typS1 == GeomAbs_SurfaceOfExtrusion ||
 351           typS1 == GeomAbs_OffsetSurface ||
 352           typS1 == GeomAbs_SurfaceOfRevolution ||
 353           typS2 == GeomAbs_SurfaceOfExtrusion ||
 354           typS2 == GeomAbs_OffsetSurface ||
 355           typS2 == GeomAbs_SurfaceOfRevolution) {
 356             sline.Append(newc);
 357             slineS1.Append(H1);
 358             slineS2.Append(H1);
 359             continue;
 360         }
 361         Standard_Boolean bFNIt, bLPIt;
 362         Standard_Real aTestPrm, dT=100.;
 363         Standard_Real u1, v1, u2, v2, TolX;
 364         //
 365         bFNIt=Precision::IsNegativeInfinite(fprm);
 366         bLPIt=Precision::IsPositiveInfinite(lprm);
 367
 368         aTestPrm=0.;
 369
 370         if (bFNIt && !bLPIt) {
 371           aTestPrm=lprm-dT;
 372         }
 373         else if (!bFNIt && bLPIt) {
 374           aTestPrm=fprm+dT;
 375         }
 376         //
 377         gp_Pnt ptref(newc->Value(aTestPrm));
 378         //
 379         TolX = Precision::Confusion();
 380         Parameters(myHS1, myHS2, ptref,  u1, v1, u2, v2);
 381         ok = (dom1->Classify(gp_Pnt2d(u1, v1), TolX) != TopAbs_OUT);
 382         if(ok) {
 383           ok = (dom2->Classify(gp_Pnt2d(u2,v2),TolX) != TopAbs_OUT);
 384         }
 385         if (ok) {
 386           sline.Append(newc);
 387           slineS1.Append(H1);
 388           slineS2.Append(H1);
 389         }
 390       }
 391     }// end of for (i=1; i<=myLConstruct.NbParts(); i++)
 392    }// case IntPatch_Lin:  case IntPatch_Parabola:  case IntPatch_Hyperbola:
 393    break;
 394
 395   //########################################
 396   // Circle and Ellipse
 397   //########################################
 398   case IntPatch_Circle:
 399   case IntPatch_Ellipse: {
 400
 401     if (typl == IntPatch_Circle) {
 402       newc = new Geom_Circle
 403         (Handle(IntPatch_GLine)::DownCast(L)->Circle());
 404     }
 405     else {
 406       newc = new Geom_Ellipse
 407         (Handle(IntPatch_GLine)::DownCast(L)->Ellipse());
 408     }
 409     //
 410     Standard_Real aPeriod, aRealEpsilon;
 411     //
 412     aRealEpsilon=RealEpsilon();
 413     aPeriod=M_PI+M_PI;
 414     //
 415     aNbParts=myLConstruct.NbParts();
 416     //
 417     for (i=1; i<=aNbParts; i++) {
 418       myLConstruct.Part(i, fprm, lprm);
 419       //
 420       if (Abs(fprm) > aRealEpsilon || Abs(lprm-aPeriod) > aRealEpsilon) {
 421         //==============================================
 422         Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
 423         //
 424         sline.Append(aTC3D);
 425         //
 426         fprm=aTC3D->FirstParameter();
 427         lprm=aTC3D->LastParameter ();
 428         ////
 429         if(myApprox1) {
 430           Handle (Geom2d_Curve) C2d;
 431           BuildPCurves(fprm,lprm,Tolpc,myHS1->ChangeSurface().Surface(),newc,C2d);
 432           if(Tolpc>myTolReached2d || myTolReached2d==0.) {
 433             myTolReached2d=Tolpc;
 434           }
 435           slineS1.Append(C2d);
 436         }
 437         else { ////
 438           slineS1.Append(H1);
 439         }
 440         //
 441         if(myApprox2) {
 442           Handle (Geom2d_Curve) C2d;
 443           BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
 444           if(Tolpc>myTolReached2d || myTolReached2d==0) {
 445             myTolReached2d=Tolpc;
 446           }
 447           slineS2.Append(C2d);
 448         }
 449         else {
 450           slineS2.Append(H1);
 451         }
 452         //==============================================
 453       } //if (Abs(fprm) > RealEpsilon() || Abs(lprm-2.*M_PI) > RealEpsilon())
 454       //
 455       else {//  on regarde si on garde
 456         //
 457         if (aNbParts==1) {
 458           if (Abs(fprm) < RealEpsilon() &&  Abs(lprm-2.*M_PI) < RealEpsilon()) {
 459             Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
 460             //
 461             sline.Append(aTC3D);
 462             fprm=aTC3D->FirstParameter();
 463             lprm=aTC3D->LastParameter ();
 464
 465             if(myApprox1) {
 466               Handle (Geom2d_Curve) C2d;
 467               BuildPCurves(fprm,lprm,Tolpc,myHS1->ChangeSurface().Surface(),newc,C2d);
 468               if(Tolpc>myTolReached2d || myTolReached2d==0) {
 469                 myTolReached2d=Tolpc;
 470               }
 471               slineS1.Append(C2d);
 472             }
 473             else { ////
 474               slineS1.Append(H1);
 475             }
 476
 477             if(myApprox2) {
 478               Handle (Geom2d_Curve) C2d;
 479               BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
 480               if(Tolpc>myTolReached2d || myTolReached2d==0) {
 481                 myTolReached2d=Tolpc;
 482               }
 483               slineS2.Append(C2d);
 484             }
 485             else {
 486               slineS2.Append(H1);
 487             }
 488             break;
 489           }
 490         }
 491         //
 492         Standard_Real aTwoPIdiv17, u1, v1, u2, v2, TolX;
 493         //
 494         aTwoPIdiv17=2.*M_PI/17.;
 495         //
 496         for (j=0; j<=17; j++) {
 497           gp_Pnt ptref (newc->Value (j*aTwoPIdiv17));
 498           TolX = Precision::Confusion();
 499
 500           Parameters(myHS1, myHS2, ptref, u1, v1, u2, v2);
 501           ok = (dom1->Classify(gp_Pnt2d(u1,v1),TolX) != TopAbs_OUT);
 502           if(ok) {
 503             ok = (dom2->Classify(gp_Pnt2d(u2,v2),TolX) != TopAbs_OUT);
 504           }
 505           if (ok) {
 506             sline.Append(newc);
 507             //==============================================
 508             if(myApprox1) {
 509               Handle (Geom2d_Curve) C2d;
 510               BuildPCurves(fprm, lprm, Tolpc, myHS1->ChangeSurface().Surface(), newc, C2d);
 511               if(Tolpc>myTolReached2d || myTolReached2d==0) {
 512                 myTolReached2d=Tolpc;
 513               }
 514               slineS1.Append(C2d);
 515             }
 516             else {
 517               slineS1.Append(H1);
 518             }
 519
 520             if(myApprox2) {
 521               Handle (Geom2d_Curve) C2d;
 522               BuildPCurves(fprm, lprm, Tolpc,myHS2->ChangeSurface().Surface(), newc, C2d);
 523               if(Tolpc>myTolReached2d || myTolReached2d==0) {
 524                 myTolReached2d=Tolpc;
 525               }
 526               slineS2.Append(C2d);
 527             }
 528             else {
 529               slineS2.Append(H1);
 530             }
 531             break;
 532           }//  end of if (ok) {
 533         }//  end of for (Standard_Integer j=0; j<=17; j++)
 534       }//  end of else { on regarde si on garde
 535     }// for (i=1; i<=myLConstruct.NbParts(); i++)
 536   }// IntPatch_Circle: IntPatch_Ellipse:
 537     break;
 538
 539     //########################################
 540     // Analytic
 541     //########################################
 542   case IntPatch_Analytic: {
 543     IntSurf_Quadric quad1,quad2;
 544     //
 545     GetQuadric(myHS1, quad1);
 546     GetQuadric(myHS2, quad2);
 547     //=========
 548     IntPatch_ALineToWLine convert (quad1, quad2);
 549
 550     if (!myApprox) {
 551       Handle(Geom2d_BSplineCurve) aH1, aH2;
 552       //
 553       aNbParts=myLConstruct.NbParts();
 554       for (i=1; i<=aNbParts; i++) {
 555         myLConstruct.Part(i, fprm, lprm);
 556         Handle(IntPatch_WLine) WL =
 557           convert.MakeWLine(Handle(IntPatch_ALine)::DownCast(L), fprm, lprm);
 558         //
 559         if(myApprox1) {
 560           aH1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
 561         }
 562
 563         if(myApprox2) {
 564           aH2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
 565         }
 566         sline.Append(MakeBSpline(WL,1,WL->NbPnts()));
 567         slineS1.Append(aH1);
 568         slineS2.Append(aH2);
 569       }
 570     } // if (!myApprox)
 571
 572     else { // myApprox=TRUE
 573       GeomInt_WLApprox theapp3d;
 574       Standard_Real tol2d = myTolApprox;
 575       //
 576       theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
 577
 578       aNbParts=myLConstruct.NbParts();
 579       for (i=1; i<=aNbParts; i++) {
 580         myLConstruct.Part(i, fprm, lprm);
 581         Handle(IntPatch_WLine) WL =
 582           convert.MakeWLine(Handle(IntPatch_ALine):: DownCast(L),fprm,lprm);
 583
 584         theapp3d.Perform(myHS1,myHS2,WL,Standard_True,myApprox1,myApprox2, 1, WL->NbPnts());
 585         if (!theapp3d.IsDone()) {
 586           //
 587           Handle(Geom2d_BSplineCurve) aH1, aH2;
 588
 589           if(myApprox1) {
 590             aH1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
 591           }
 592
 593           if(myApprox2) {
 594             aH2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
 595           }
 596           sline.Append(MakeBSpline(WL,1,WL->NbPnts()));
 597           slineS1.Append(aH1);
 598           slineS2.Append(aH1);
 599         }
 600         //
 601         else {
 602           if(myApprox1 || myApprox2) {
 603             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0) {
 604               myTolReached2d = theapp3d.TolReached2d();
 605             }
 606           }
 607
 608           if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0) {
 609             myTolReached3d = theapp3d.TolReached3d();
 610           }
 611
 612           Standard_Integer aNbMultiCurves, nbpoles;
 613           aNbMultiCurves=theapp3d.NbMultiCurves();
 614           for (j=1; j<=aNbMultiCurves; j++) {
 615             const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
 616
 617             nbpoles = mbspc.NbPoles();
 618             TColgp_Array1OfPnt tpoles(1, nbpoles);
 619             mbspc.Curve(1, tpoles);
 620             Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
 621                                                                mbspc.Knots(),
 622                                                                mbspc.Multiplicities(),
 623                                                                mbspc.Degree());
 624
 625             GeomLib_CheckBSplineCurve Check(BS, myTolCheck, myTolAngCheck);
 626             Check.FixTangent(Standard_True,Standard_True);
 627             //
 628             sline.Append(BS);
 629             //
 630             if(myApprox1) {
 631               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
 632               mbspc.Curve(2,tpoles2d);
 633               Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
 634                                                                       mbspc.Knots(),
 635                                                                       mbspc.Multiplicities(),
 636                                                                       mbspc.Degree());
 637
 638               GeomLib_Check2dBSplineCurve newCheck(BS2,myTolCheck,myTolAngCheck);
 639               newCheck.FixTangent(Standard_True,Standard_True);
 640               slineS1.Append(BS2);
 641             }
 642             else {
 643               slineS1.Append(H1);
 644             }
 645
 646             if(myApprox2) {
 647               TColgp_Array1OfPnt2d tpoles2d(1, nbpoles);
 648               Standard_Integer TwoOrThree;
 649               TwoOrThree=myApprox1 ? 3 : 2;
 650               mbspc.Curve(TwoOrThree, tpoles2d);
 651               Handle(Geom2d_BSplineCurve) BS2 =new Geom2d_BSplineCurve(tpoles2d,
 652                                                                        mbspc.Knots(),
 653                                                                        mbspc.Multiplicities(),
 654                                                                        mbspc.Degree());
 655
 656               GeomLib_Check2dBSplineCurve newCheck(BS2,myTolCheck,myTolAngCheck);
 657               newCheck.FixTangent(Standard_True,Standard_True);
 658               //
 659               slineS2.Append(BS2);
 660             }
 661             else {
 662               slineS2.Append(H1);
 663             }
 664             //
 665           }// for (j=1; j<=aNbMultiCurves; j++) {
 666         }// else from if (!theapp3d.IsDone())
 667       }// for (i=1; i<=aNbParts; i++) {
 668     }// else { // myApprox=TRUE
 669   }// case IntPatch_Analytic:
 670     break;
 671
 672     //########################################
 673     // Walking
 674     //########################################
 675   case IntPatch_Walking:{
 676     Handle(IntPatch_WLine) WL =
 677       Handle(IntPatch_WLine)::DownCast(L);
 678     //
 679     Standard_Integer ifprm, ilprm;
 680     //
 681     if (!myApprox) {
 682       aNbParts=myLConstruct.NbParts();
 683       for (i=1; i<=aNbParts; i++) {
 684         myLConstruct.Part(i, fprm, lprm);
 685         ifprm=(Standard_Integer)fprm;
 686         ilprm=(Standard_Integer)lprm;
 687         //
 688         Handle(Geom2d_BSplineCurve) aH1, aH2;
 689
 690         if(myApprox1) {
 691           aH1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
 692         }
 693         if(myApprox2) {
 694           aH2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
 695         }
 696         //
 697         Handle(Geom_Curve) aBSp=MakeBSpline(WL, ifprm, ilprm);
 698         //
 699         sline.Append(aBSp);
 700         slineS1.Append(aH1);
 701         slineS2.Append(aH2);
 702       }
 703     }
 704     //
 705     else {
 706       Standard_Boolean bIsDecomposited;
 707       Standard_Integer nbiter, aNbSeqOfL;
 708       GeomInt_WLApprox theapp3d;
 709       IntPatch_SequenceOfLine aSeqOfL;
 710       Standard_Real tol2d, aTolSS;
 711       //
 712       tol2d = myTolApprox;
 713       aTolSS=2.e-7;
 714       if(myHS1 == myHS2) {
 715         theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False);
 716       }
 717       else {
 718         theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
 719       }
 720       //
 721       bIsDecomposited =
 722         DecompositionOfWLine(WL, myHS1, myHS2, aTolSS, myLConstruct, aSeqOfL);
 723       //
 724       aNbParts=myLConstruct.NbParts();
 725       aNbSeqOfL=aSeqOfL.Length();
 726       //
 727       nbiter = (bIsDecomposited) ? aNbSeqOfL : aNbParts;
 728       //
 729       for(i = 1; i <= nbiter; i++) {
 730         if(bIsDecomposited) {
 731           WL = Handle(IntPatch_WLine)::DownCast(aSeqOfL.Value(i));
 732           ifprm = 1;
 733           ilprm = WL->NbPnts();
 734         }
 735         else {
 736           myLConstruct.Part(i, fprm, lprm);
 737           ifprm = (Standard_Integer)fprm;
 738           ilprm = (Standard_Integer)lprm;
 739         }
 740         //-- lbr :
 741         //-- Si une des surfaces est un plan , on approxime en 2d
 742         //-- sur cette surface et on remonte les points 2d en 3d.
 743         GeomAbs_SurfaceType typs1, typs2;
 744         typs1 = myHS1->Surface().GetType();
 745         typs2 = myHS2->Surface().GetType();
 746         //
 747         if(typs1 == GeomAbs_Plane) {
 748           theapp3d.Perform(myHS1, myHS2, WL, Standard_False,
 749                            Standard_True, myApprox2,
 750                            ifprm,  ilprm);
 751         }
 752         else if(typs2 == GeomAbs_Plane) {
 753           theapp3d.Perform(myHS1,myHS2,WL,Standard_False,
 754                            myApprox1,Standard_True,
 755                            ifprm,  ilprm);
 756         }
 757         else {
 758           //
 759           if (myHS1 != myHS2){
 760             if ((typs1==GeomAbs_BezierSurface || typs1==GeomAbs_BSplineSurface) &&
 761                 (typs2==GeomAbs_BezierSurface || typs2==GeomAbs_BSplineSurface)) {
 762
 763               theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
 764               //Standard_Boolean bUseSurfaces;
 765               //bUseSurfaces=NotUseSurfacesForApprox(myFace1, myFace2, WL, ifprm,  ilprm);
 766               //if (bUseSurfaces) {
 767                 //theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False);
 768               //}
 769             }
 770           }
 771           //
 772           theapp3d.Perform(myHS1,myHS2,WL,Standard_True,
 773                            myApprox1,myApprox2,
 774                            ifprm,  ilprm);
 775         }
 776
 777         if (!theapp3d.IsDone()) {
 778           //
 779           Handle(Geom2d_BSplineCurve) aH1, aH2;
 780           //
 781           Handle(Geom_Curve) aBSp=MakeBSpline(WL, ifprm, ilprm);
 782           if(myApprox1) {
 783             aH1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
 784           }
 785           if(myApprox2) {
 786             aH2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
 787           }
 788           //
 789           sline.Append(aBSp);
 790           slineS1.Append(aH1);
 791           slineS2.Append(aH2);
 792         }//if (!theapp3d.IsDone())
 793
 794         else {
 795           if(myApprox1 || myApprox2 || (typs1==GeomAbs_Plane || typs2==GeomAbs_Plane)) {
 796             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0.) {
 797               myTolReached2d = theapp3d.TolReached2d();
 798             }
 799           }
 800           if(typs1==GeomAbs_Plane || typs2==GeomAbs_Plane) {
 801             myTolReached3d = myTolReached2d;
 802           }
 803           else  if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0.) {
 804             myTolReached3d = theapp3d.TolReached3d();
 805           }
 806
 807           Standard_Integer aNbMultiCurves, nbpoles;
 808           //
 809           aNbMultiCurves=theapp3d.NbMultiCurves();
 810           for (j=1; j<=aNbMultiCurves; j++) {
 811             if(typs1 == GeomAbs_Plane) {
 812               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
 813               nbpoles = mbspc.NbPoles();
 814
 815               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
 816               TColgp_Array1OfPnt   tpoles(1,nbpoles);
 817
 818               mbspc.Curve(1,tpoles2d);
 819               const gp_Pln&  Pln = myHS1->Surface().Plane();
 820               //
 821               Standard_Integer ik;
 822               for(ik = 1; ik<= nbpoles; ik++) {
 823                 tpoles.SetValue(ik,
 824                                 ElSLib::Value(tpoles2d.Value(ik).X(),
 825                                               tpoles2d.Value(ik).Y(),
 826                                               Pln));
 827               }
 828               //
 829               Handle(Geom_BSplineCurve) BS =
 830                 new Geom_BSplineCurve(tpoles,
 831                                       mbspc.Knots(),
 832                                       mbspc.Multiplicities(),
 833                                       mbspc.Degree());
 834               GeomLib_CheckBSplineCurve Check(BS,myTolCheck,myTolAngCheck);
 835               Check.FixTangent(Standard_True, Standard_True);
 836               //
 837               sline.Append(BS);
 838               //
 839               if(myApprox1) {
 840                 Handle(Geom2d_BSplineCurve) BS1 =
 841                   new Geom2d_BSplineCurve(tpoles2d,
 842                                           mbspc.Knots(),
 843                                           mbspc.Multiplicities(),
 844                                           mbspc.Degree());
 845                 GeomLib_Check2dBSplineCurve Check1(BS1,myTolCheck,myTolAngCheck);
 846                 Check1.FixTangent(Standard_True,Standard_True);
 847                 //
 848                 AdjustUPeriodic (aS1, BS1);
 849                 //
 850                 slineS1.Append(BS1);
 851               }
 852               else {
 853                 slineS1.Append(H1);
 854               }
 855
 856               if(myApprox2) {
 857                 mbspc.Curve(2, tpoles2d);
 858
 859                 Handle(Geom2d_BSplineCurve) BS2 = new Geom2d_BSplineCurve(tpoles2d,
 860                                                                           mbspc.Knots(),
 861                                                                           mbspc.Multiplicities(),
 862                                                                           mbspc.Degree());
 863                 GeomLib_Check2dBSplineCurve newCheck(BS2,myTolCheck,myTolAngCheck);
 864                 newCheck.FixTangent(Standard_True,Standard_True);
 865                 //
 866                 AdjustUPeriodic (aS2, BS2);
 867                 //
 868                 slineS2.Append(BS2);
 869               }
 870               else {
 871                 slineS2.Append(H1);
 872               }
 873             }//if(typs1 == GeomAbs_Plane)
 874             //
 875             else if(typs2 == GeomAbs_Plane) {
 876               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
 877               nbpoles = mbspc.NbPoles();
 878
 879               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
 880               TColgp_Array1OfPnt   tpoles(1,nbpoles);
 881               mbspc.Curve((myApprox1==Standard_True)? 2 : 1,tpoles2d);
 882               const gp_Pln&  Pln = myHS2->Surface().Plane();
 883               //
 884               Standard_Integer ik;
 885               for(ik = 1; ik<= nbpoles; ik++) {
 886                 tpoles.SetValue(ik,
 887                                 ElSLib::Value(tpoles2d.Value(ik).X(),
 888                                               tpoles2d.Value(ik).Y(),
 889                                               Pln));
 890
 891               }
 892               //
 893               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
 894                                                                  mbspc.Knots(),
 895                                                                  mbspc.Multiplicities(),
 896                                                                  mbspc.Degree());
 897               GeomLib_CheckBSplineCurve Check(BS,myTolCheck,myTolAngCheck);
 898               Check.FixTangent(Standard_True,Standard_True);
 899               //
 900               sline.Append(BS);
 901               //
 902               if(myApprox2) {
 903                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
 904                                                                         mbspc.Knots(),
 905                                                                         mbspc.Multiplicities(),
 906                                                                         mbspc.Degree());
 907                 GeomLib_Check2dBSplineCurve Check1(BS1,myTolCheck,myTolAngCheck);
 908                 Check1.FixTangent(Standard_True,Standard_True);
 909                 //
 910                 //
 911                 AdjustUPeriodic (aS2, BS1);
 912                 //
 913                 slineS2.Append(BS1);
 914               }
 915               else {
 916                 slineS2.Append(H1);
 917               }
 918
 919               if(myApprox1) {
 920                 mbspc.Curve(1,tpoles2d);
 921                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
 922                                                                         mbspc.Knots(),
 923                                                                         mbspc.Multiplicities(),
 924                                                                         mbspc.Degree());
 925                 GeomLib_Check2dBSplineCurve Check2(BS2,myTolCheck,myTolAngCheck);
 926                 Check2.FixTangent(Standard_True,Standard_True);
 927                 //
 928                 //
 929                 AdjustUPeriodic (aS1, BS2);
 930                 //
 931                 slineS1.Append(BS2);
 932               }
 933               else {
 934                 slineS1.Append(H1);
 935               }
 936             } // else if(typs2 == GeomAbs_Plane)
 937             //
 938             else { // typs1!=GeomAbs_Plane && typs2!=GeomAbs_Plane
 939               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
 940               nbpoles = mbspc.NbPoles();
 941               TColgp_Array1OfPnt tpoles(1,nbpoles);
 942               mbspc.Curve(1,tpoles);
 943               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
 944                                                                  mbspc.Knots(),
 945                                                                  mbspc.Multiplicities(),
 946                                                                  mbspc.Degree());
 947               GeomLib_CheckBSplineCurve Check(BS,myTolCheck,myTolAngCheck);
 948               Check.FixTangent(Standard_True,Standard_True);
 949               //
 950         //Check IsClosed()
 951         Standard_Real aDist = Max(BS->StartPoint().XYZ().SquareModulus(),
 952                                   BS->EndPoint().XYZ().SquareModulus());
 953         Standard_Real eps = Epsilon(aDist);
 954         if(BS->StartPoint().SquareDistance(BS->EndPoint()) < 2.*eps &&
 955           !BS->IsClosed() && !BS->IsPeriodic())
 956         {
 957           //force Closed()
 958           gp_Pnt aPm((BS->Pole(1).XYZ() + BS->Pole(BS->NbPoles()).XYZ()) / 2.);
 959           BS->SetPole(1, aPm);
 960           BS->SetPole(BS->NbPoles(), aPm);
 961         }
 962               sline.Append(BS);
 963
 964               if(myApprox1) {
 965                 TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
 966                 mbspc.Curve(2,tpoles2d);
 967                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
 968                                                                         mbspc.Knots(),
 969                                                                         mbspc.Multiplicities(),
 970                                                                         mbspc.Degree());
 971                 GeomLib_Check2dBSplineCurve newCheck(BS1,myTolCheck,myTolAngCheck);
 972                 newCheck.FixTangent(Standard_True,Standard_True);
 973                 //
 974                 AdjustUPeriodic (aS1, BS1);
 975                 //
 976                 slineS1.Append(BS1);
 977               }
 978               else {
 979                 slineS1.Append(H1);
 980               }
 981               if(myApprox2) {
 982                 TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
 983                 mbspc.Curve((myApprox1==Standard_True)? 3 : 2,tpoles2d);
 984                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
 985                                                                         mbspc.Knots(),
 986                                                                         mbspc.Multiplicities(),
 987                                                                         mbspc.Degree());
 988                 GeomLib_Check2dBSplineCurve newCheck(BS2,myTolCheck,myTolAngCheck);
 989                 newCheck.FixTangent(Standard_True,Standard_True);
 990                 //
 991                 AdjustUPeriodic (aS2, BS2);
 992                 //
 993                 slineS2.Append(BS2);
 994               }
 995               else {
 996                 slineS2.Append(H1);
 997               }
 998             }// else { // typs1!=GeomAbs_Plane && typs2!=GeomAbs_Plane
 999           }// for (j=1; j<=aNbMultiCurves; j++
1000         }
1001       }
1002     }
1003   }
1004     break;
1005
1006   case IntPatch_Restriction:
1007     break;
1008
1009   }
1010 }
1011 //
1012 //=======================================================================
1013 //function : AdjustUPeriodic
1014 //purpose  :
1015 //=======================================================================
1016  void AdjustUPeriodic (const Handle(Geom_Surface)& aS, Handle(Geom2d_Curve)& aC2D)
1017 {
1018   if (aC2D.IsNull() || !aS->IsUPeriodic())
1019     return;
1020   //
1021   const Standard_Real aEps=Precision::PConfusion();//1.e-9
1022   const Standard_Real aEpsilon=Epsilon(10.);//1.77e-15
1023   //
1024   Standard_Real umin,umax,vmin,vmax;
1025   aS->Bounds(umin,umax,vmin,vmax);
1026   const Standard_Real aPeriod = aS->UPeriod();
1027
1028   const Standard_Real aT1=aC2D->FirstParameter();
1029   const Standard_Real aT2=aC2D->LastParameter();
1030   const Standard_Real aTx=aT1+0.467*(aT2-aT1);
1031   const gp_Pnt2d aPx=aC2D->Value(aTx);
1032   //
1033   Standard_Real aUx=aPx.X();
1034   if (fabs(aUx)<aEpsilon)
1035     aUx=0.;
1036   if (fabs(aUx-aPeriod)<aEpsilon)
1037     aUx=aPeriod;
1038   //
1039   Standard_Real dU=0.;
1040   while(aUx <(umin-aEps)) {
1041     aUx+=aPeriod;
1042     dU+=aPeriod;
1043   }
1044   while(aUx>(umax+aEps)) {
1045     aUx-=aPeriod;
1046     dU-=aPeriod;
1047   }
1048   //
1049   if (dU!=0.) {
1050     gp_Vec2d aV2D(dU, 0.);
1051     aC2D->Translate(aV2D);
1052   }
1053 }
1054 //
1055 //
1056 //=======================================================================
1057 //function : BuildPCurves
1058 //purpose  :
1059 //=======================================================================
1060  void BuildPCurves (Standard_Real f, Standard_Real l,
1061                     Standard_Real& Tol,
1062                     const Handle(Geom_Surface)& S,
1063                     const Handle(Geom_Curve)& C,
1064                     Handle(Geom2d_Curve)& C2d)
1065 {
1066   if (C2d.IsNull())
1067   {
1068     C2d = GeomProjLib::Curve2d(C,f,l,S,Tol);
1069     if (S->IsUPeriodic())
1070     {
1071       Standard_Real umin,umax,vmin,vmax;
1072       S->Bounds(umin,umax,vmin,vmax);
1073       // Recadre dans le domaine UV de la face
1074       const Standard_Real period = S->UPeriod();
1075       gp_Pnt2d Pf=C2d->Value(f);
1076       //
1077       Standard_Real U0=Pf.X();
1078       //
1079       const Standard_Real aEpsilon=Epsilon(10.);//1.77e-15
1080       if (fabs(U0)<aEpsilon)
1081         U0=0.;
1082       if (fabs(U0-period)<aEpsilon)
1083         U0=period;
1084       //
1085       const Standard_Real aEps=Precision::PConfusion();//1.e-9
1086       Standard_Real du=0.;
1087       while(U0 <(umin-aEps)) {
1088         U0+=period;
1089         du+=period;
1090       }
1091       while(U0>(umax+aEps)) {
1092         U0-=period;
1093         du-=period;
1094       }
1095       if (du!=0.) {
1096         gp_Vec2d T1(du, 0.);
1097         C2d->Translate(T1);
1098       }
1099     }
1100   }
1101 }
1102 //
1103 //=======================================================================
1104 //function : GetQuadric
1105 //purpose  :
1106 //=======================================================================
1107  void GetQuadric(const Handle(GeomAdaptor_HSurface)& HS1, IntSurf_Quadric& quad1)
1108 {
1109   switch (HS1->Surface().GetType())
1110   {
1111     case GeomAbs_Plane:    quad1.SetValue(HS1->Surface().Plane()); break;
1112     case GeomAbs_Cylinder: quad1.SetValue(HS1->Surface().Cylinder()); break;
1113     case GeomAbs_Cone:     quad1.SetValue(HS1->Surface().Cone()); break;
1114     case GeomAbs_Sphere:   quad1.SetValue(HS1->Surface().Sphere()); break;
1115     case GeomAbs_Torus:    quad1.SetValue(HS1->Surface().Torus()); break;
1116     default: Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
1117   }
1118 }
1119
1120 //=======================================================================
1121 //function : Parameters
1122 //purpose  :
1123 //=======================================================================
1124  void Parameters(const Handle(GeomAdaptor_HSurface)& HS1,
1125                  const Handle(GeomAdaptor_HSurface)& HS2,
1126                  const gp_Pnt& Ptref,
1127                  Standard_Real& U1,
1128                  Standard_Real& V1,
1129                  Standard_Real& U2,
1130                  Standard_Real& V2)
1131 {
1132   IntSurf_Quadric quad1,quad2;
1133   //
1134   GetQuadric(HS1, quad1);
1135   GetQuadric(HS2, quad2);
1136   //
1137   quad1.Parameters(Ptref,U1,V1);
1138   quad2.Parameters(Ptref,U2,V2);
1139 }
1140
1141 //=======================================================================
1142 //function : MakeBSpline
1143 //purpose  :
1144 //=======================================================================
1145 Handle(Geom_Curve) MakeBSpline  (const Handle(IntPatch_WLine)& WL,
1146                                  const Standard_Integer ideb,
1147                                  const Standard_Integer ifin)
1148 {
1149   const Standard_Integer nbpnt = ifin-ideb+1;
1150   TColgp_Array1OfPnt poles(1,nbpnt);
1151   TColStd_Array1OfReal knots(1,nbpnt);
1152   TColStd_Array1OfInteger mults(1,nbpnt);
1153   Standard_Integer i = 1, ipidebm1 = ideb;
1154   for(; i<=nbpnt; ipidebm1++, i++)
1155   {
1156     poles(i) = WL->Point(ipidebm1).Value();
1157     mults(i) = 1;
1158     knots(i) = i-1;
1159   }
1160   mults(1) = mults(nbpnt) = 2;
1161   return new Geom_BSplineCurve(poles,knots,mults,1);
1162 }
1163
1164 //=======================================================================
1165 //function : MakeBSpline2d
1166 //purpose  :
1167 //=======================================================================
1168 Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
1169                                           const Standard_Integer ideb,
1170                                           const Standard_Integer ifin,
1171                                           const Standard_Boolean onFirst)
1172 {
1173   const Standard_Integer nbpnt = ifin-ideb+1;
1174   TColgp_Array1OfPnt2d poles(1,nbpnt);
1175   TColStd_Array1OfReal knots(1,nbpnt);
1176   TColStd_Array1OfInteger mults(1,nbpnt);
1177   Standard_Integer i = 1, ipidebm1 = ideb;
1178   for(; i <= nbpnt; ipidebm1++, i++)
1179   {
1180     Standard_Real U, V;
1181     if(onFirst)
1182           theWLine->Point(ipidebm1).ParametersOnS1(U, V);
1183     else
1184           theWLine->Point(ipidebm1).ParametersOnS2(U, V);
1185     poles(i).SetCoord(U, V);
1186     mults(i) = 1;
1187     knots(i) = i-1;
1188   }
1189   mults(1) = mults(nbpnt) = 2;
1190   return new Geom2d_BSplineCurve(poles,knots,mults,1);
1191 }
1192
1193 //=========================================================================
1194 // static function : ComputePurgedWLine
1195 // purpose : Removes equal points (leave one of equal points) from theWLine
1196 //           and recompute vertex parameters.
1197 //           Returns new WLine or null WLine if the number
1198 //           of the points is less than 2.
1199 //=========================================================================
1200 Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine)
1201 {
1202   Handle(IntPatch_WLine) aResult;
1203   Handle(IntPatch_WLine) aLocalWLine;
1204   Handle(IntPatch_WLine) aTmpWLine = theWLine;
1205
1206   Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
1207   aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
1208   Standard_Integer i, k, v, nb, nbvtx;
1209   nbvtx = theWLine->NbVertex();
1210   nb = theWLine->NbPnts();
1211
1212   for(i = 1; i <= nb; i++) {
1213     aLineOn2S->Add(theWLine->Point(i));
1214   }
1215
1216   for(v = 1; v <= nbvtx; v++) {
1217     aLocalWLine->AddVertex(theWLine->Vertex(v));
1218   }
1219
1220   for(i = 1; i <= aLineOn2S->NbPoints(); i++) {
1221     Standard_Integer aStartIndex = i + 1;
1222     Standard_Integer anEndIndex = i + 5;
1223     nb = aLineOn2S->NbPoints();
1224     anEndIndex = (anEndIndex > nb) ? nb : anEndIndex;
1225
1226     if((aStartIndex >= nb) || (anEndIndex <= 1)) {
1227       continue;
1228     }
1229     k = aStartIndex;
1230
1231     while(k <= anEndIndex) {
1232
1233       if(i != k) {
1234         IntSurf_PntOn2S p1 = aLineOn2S->Value(i);
1235         IntSurf_PntOn2S p2 = aLineOn2S->Value(k);
1236
1237         if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) {
1238           aTmpWLine = aLocalWLine;
1239           aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
1240
1241           for(v = 1; v <= aTmpWLine->NbVertex(); v++) {
1242             IntPatch_Point aVertex = aTmpWLine->Vertex(v);
1243             Standard_Integer avertexindex = (Standard_Integer)aVertex.ParameterOnLine();
1244
1245             if(avertexindex >= k) {
1246               aVertex.SetParameter(aVertex.ParameterOnLine() - 1.);
1247             }
1248             aLocalWLine->AddVertex(aVertex);
1249           }
1250           aLineOn2S->RemovePoint(k);
1251           anEndIndex--;
1252           continue;
1253         }
1254       }
1255       k++;
1256     }
1257   }
1258
1259   if(aLineOn2S->NbPoints() > 1) {
1260     aResult = aLocalWLine;
1261   }
1262   return aResult;
1263 }
1264 //=======================================================================
1265 //function : DecompositionOfWLine
1266 //purpose  :
1267 //=======================================================================
1268 Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
1269                                       const Handle(GeomAdaptor_HSurface)&            theSurface1,
1270                                       const Handle(GeomAdaptor_HSurface)&            theSurface2,
1271                                       const Standard_Real aTolSum,
1272                                       const GeomInt_LineConstructor&                 theLConstructor,
1273                                       IntPatch_SequenceOfLine&                       theNewLines)
1274 {
1275   typedef NCollection_List<Standard_Integer> ListOfInteger;
1276   //have to use std::vector, not NCollection_Vector in order to use copy constructor of
1277   //ListOfInteger which will be created with specific allocator instance
1278   typedef std::vector<ListOfInteger, NCollection_StdAllocator<
1279       ListOfInteger> > ArrayOfListOfInteger;
1280
1281   Standard_Boolean bIsPrevPointOnBoundary, bIsCurrentPointOnBoundary;
1282   Standard_Integer nblines, aNbPnts, aNbParts, pit, i, j, aNbListOfPointIndex;
1283   Standard_Real aTol, umin, umax, vmin, vmax;
1284
1285   //an inc allocator, it will contain wasted space (upon list's Clear()) but it should
1286   //still be faster than the standard allocator, and wasted memory should not be
1287   //significant and will be limited by time span of this function;
1288   //this is a separate allocator from the anIncAlloc below what provides better data
1289   //locality in the latter (by avoiding wastes which will only be in anIdxAlloc)
1290   Handle(NCollection_IncAllocator) anIdxAlloc = new NCollection_IncAllocator();
1291   ListOfInteger aListOfPointIndex (anIdxAlloc);
1292
1293   //GeomAPI_ProjectPointOnSurf aPrj1, aPrj2;
1294   ProjectPointOnSurf aPrj1, aPrj2;
1295   Handle(Geom_Surface) aSurf1,  aSurf2;
1296   //
1297   aNbParts=theLConstructor.NbParts();
1298   aNbPnts=theWLine->NbPnts();
1299   //
1300   if((!aNbPnts) || (!aNbParts)){
1301     return Standard_False;
1302   }
1303   //
1304   Handle(NCollection_IncAllocator) anIncAlloc = new NCollection_IncAllocator();
1305   NCollection_StdAllocator<ListOfInteger> anAlloc (anIncAlloc);
1306   const ListOfInteger aDummy (anIncAlloc); //empty list to be copy constructed from
1307   ArrayOfListOfInteger anArrayOfLines (aNbPnts + 1, aDummy, anAlloc);
1308
1309   NCollection_LocalArray<Standard_Integer> anArrayOfLineTypeArr (aNbPnts + 1);
1310   Standard_Integer* anArrayOfLineType = anArrayOfLineTypeArr;
1311   //
1312   nblines = 0;
1313   aTol = Precision::Confusion();
1314   //
1315   aSurf1 = theSurface1->ChangeSurface().Surface();
1316   aSurf1->Bounds(umin, umax, vmin, vmax);
1317   aPrj1.Init(aSurf1, umin, umax, vmin, vmax);
1318   //
1319   aSurf2 = theSurface2->ChangeSurface().Surface();
1320   aSurf2->Bounds(umin, umax, vmin, vmax);
1321   aPrj2.Init(aSurf2, umin, umax, vmin, vmax);
1322   //
1323   //
1324   bIsPrevPointOnBoundary=Standard_False;
1325   for(pit=1; pit<=aNbPnts; pit++) {
1326     const IntSurf_PntOn2S& aPoint = theWLine->Point(pit);
1327     bIsCurrentPointOnBoundary=Standard_False;
1328     //
1329     // whether aPoint is on boundary or not
1330     //
1331     for(i=0; i<2; i++) {// exploration Surface 1,2
1332       Handle(GeomAdaptor_HSurface) aGASurface = (!i) ? theSurface1 : theSurface2;
1333       aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
1334       //
1335       for(j=0; j<2; j++) {// exploration of coordinate U,V
1336         Standard_Boolean isperiodic;
1337         //
1338         isperiodic = (!j) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
1339         if(!isperiodic) {
1340           continue;
1341         }
1342         //
1343         Standard_Real aResolution, aPeriod, alowerboundary, aupperboundary, U, V;
1344         Standard_Real aParameter, anoffset, anAdjustPar;
1345         Standard_Boolean bIsOnFirstBoundary, bIsPointOnBoundary;
1346         //
1347         aResolution = (!j) ? aGASurface->UResolution(aTol) : aGASurface->VResolution(aTol);
1348         aPeriod     = (!j) ? aGASurface->UPeriod() : aGASurface->VPeriod();
1349         alowerboundary = (!j) ? umin : vmin;
1350         aupperboundary = (!j) ? umax : vmax;
1351         U=0.;V=0.;//?
1352         //
1353         if(!i){
1354           aPoint.ParametersOnS1(U, V);
1355         }
1356         else{
1357           aPoint.ParametersOnS2(U, V);
1358         }
1359         //
1360         aParameter = (!j) ? U : V;
1361         anoffset=0.;
1362         anAdjustPar=AdjustPeriodic(aParameter, alowerboundary, aupperboundary, aPeriod, anoffset);
1363         //
1364         bIsOnFirstBoundary=Standard_True;
1365         //
1366         bIsPointOnBoundary=
1367           IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
1368
1369         if(bIsPointOnBoundary) {
1370           bIsCurrentPointOnBoundary = Standard_True;
1371           break;
1372         }
1373       }// for(j=0; j<2; j++)
1374
1375       if(bIsCurrentPointOnBoundary){
1376         break;
1377       }
1378     }// for(i=0; i<2; i++)
1379     //
1380     if((bIsCurrentPointOnBoundary != bIsPrevPointOnBoundary)) {
1381
1382       if(!aListOfPointIndex.IsEmpty()) {
1383         nblines++;
1384         anArrayOfLines[nblines] = aListOfPointIndex;
1385         anArrayOfLineType[nblines] = bIsPrevPointOnBoundary;
1386         aListOfPointIndex.Clear();
1387       }
1388       bIsPrevPointOnBoundary = bIsCurrentPointOnBoundary;
1389     }
1390     aListOfPointIndex.Append(pit);
1391   } // for(pit=1; pit<=aNbPnts; pit++)
1392   //
1393   aNbListOfPointIndex=aListOfPointIndex.Extent();
1394   if(aNbListOfPointIndex) {
1395     nblines++;
1396     anArrayOfLines[nblines] = aListOfPointIndex;
1397     anArrayOfLineType[nblines] = bIsPrevPointOnBoundary;
1398     aListOfPointIndex.Clear();
1399   }
1400   //
1401   if(nblines <= 1){
1402     return Standard_False;
1403   }
1404   //
1405   // Correct wlines.begin
1406   Standard_Integer aLineType;
1407   TColStd_Array1OfListOfInteger anArrayOfLineEnds(1, nblines);
1408   Handle(IntSurf_LineOn2S) aSeqOfPntOn2S = new IntSurf_LineOn2S (new NCollection_IncAllocator());
1409   //
1410   for(i = 1; i <= nblines; i++) {
1411     aLineType=anArrayOfLineType[i];
1412     if(aLineType) {
1413       continue;
1414     }
1415     //
1416     const ListOfInteger& aListOfIndex = anArrayOfLines[i];
1417     if(aListOfIndex.Extent() < 2) {
1418       continue;
1419     }
1420     //
1421     TColStd_ListOfInteger aListOfFLIndex;
1422     Standard_Integer aneighbourindex, aLineTypeNeib;
1423     //
1424     for(j = 0; j < 2; j++) {// neighbour line choice
1425       aneighbourindex = (!j) ? (i-1) : (i+1);
1426       if((aneighbourindex < 1) || (aneighbourindex > nblines)){
1427         continue;
1428       }
1429       //
1430       aLineTypeNeib=anArrayOfLineType[aneighbourindex];
1431       if(!aLineTypeNeib){
1432         continue;
1433       }
1434       //
1435       const ListOfInteger& aNeighbour = anArrayOfLines[aneighbourindex];
1436       Standard_Integer anIndex = (!j) ? aNeighbour.Last() : aNeighbour.First();
1437       const IntSurf_PntOn2S& aPoint = theWLine->Point(anIndex);
1438       // check if need use derivative.begin .end [absence]
1439       //
1440       IntSurf_PntOn2S aNewP = aPoint;
1441       Standard_Integer surfit, parit;
1442       //
1443       for(surfit = 0; surfit < 2; ++surfit) {
1444
1445         Handle(GeomAdaptor_HSurface) aGASurface = (!surfit) ? theSurface1 : theSurface2;
1446
1447         aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
1448         Standard_Real U=0., V=0.;
1449
1450         if(!surfit) {
1451           aNewP.ParametersOnS1(U, V);
1452         }
1453         else {
1454           aNewP.ParametersOnS2(U, V);
1455         }
1456         //
1457         Standard_Integer nbboundaries = 0;
1458         Standard_Integer bIsUBoundary = Standard_False; // use if nbboundaries == 1
1459         Standard_Integer bIsFirstBoundary = Standard_False; // use if nbboundaries == 1
1460         //
1461         for(parit = 0; parit < 2; parit++) {
1462           Standard_Boolean isperiodic = (!parit) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
1463
1464           Standard_Real aResolution = (!parit) ? aGASurface->UResolution(aTol) : aGASurface->VResolution(aTol);
1465           Standard_Real alowerboundary = (!parit) ? umin : vmin;
1466           Standard_Real aupperboundary = (!parit) ? umax : vmax;
1467
1468           Standard_Real aParameter = (!parit) ? U : V;
1469           Standard_Boolean bIsOnFirstBoundary = Standard_True;
1470
1471           if(!isperiodic) {
1472             if(IsPointOnBoundary(aParameter, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary)) {
1473               bIsUBoundary = (!parit);
1474               bIsFirstBoundary = bIsOnFirstBoundary;
1475               nbboundaries++;
1476             }
1477           }
1478           else {
1479             Standard_Real aPeriod     = (!parit) ? aGASurface->UPeriod() : aGASurface->VPeriod();
1480             Standard_Real anoffset = 0.;
1481             Standard_Real anAdjustPar = AdjustPeriodic(aParameter, alowerboundary, aupperboundary, aPeriod, anoffset);
1482
1483             if(IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary)) {
1484               bIsUBoundary = (parit == 0);
1485               bIsFirstBoundary = bIsOnFirstBoundary;
1486               nbboundaries++;
1487             }
1488           }
1489         }
1490         //
1491         Standard_Boolean bComputeLineEnd = Standard_False;
1492
1493         if(nbboundaries == 2) {
1494           bComputeLineEnd = Standard_True;
1495         }
1496         else if(nbboundaries == 1) {
1497           Standard_Boolean isperiodic = (bIsUBoundary) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
1498
1499           if(isperiodic) {
1500             Standard_Real alowerboundary = (bIsUBoundary) ? umin : vmin;
1501             Standard_Real aupperboundary = (bIsUBoundary) ? umax : vmax;
1502             Standard_Real aPeriod     = (bIsUBoundary) ? aGASurface->UPeriod() : aGASurface->VPeriod();
1503             Standard_Real aParameter = (bIsUBoundary) ? U : V;
1504             Standard_Real anoffset = 0.;
1505             Standard_Real anAdjustPar = AdjustPeriodic(aParameter, alowerboundary, aupperboundary, aPeriod, anoffset);
1506
1507             Standard_Real adist = (bIsFirstBoundary) ? fabs(anAdjustPar - alowerboundary) : fabs(anAdjustPar - aupperboundary);
1508             Standard_Real anotherPar = (bIsFirstBoundary) ? (aupperboundary - adist) : (alowerboundary + adist);
1509             anotherPar += anoffset;
1510             Standard_Integer aneighbourpointindex = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
1511             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex);
1512             Standard_Real nU1, nV1;
1513
1514             if(surfit == 0)
1515               aNeighbourPoint.ParametersOnS1(nU1, nV1);
1516             else
1517               aNeighbourPoint.ParametersOnS2(nU1, nV1);
1518
1519             Standard_Real adist1 = (bIsUBoundary) ? fabs(nU1 - U) : fabs(nV1 - V);
1520             Standard_Real adist2 = (bIsUBoundary) ? fabs(nU1 - anotherPar) : fabs(nV1 - anotherPar);
1521             bComputeLineEnd = Standard_True;
1522             Standard_Boolean bCheckAngle1 = Standard_False;
1523             Standard_Boolean bCheckAngle2 = Standard_False;
1524             gp_Vec2d aNewVec;
1525             Standard_Real anewU = (bIsUBoundary) ? anotherPar : U;
1526             Standard_Real anewV = (bIsUBoundary) ? V : anotherPar;
1527             //
1528             if(((adist1 - adist2) > Precision::PConfusion()) &&
1529                (adist2 < (aPeriod / 4.))) {
1530               bCheckAngle1 = Standard_True;
1531               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(anewU, anewV));
1532
1533               if(aNewVec.SquareMagnitude() < (gp::Resolution() * gp::Resolution())) {
1534                 aNewP.SetValue((surfit == 0), anewU, anewV);
1535                 bCheckAngle1 = Standard_False;
1536               }
1537             }
1538             else if(adist1 < (aPeriod / 4.)) {
1539               bCheckAngle2 = Standard_True;
1540               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(U, V));
1541
1542               if(aNewVec.SquareMagnitude() < (gp::Resolution() * gp::Resolution())) {
1543                 bCheckAngle2 = Standard_False;
1544               }
1545             }
1546             //
1547             if(bCheckAngle1 || bCheckAngle2) {
1548               // assume there are at least two points in line (see "if" above)
1549               Standard_Integer anindexother = aneighbourpointindex;
1550
1551               while((anindexother <= aListOfIndex.Last()) && (anindexother >= aListOfIndex.First())) {
1552                 anindexother = (j == 0) ? (anindexother + 1) : (anindexother - 1);
1553                 const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(anindexother);
1554                 Standard_Real nU2, nV2;
1555
1556                 if(surfit == 0)
1557                   aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
1558                 else
1559                   aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
1560                 gp_Vec2d aVecOld(gp_Pnt2d(nU2, nV2), gp_Pnt2d(nU1, nV1));
1561
1562                 if(aVecOld.SquareMagnitude() <= (gp::Resolution() * gp::Resolution())) {
1563                   continue;
1564                 }
1565                 else {
1566                   Standard_Real anAngle = aNewVec.Angle(aVecOld);
1567
1568                   if((fabs(anAngle) < (M_PI * 0.25)) && (aNewVec.Dot(aVecOld) > 0.)) {
1569
1570                     if(bCheckAngle1) {
1571                       Standard_Real U1, U2, V1, V2;
1572                       IntSurf_PntOn2S atmppoint = aNewP;
1573                       atmppoint.SetValue((surfit == 0), anewU, anewV);
1574                       atmppoint.Parameters(U1, V1, U2, V2);
1575                       gp_Pnt P1 = theSurface1->Value(U1, V1);
1576                       gp_Pnt P2 = theSurface2->Value(U2, V2);
1577                       gp_Pnt P0 = aPoint.Value();
1578
1579                       if(P0.IsEqual(P1, aTol) &&
1580                          P0.IsEqual(P2, aTol) &&
1581                          P1.IsEqual(P2, aTol)) {
1582                         bComputeLineEnd = Standard_False;
1583                         aNewP.SetValue((surfit == 0), anewU, anewV);
1584                       }
1585                     }
1586
1587                     if(bCheckAngle2) {
1588                       bComputeLineEnd = Standard_False;
1589                     }
1590                   }
1591                   break;
1592                 }
1593               } // end while(anindexother...)
1594             }
1595           }
1596         }
1597         //
1598         if(bComputeLineEnd) {
1599           Standard_Integer aneighbourpointindex1 = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
1600           const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
1601           Standard_Real nU1, nV1;
1602
1603           if(surfit == 0)
1604             aNeighbourPoint.ParametersOnS1(nU1, nV1);
1605           else
1606             aNeighbourPoint.ParametersOnS2(nU1, nV1);
1607           gp_Pnt2d ap1(nU1, nV1);
1608           gp_Pnt2d ap2(nU1, nV1);
1609           Standard_Integer aneighbourpointindex2 = aneighbourpointindex1;
1610
1611           while((aneighbourpointindex2 <= aListOfIndex.Last()) && (aneighbourpointindex2 >= aListOfIndex.First())) {
1612             aneighbourpointindex2 = (j == 0) ? (aneighbourpointindex2 + 1) : (aneighbourpointindex2 - 1);
1613             const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(aneighbourpointindex2);
1614             Standard_Real nU2, nV2;
1615
1616             if(surfit == 0)
1617               aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
1618             else
1619               aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
1620             ap2.SetX(nU2);
1621             ap2.SetY(nV2);
1622
1623             if(ap1.SquareDistance(ap2) > (gp::Resolution() * gp::Resolution())) {
1624               break;
1625             }
1626           }
1627           gp_Pnt2d anewpoint;
1628           Standard_Boolean found = FindPoint(ap2, ap1, umin, umax, vmin, vmax, anewpoint);
1629
1630           if(found) {
1631             // check point
1632             Standard_Real aCriteria =aTolSum;// BRep_Tool::Tolerance(theFace1) + BRep_Tool::Tolerance(theFace2);
1633             //GeomAPI_ProjectPointOnSurf& aProjector = (surfit == 0) ? aPrj2 : aPrj1;
1634             ProjectPointOnSurf& aProjector = (surfit == 0) ? aPrj2 : aPrj1;
1635             Handle(GeomAdaptor_HSurface) aSurface = (surfit == 0) ? theSurface1 : theSurface2;
1636
1637             gp_Pnt aP3d = aSurface->Value(anewpoint.X(), anewpoint.Y());
1638             aProjector.Perform(aP3d);
1639
1640             if(aProjector.IsDone()) {
1641               if(aProjector.LowerDistance() < aCriteria) {
1642                 Standard_Real foundU = U, foundV = V;
1643                 aProjector.LowerDistanceParameters(foundU, foundV);
1644
1645                 if(surfit == 0)
1646                   aNewP.SetValue(aP3d, anewpoint.X(), anewpoint.Y(), foundU, foundV);
1647                 else
1648                   aNewP.SetValue(aP3d, foundU, foundV, anewpoint.X(), anewpoint.Y());
1649               }
1650             }
1651           }
1652         }
1653       }
1654       aSeqOfPntOn2S->Add(aNewP);
1655       aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
1656     }
1657     anArrayOfLineEnds.SetValue(i, aListOfFLIndex);
1658   }
1659   // Correct wlines.end
1660
1661   // Split wlines.begin
1662   for(j = 1; j <= theLConstructor.NbParts(); j++) {
1663     Standard_Real fprm = 0., lprm = 0.;
1664     theLConstructor.Part(j, fprm, lprm);
1665     Standard_Integer ifprm = (Standard_Integer)fprm;
1666     Standard_Integer ilprm = (Standard_Integer)lprm;
1667     //
1668     Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
1669     //
1670     for(i = 1; i <= nblines; i++) {
1671       if(anArrayOfLineType[i] != 0) {
1672         continue;
1673       }
1674       const ListOfInteger& aListOfIndex = anArrayOfLines[i];
1675
1676       if(aListOfIndex.Extent() < 2) {
1677         continue;
1678       }
1679       const TColStd_ListOfInteger& aListOfFLIndex = anArrayOfLineEnds.Value(i);
1680       Standard_Boolean bhasfirstpoint = (aListOfFLIndex.Extent() == 2);
1681       Standard_Boolean bhaslastpoint = (aListOfFLIndex.Extent() == 2);
1682
1683       if(!bhasfirstpoint && !aListOfFLIndex.IsEmpty()) {
1684         bhasfirstpoint = (i != 1);
1685       }
1686
1687       if(!bhaslastpoint && !aListOfFLIndex.IsEmpty()) {
1688         bhaslastpoint = (i != nblines);
1689       }
1690       Standard_Boolean bIsFirstInside = ((ifprm >= aListOfIndex.First()) && (ifprm <= aListOfIndex.Last()));
1691       Standard_Boolean bIsLastInside =  ((ilprm >= aListOfIndex.First()) && (ilprm <= aListOfIndex.Last()));
1692
1693       if(!bIsFirstInside && !bIsLastInside) {
1694         if((ifprm < aListOfIndex.First()) && (ilprm > aListOfIndex.Last())) {
1695           // append whole line, and boundaries if neccesary
1696           if(bhasfirstpoint) {
1697             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.First());
1698             aLineOn2S->Add(aP);
1699           }
1700           ListOfInteger::Iterator anIt(aListOfIndex);
1701
1702           for(; anIt.More(); anIt.Next()) {
1703             const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
1704             aLineOn2S->Add(aP);
1705           }
1706
1707           if(bhaslastpoint) {
1708             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.Last());
1709             aLineOn2S->Add(aP);
1710           }
1711
1712           // check end of split line (end is almost always)
1713           Standard_Integer aneighbour = i + 1;
1714           Standard_Boolean bIsEndOfLine = Standard_True;
1715
1716           if(aneighbour <= nblines) {
1717             const ListOfInteger& aListOfNeighbourIndex = anArrayOfLines[aneighbour];
1718
1719             if((anArrayOfLineType[aneighbour] != 0) &&
1720                (aListOfNeighbourIndex.IsEmpty())) {
1721               bIsEndOfLine = Standard_False;
1722             }
1723           }
1724
1725           if(bIsEndOfLine) {
1726             if(aLineOn2S->NbPoints() > 1) {
1727               Handle(IntPatch_WLine) aNewWLine =
1728                 new IntPatch_WLine(aLineOn2S, Standard_False);
1729               theNewLines.Append(aNewWLine);
1730             }
1731             aLineOn2S = new IntSurf_LineOn2S();
1732           }
1733         }
1734         continue;
1735       }
1736       // end if(!bIsFirstInside && !bIsLastInside)
1737
1738       if(bIsFirstInside && bIsLastInside) {
1739         // append inside points between ifprm and ilprm
1740         ListOfInteger::Iterator anIt(aListOfIndex);
1741
1742         for(; anIt.More(); anIt.Next()) {
1743           if((anIt.Value() < ifprm) || (anIt.Value() > ilprm))
1744             continue;
1745           const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
1746           aLineOn2S->Add(aP);
1747         }
1748       }
1749       else {
1750
1751         if(bIsFirstInside) {
1752           // append points from ifprm to last point + boundary point
1753           ListOfInteger::Iterator anIt(aListOfIndex);
1754
1755           for(; anIt.More(); anIt.Next()) {
1756             if(anIt.Value() < ifprm)
1757               continue;
1758             const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
1759             aLineOn2S->Add(aP);
1760           }
1761
1762           if(bhaslastpoint) {
1763             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.Last());
1764             aLineOn2S->Add(aP);
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         // end if(bIsFirstInside)
1789
1790         if(bIsLastInside) {
1791           // append points from first boundary point to ilprm
1792           if(bhasfirstpoint) {
1793             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.First());
1794             aLineOn2S->Add(aP);
1795           }
1796           ListOfInteger::Iterator anIt(aListOfIndex);
1797
1798           for(; anIt.More(); anIt.Next()) {
1799             if(anIt.Value() > ilprm)
1800               continue;
1801             const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
1802             aLineOn2S->Add(aP);
1803           }
1804         }
1805         //end if(bIsLastInside)
1806       }
1807     }
1808
1809     if(aLineOn2S->NbPoints() > 1) {
1810       Handle(IntPatch_WLine) aNewWLine =
1811         new IntPatch_WLine(aLineOn2S, Standard_False);
1812       theNewLines.Append(aNewWLine);
1813     }
1814   }
1815   // Split wlines.end
1816   //
1817   // cda002/I3
1818   Standard_Real fprm, lprm;
1819   Standard_Integer ifprm, ilprm, aNbPoints, aIndex;
1820   //
1821   aNbParts=theLConstructor.NbParts();
1822   //
1823   for(j = 1; j <= aNbParts; j++) {
1824     theLConstructor.Part(j, fprm, lprm);
1825     ifprm=(Standard_Integer)fprm;
1826     ilprm=(Standard_Integer)lprm;
1827     //
1828     if ((ilprm-ifprm)==1) {
1829       for(i = 1; i <= nblines; i++) {
1830         aLineType=anArrayOfLineType[i];
1831         if(aLineType) {
1832           continue;
1833         }
1834         //
1835         const ListOfInteger& aListOfIndex = anArrayOfLines[i];
1836         aNbPoints=aListOfIndex.Extent();
1837         if(aNbPoints==1) {
1838           aIndex=aListOfIndex.First();
1839           if (aIndex==ifprm || aIndex==ilprm) {
1840             Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
1841             const IntSurf_PntOn2S& aP1 = theWLine->Point(ifprm);
1842             const IntSurf_PntOn2S& aP2 = theWLine->Point(ilprm);
1843             aLineOn2S->Add(aP1);
1844             aLineOn2S->Add(aP2);
1845             Handle(IntPatch_WLine) aNewWLine =
1846               new IntPatch_WLine(aLineOn2S, Standard_False);
1847             theNewLines.Append(aNewWLine);
1848           }
1849         }
1850       }
1851     }
1852   }
1853   //
1854   return Standard_True;
1855 }
1856
1857 //=======================================================================
1858 //function : AdjustPeriodic
1859 //purpose  :
1860 //=======================================================================
1861 Standard_Real AdjustPeriodic(const Standard_Real theParameter,
1862                              const Standard_Real parmin,
1863                              const Standard_Real parmax,
1864                              const Standard_Real thePeriod,
1865                              Standard_Real&      theOffset)
1866 {
1867   Standard_Real aresult = theParameter;
1868   theOffset = 0.;
1869   while(aresult < parmin) {
1870     aresult += thePeriod;
1871     theOffset += thePeriod;
1872   }
1873   while(aresult > parmax) {
1874     aresult -= thePeriod;
1875     theOffset -= thePeriod;
1876   }
1877   return aresult;
1878 }
1879
1880 //=======================================================================
1881 //function : IsPointOnBoundary
1882 //purpose  :
1883 //=======================================================================
1884 Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
1885                                    const Standard_Real theFirstBoundary,
1886                                    const Standard_Real theSecondBoundary,
1887                                    const Standard_Real theResolution,
1888                                    Standard_Boolean&   IsOnFirstBoundary)
1889 {
1890   IsOnFirstBoundary = Standard_True;
1891   if(fabs(theParameter - theFirstBoundary) < theResolution)
1892     return Standard_True;
1893   if(fabs(theParameter - theSecondBoundary) < theResolution)
1894   {
1895     IsOnFirstBoundary = Standard_False;
1896     return Standard_True;
1897   }
1898   return Standard_False;
1899 }
1900 //=======================================================================
1901 //function : FindPoint
1902 //purpose  :
1903 //=======================================================================
1904 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
1905                            const gp_Pnt2d&     theLastPoint,
1906                            const Standard_Real theUmin,
1907                            const Standard_Real theUmax,
1908                            const Standard_Real theVmin,
1909                            const Standard_Real theVmax,
1910                            gp_Pnt2d&           theNewPoint)
1911 {
1912   gp_Vec2d aVec(theFirstPoint, theLastPoint);
1913   Standard_Integer i = 0, j = 0;
1914
1915   for(i = 0; i < 4; i++) {
1916     gp_Vec2d anOtherVec;
1917     gp_Vec2d anOtherVecNormal;
1918     gp_Pnt2d aprojpoint = theLastPoint;
1919
1920     if((i % 2) == 0) {
1921       anOtherVec.SetX(0.);
1922       anOtherVec.SetY(1.);
1923       anOtherVecNormal.SetX(1.);
1924       anOtherVecNormal.SetY(0.);
1925
1926       if(i < 2)
1927         aprojpoint.SetX(theUmin);
1928       else
1929         aprojpoint.SetX(theUmax);
1930     }
1931     else {
1932       anOtherVec.SetX(1.);
1933       anOtherVec.SetY(0.);
1934       anOtherVecNormal.SetX(0.);
1935       anOtherVecNormal.SetY(1.);
1936
1937       if(i < 2)
1938         aprojpoint.SetY(theVmin);
1939       else
1940         aprojpoint.SetY(theVmax);
1941     }
1942     gp_Vec2d anormvec = aVec;
1943     anormvec.Normalize();
1944     Standard_Real adot1 = anormvec.Dot(anOtherVecNormal);
1945
1946     if(fabs(adot1) < Precision::Angular())
1947       continue;
1948     Standard_Real adist = 0.;
1949
1950     if((i % 2) == 0) {
1951       adist = (i < 2) ? fabs(theLastPoint.X() - theUmin) : fabs(theLastPoint.X() - theUmax);
1952     }
1953     else {
1954       adist = (i < 2) ? fabs(theLastPoint.Y() - theVmin) : fabs(theLastPoint.Y() - theVmax);
1955     }
1956     Standard_Real anoffset = adist * anOtherVec.Dot(anormvec) / adot1;
1957
1958     for(j = 0; j < 2; j++) {
1959       anoffset = (j == 0) ? anoffset : -anoffset;
1960       gp_Pnt2d acurpoint(aprojpoint.XY() + (anOtherVec.XY()*anoffset));
1961       gp_Vec2d acurvec(theLastPoint, acurpoint);
1962
1963       //
1964       Standard_Real aDotX, anAngleX, aPC;
1965       //
1966       aDotX=aVec.Dot(acurvec);
1967       anAngleX=aVec.Angle(acurvec);
1968       aPC=Precision::PConfusion();
1969       //
1970       if(aDotX > 0. && fabs(anAngleX) < aPC) {
1971       //
1972         if((i % 2) == 0) {
1973           if((acurpoint.Y() >= theVmin) &&
1974              (acurpoint.Y() <= theVmax)) {
1975             theNewPoint = acurpoint;
1976             return Standard_True;
1977           }
1978         }
1979         else {
1980           if((acurpoint.X() >= theUmin) &&
1981              (acurpoint.X() <= theUmax)) {
1982             theNewPoint = acurpoint;
1983             return Standard_True;
1984           }
1985         }
1986       }
1987     }
1988   }
1989   return Standard_False;
1990 }
1991
1992 /*
1993 static
1994   void DumpWLine (const Handle(IntPatch_WLine)& theWLine);
1995 //=======================================================================
1996 //function : DumpWLine
1997 //purpose  :
1998 //=======================================================================
1999   void DumpWLine (const Handle(IntPatch_WLine)& theWLine)
2000 {
2001   Standard_Integer i, nbp;
2002   Standard_Real u1,v1,u2,v2;
2003   //
2004   nbp=theWLine->NbPnts();
2005   printf(" nbp=%d\n", nbp);
2006   for(i=1;i<=nbp;i++) {
2007     const IntSurf_PntOn2S& aPoint = theWLine->Point(i);
2008     const gp_Pnt& aP=aPoint.Value();
2009     aPoint.Parameters(u1,v1,u2,v2);
2010     printf("point i%d %lf %lf %lf [ %lf %lf ] [ %lf %lf ]\n",
2011            i, aP.X(), aP.Y(), aP.Z(), u1, v1, u2, v2);
2012   }
2013 }
2014 */