src/IntPatch/IntPatch_ImpPrmIntersection.cxx

   1 // Created on: 1992-05-07
   2 // Created by: Jacques GOUSSARD
   3 // Copyright (c) 1992-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 #include <IntPatch_ImpPrmIntersection.hxx>
  18
  19 #include <Adaptor3d_Surface.hxx>
  20 #include <Adaptor3d_TopolTool.hxx>
  21 #include <ElCLib.hxx>
  22 #include <ElSLib.hxx>
  23 #include <IntPatch_ArcFunction.hxx>
  24 #include <IntPatch_PointLine.hxx>
  25 #include <IntPatch_RLine.hxx>
  26 #include <IntPatch_RstInt.hxx>
  27 #include <IntPatch_SpecialPoints.hxx>
  28 #include <IntPatch_TheIWLineOfTheIWalking.hxx>
  29 #include <IntPatch_TheIWalking.hxx>
  30 #include <IntPatch_TheSurfFunction.hxx>
  31 #include <IntPatch_WLine.hxx>
  32 #include <IntSurf.hxx>
  33 #include <IntSurf_Quadric.hxx>
  34 #include <IntSurf_QuadricTool.hxx>
  35 #include <IntSurf_SequenceOfPathPoint.hxx>
  36 #include <TColStd_Array1OfInteger.hxx>
  37 #include <TopAbs_Orientation.hxx>
  38 #include <TopTrans_CurveTransition.hxx>
  39 #include <math_Matrix.hxx>
  40 #include <math_Vector.hxx>
  41
  42 #ifndef OCCT_DEBUG
  43 #define No_Standard_RangeError
  44 #define No_Standard_OutOfRange
  45 #endif
  46
  47 static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLine,
  48                                         const Standard_Boolean       IsReversed,
  49                                         const IntSurf_Quadric&       theQuad,
  50                                         const Handle(Adaptor3d_TopolTool)& thePDomain,
  51                                         const Handle(Adaptor3d_Surface)&  theQSurf,
  52                                         const Handle(Adaptor3d_Surface)&  theOtherSurf,
  53                                         const Standard_Real                theArcTol,
  54                                         const Standard_Real                theTolTang,
  55                                         IntPatch_SequenceOfLine&           theLines);
  56 static
  57   void ComputeTangency (const IntPatch_TheSOnBounds& solrst,
  58   IntSurf_SequenceOfPathPoint& seqpdep,
  59   const Handle(Adaptor3d_TopolTool)& Domain,
  60   IntPatch_TheSurfFunction& Func,
  61   const Handle(Adaptor3d_Surface)& PSurf,
  62   TColStd_Array1OfInteger& Destination);
  63 static
  64   void Recadre(const Standard_Boolean ,
  65   GeomAbs_SurfaceType typeS1,
  66   GeomAbs_SurfaceType typeS2,
  67   IntPatch_Point&  pt,
  68   const Handle(IntPatch_TheIWLineOfTheIWalking)& iwline,
  69   Standard_Integer Param,
  70   Standard_Real U1,
  71   Standard_Real V1,
  72   Standard_Real U2,
  73   Standard_Real V2);
  74
  75 static
  76   Standard_Boolean IsCoincide(IntPatch_TheSurfFunction& theFunc,
  77                               const Handle(IntPatch_PointLine)& theLine,
  78                               const Handle(Adaptor2d_Curve2d)& theArc,
  79                               const Standard_Boolean isTheSurface1Using,
  80                               const Standard_Real theToler3D,
  81                               const Standard_Real theToler2D,
  82                               const Standard_Real thePeriod);
  83
  84 static
  85   Standard_Real GetLocalStep(const Handle(Adaptor3d_Surface)& theSurf,
  86                              const Standard_Real theStep);
  87
  88
  89 //=======================================================================
  90 //function : IsSeamOrPole
  91 //purpose  :
  92 //=======================================================================
  93 static IntPatch_SpecPntType IsSeamOrPole(const Handle(Adaptor3d_Surface)& theQSurf,
  94                                          const Handle(IntSurf_LineOn2S)& theLine,
  95                                          const Standard_Boolean IsReversed,
  96                                          const Standard_Integer theRefIndex,
  97                                          const Standard_Real theTol3D,
  98                                          const Standard_Real theDeltaMax)
  99 {
 100   if((theRefIndex < 1) || (theRefIndex >= theLine->NbPoints()))
 101     return IntPatch_SPntNone;
 102
 103   //Parameters on Quadric and on parametric for reference point
 104   Standard_Real aUQRef, aVQRef, aUPRef, aVPRef;
 105   Standard_Real aUQNext, aVQNext, aUPNext, aVPNext;
 106
 107   const gp_Pnt &aP3d = theLine->Value(theRefIndex + 1).Value();
 108
 109   if(IsReversed)
 110   {
 111     theLine->Value(theRefIndex).Parameters  (aUPRef, aVPRef, aUQRef, aVQRef);
 112     theLine->Value(theRefIndex+1).Parameters(aUPNext, aVPNext, aUQNext, aVQNext);
 113   }
 114   else
 115   {
 116     theLine->Value(theRefIndex).Parameters  (aUQRef, aVQRef, aUPRef, aVPRef);
 117     theLine->Value(theRefIndex+1).Parameters(aUQNext, aVQNext, aUPNext, aVPNext);
 118   }
 119
 120   const GeomAbs_SurfaceType aType = theQSurf->GetType();
 121
 122   if ((aType == GeomAbs_Cone) &&
 123       (theQSurf->Cone().Apex().SquareDistance(aP3d) < theTol3D*theTol3D))
 124   {
 125     return IntPatch_SPntPoleSeamU;
 126   }
 127   else if (aType == GeomAbs_Sphere)
 128   {
 129     const Standard_Real aSqTol = theTol3D*theTol3D;
 130     gp_Pnt aP(ElSLib::Value(0.0, M_PI_2, theQSurf->Sphere()));
 131     if (aP.SquareDistance(aP3d) < aSqTol)
 132     {
 133       return IntPatch_SPntPoleSeamU;
 134     }
 135
 136     aP = ElSLib::Value(0.0, -M_PI_2, theQSurf->Sphere());
 137     if (aP.SquareDistance(aP3d) < aSqTol)
 138     {
 139       return IntPatch_SPntPoleSeamU;
 140     }
 141   }
 142
 143
 144   const Standard_Real aDeltaU = Abs(aUQRef - aUQNext);
 145
 146   if((aType != GeomAbs_Torus) && (aDeltaU < theDeltaMax))
 147     return IntPatch_SPntNone;
 148
 149   switch(aType)
 150   {
 151   case GeomAbs_Cylinder:
 152     return IntPatch_SPntSeamU;
 153
 154   case GeomAbs_Torus:
 155     {
 156       const Standard_Real aDeltaV = Abs(aVQRef - aVQNext);
 157
 158       if((aDeltaU >= theDeltaMax) && (aDeltaV >= theDeltaMax))
 159         return IntPatch_SPntSeamUV;
 160
 161       if(aDeltaU >= theDeltaMax)
 162         return IntPatch_SPntSeamU;
 163
 164       if(aDeltaV >= theDeltaMax)
 165         return IntPatch_SPntSeamV;
 166     }
 167
 168     break;
 169   case GeomAbs_Sphere:
 170   case GeomAbs_Cone:
 171     return IntPatch_SPntPoleSeamU;
 172   default:
 173     break;
 174   }
 175
 176   return IntPatch_SPntNone;
 177 }
 178
 179 //=======================================================================
 180 //function : IntPatch_ImpPrmIntersection
 181 //purpose  :
 182 //=======================================================================
 183 IntPatch_ImpPrmIntersection::IntPatch_ImpPrmIntersection ()
 184   : done(Standard_False),
 185   empt(Standard_False),
 186   myIsStartPnt(Standard_False),
 187   myUStart(0.0),
 188   myVStart(0.0)
 189 { }
 190
 191
 192 //=======================================================================
 193 //function : IntPatch_ImpPrmIntersection
 194 //purpose  :
 195 //=======================================================================
 196
 197 IntPatch_ImpPrmIntersection::IntPatch_ImpPrmIntersection
 198   (const Handle(Adaptor3d_Surface)&    Surf1,
 199   const Handle(Adaptor3d_TopolTool)&   D1,
 200   const Handle(Adaptor3d_Surface)&    Surf2,
 201   const Handle(Adaptor3d_TopolTool)&   D2,
 202   const Standard_Real    TolArc,
 203   const Standard_Real    TolTang,
 204   const Standard_Real    Fleche,
 205   const Standard_Real    Pas)
 206   : done(Standard_False),
 207   empt(Standard_False),
 208   myIsStartPnt(Standard_False),
 209   myUStart(0.0),
 210   myVStart(0.0)
 211 {
 212   Perform(Surf1,D1,Surf2,D2,TolArc,TolTang,Fleche,Pas);
 213 }
 214
 215
 216 //=======================================================================
 217 //function : SetStartPoint
 218 //purpose  :
 219 //=======================================================================
 220
 221 void IntPatch_ImpPrmIntersection::SetStartPoint(const Standard_Real U,
 222   const Standard_Real V)
 223 {
 224   myIsStartPnt = Standard_True;
 225   myUStart = U; myVStart = V;
 226 }
 227
 228 //=======================================================================
 229 //function : ComputeTangency
 230 //purpose  :
 231 //=======================================================================
 232 void ComputeTangency (const IntPatch_TheSOnBounds& solrst,
 233   IntSurf_SequenceOfPathPoint& seqpdep,
 234   const Handle(Adaptor3d_TopolTool)& Domain,
 235   IntPatch_TheSurfFunction& Func,
 236   const Handle(Adaptor3d_Surface)& PSurf,
 237   TColStd_Array1OfInteger& Destination)
 238 {
 239   Standard_Integer i,k, NbPoints, seqlength;
 240   Standard_Real theparam,test;
 241   Standard_Boolean fairpt, ispassing;
 242   TopAbs_Orientation arcorien,vtxorien;
 243   Handle(Adaptor2d_Curve2d) thearc;
 244   Handle(Adaptor3d_HVertex) vtx,vtxbis;
 245   //Standard_Boolean ispassing;
 246   IntPatch_ThePathPointOfTheSOnBounds PStart;
 247   IntSurf_PathPoint PPoint;
 248   gp_Vec vectg;
 249   gp_Dir2d dirtg;
 250   gp_Pnt ptbid;
 251   gp_Vec d1u,d1v,v1,v2;
 252   gp_Pnt2d p2d;
 253   gp_Vec2d d2d;
 254   //
 255   double aX[2], aF[1], aD[1][2];
 256   math_Vector X(aX, 1, 2);
 257   math_Vector F(aF, 1, 1);
 258   math_Matrix D(aD, 1, 1, 1, 2);
 259   //
 260   seqlength = 0;
 261   NbPoints = solrst.NbPoints();
 262   for (i=1; i<= NbPoints; i++) {
 263     if (Destination(i) == 0) {
 264       PStart = solrst.Point(i);
 265       thearc   = PStart.Arc();
 266       theparam = PStart.Parameter();
 267       arcorien = Domain->Orientation(thearc);
 268       ispassing = (arcorien == TopAbs_INTERNAL ||
 269         arcorien == TopAbs_EXTERNAL);
 270
 271       thearc->D0(theparam,p2d);
 272       X(1) = p2d.X();
 273       X(2) = p2d.Y();
 274       PPoint.SetValue(PStart.Value(),X(1),X(2));
 275
 276       Func.Values(X,F,D);
 277       if (Func.IsTangent()) {
 278         PPoint.SetTangency(Standard_True);
 279         Destination(i) = seqlength+1;
 280         if (!PStart.IsNew()) {
 281           vtx = PStart.Vertex();
 282           for (k=i+1; k<=NbPoints; k++) {
 283             if (Destination(k) ==0) {
 284               PStart = solrst.Point(k);
 285               if (!PStart.IsNew()) {
 286                 vtxbis = PStart.Vertex();
 287                 if (Domain->Identical(vtx,vtxbis)) {
 288                   thearc   = PStart.Arc();
 289                   theparam = PStart.Parameter();
 290                   arcorien = Domain->Orientation(thearc);
 291                   ispassing = ispassing && (arcorien == TopAbs_INTERNAL ||
 292                     arcorien == TopAbs_EXTERNAL);
 293
 294                   thearc->D0(theparam,p2d);
 295                   PPoint.AddUV(p2d.X(),p2d.Y());
 296                   Destination(k) = seqlength+1;
 297                 }
 298               }
 299             }
 300           }
 301         }
 302         PPoint.SetPassing(ispassing);
 303         seqpdep.Append(PPoint);
 304         seqlength++;
 305       }
 306       else { // on a un point de depart potentiel
 307
 308         vectg = Func.Direction3d();
 309         dirtg = Func.Direction2d();
 310
 311         PSurf->D1(X(1),X(2),ptbid,d1u,d1v);
 312         thearc->D1(theparam,p2d,d2d);
 313         v2.SetLinearForm(d2d.X(),d1u,d2d.Y(),d1v);
 314         v1 = d1u.Crossed(d1v);
 315
 316         test = vectg.Dot(v1.Crossed(v2));
 317         if (PStart.IsNew()) {
 318           if ((test < 0. && arcorien == TopAbs_FORWARD) ||
 319             (test > 0. && arcorien == TopAbs_REVERSED)) {
 320               vectg.Reverse();
 321               dirtg.Reverse();
 322           }
 323           PPoint.SetDirections(vectg,dirtg);
 324           PPoint.SetPassing(ispassing);
 325           Destination(i) = seqlength+1;
 326           seqpdep.Append(PPoint);
 327           seqlength++;
 328         }
 329         else { // traiter la transition complexe
 330           gp_Dir bidnorm(1.,1.,1.);
 331           Standard_Real tole = 1.e-8;
 332           TopAbs_Orientation LocTrans;
 333           TopTrans_CurveTransition comptrans;
 334           comptrans.Reset(vectg,bidnorm,0.);
 335           if (arcorien == TopAbs_FORWARD ||
 336             arcorien == TopAbs_REVERSED) {
 337               // pour essai
 338
 339               vtx = PStart.Vertex();
 340               vtxorien = Domain->Orientation(vtx);
 341               if (Abs(test) <= tole) {
 342                 LocTrans = TopAbs_EXTERNAL; // et pourquoi pas INTERNAL
 343               }
 344               else {
 345                 if (((test > 0.)&& arcorien == TopAbs_FORWARD) ||
 346                   ((test < 0.)&& arcorien == TopAbs_REVERSED)){
 347                     LocTrans = TopAbs_FORWARD;
 348                 }
 349                 else {
 350                   LocTrans = TopAbs_REVERSED;
 351                 }
 352                 if (arcorien == TopAbs_REVERSED) {v2.Reverse();}
 353               }
 354
 355               comptrans.Compare(tole,v2,bidnorm,0.,LocTrans,vtxorien);
 356           }
 357           Destination(i) = seqlength+1;
 358           for (k= i+1; k<=NbPoints; k++) {
 359             if (Destination(k) == 0) {
 360               PStart = solrst.Point(k);
 361               if (!PStart.IsNew()) {
 362                 vtxbis = PStart.Vertex();
 363                 if (Domain->Identical(vtx,vtxbis)) {
 364                   thearc   = PStart.Arc();
 365                   theparam = PStart.Parameter();
 366                   arcorien = Domain->Orientation(thearc);
 367
 368                   PPoint.AddUV(X(1),X(2));
 369
 370                   thearc->D1(theparam,p2d,d2d);
 371                   PPoint.AddUV(p2d.X(),p2d.Y());
 372
 373                   if (arcorien == TopAbs_FORWARD ||
 374                     arcorien == TopAbs_REVERSED) {
 375                       ispassing = Standard_False;
 376                       v2.SetLinearForm(d2d.X(),d1u,d2d.Y(),d1v);
 377
 378                       test = vectg.Dot(v1.Crossed(v2));
 379                       vtxorien = Domain->Orientation(PStart.Vertex());
 380                       if (Abs(test) <= tole) {
 381                         LocTrans = TopAbs_EXTERNAL; // et pourquoi pas INTERNAL
 382                       }
 383                       else {
 384                         if (((test > 0.)&& arcorien == TopAbs_FORWARD) ||
 385                           ((test < 0.)&& arcorien == TopAbs_REVERSED)){
 386                             LocTrans = TopAbs_FORWARD;
 387                         }
 388                         else {
 389                           LocTrans = TopAbs_REVERSED;
 390                         }
 391                         if (arcorien == TopAbs_REVERSED) {v2.Reverse();}
 392                       }
 393
 394                       comptrans.Compare(tole,v2,bidnorm,0.,LocTrans,vtxorien);
 395                   }
 396                   Destination(k) = seqlength+1;
 397                 }
 398               }
 399             }
 400           }
 401           fairpt = Standard_True;
 402           if (!ispassing) {
 403             TopAbs_State Before = comptrans.StateBefore();
 404             TopAbs_State After  = comptrans.StateAfter();
 405             if ((Before == TopAbs_UNKNOWN)||(After == TopAbs_UNKNOWN)) {
 406               fairpt = Standard_False;
 407             }
 408             else if (Before == TopAbs_IN) {
 409               if (After == TopAbs_IN) {
 410                 ispassing = Standard_True;
 411               }
 412               else {
 413                 vectg.Reverse();
 414                 dirtg.Reverse();
 415               }
 416             }
 417             else {
 418               if (After !=TopAbs_IN) {
 419                 fairpt = Standard_False;
 420               }
 421             }
 422           }
 423           if (fairpt) {
 424             PPoint.SetDirections(vectg,dirtg);
 425             PPoint.SetPassing(ispassing);
 426             seqpdep.Append(PPoint);
 427             seqlength++;
 428           }
 429           else { // il faut remettre en "ordre" si on ne garde pas le point.
 430             for (k=i; k <=NbPoints ; k++) {
 431               if (Destination(k)==seqlength + 1) {
 432                 Destination(k) = -Destination(k);
 433               }
 434             }
 435           }
 436         }
 437       }
 438     }
 439   }
 440 }
 441 //=======================================================================
 442 //function : Recadre
 443 //purpose  :
 444 //=======================================================================
 445 void Recadre(const Standard_Boolean ,
 446   GeomAbs_SurfaceType typeS1,
 447   GeomAbs_SurfaceType typeS2,
 448   IntPatch_Point&  pt,
 449   const Handle(IntPatch_TheIWLineOfTheIWalking)& iwline,
 450   Standard_Integer Param,
 451   Standard_Real U1,
 452   Standard_Real V1,
 453   Standard_Real U2,
 454   Standard_Real V2)
 455 {
 456   Standard_Real U1p,V1p,U2p,V2p;
 457   iwline->Line()->Value(Param).Parameters(U1p,V1p,U2p,V2p);
 458   switch(typeS1)
 459   {
 460   case GeomAbs_Torus:
 461     while(V1<(V1p-1.5*M_PI)) V1+=M_PI+M_PI;
 462     while(V1>(V1p+1.5*M_PI)) V1-=M_PI+M_PI;
 463     Standard_FALLTHROUGH
 464   case GeomAbs_Cylinder:
 465   case GeomAbs_Cone:
 466   case GeomAbs_Sphere:
 467     while(U1<(U1p-1.5*M_PI)) U1+=M_PI+M_PI;
 468     while(U1>(U1p+1.5*M_PI)) U1-=M_PI+M_PI;
 469   default:
 470     break;
 471   }
 472   switch(typeS2)
 473   {
 474   case GeomAbs_Torus:
 475     while(V2<(V2p-1.5*M_PI)) V2+=M_PI+M_PI;
 476     while(V2>(V2p+1.5*M_PI)) V2-=M_PI+M_PI;
 477     Standard_FALLTHROUGH
 478   case GeomAbs_Cylinder:
 479   case GeomAbs_Cone:
 480   case GeomAbs_Sphere:
 481     while(U2<(U2p-1.5*M_PI)) U2+=M_PI+M_PI;
 482     while(U2>(U2p+1.5*M_PI)) U2-=M_PI+M_PI;
 483   default:
 484     break;
 485   }
 486   pt.SetParameters(U1,V1,U2,V2);
 487 }
 488
 489 //=======================================================================
 490 //function : GetLocalStep
 491 //purpose  :
 492 //=======================================================================
 493 Standard_Real GetLocalStep(const Handle(Adaptor3d_Surface)& theSurf,
 494                                   const Standard_Real theStep)
 495 {
 496   Standard_Real aLocalStep = theStep;
 497   if (theSurf->UContinuity() > GeomAbs_C0 && theSurf->VContinuity() > GeomAbs_C0)
 498   {
 499     GeomAbs_SurfaceType aSType = theSurf->GetType();
 500
 501     if (aSType == GeomAbs_BezierSurface || aSType == GeomAbs_BSplineSurface)
 502     {
 503       Standard_Real aMinRes = Precision::Infinite();
 504       Standard_Integer aMaxDeg = 0;
 505       const Standard_Real aLimRes = 1.e-10;
 506
 507       aMinRes = Min(theSurf->UResolution(Precision::Confusion()),
 508                     theSurf->VResolution(Precision::Confusion()));
 509       aMaxDeg = Max(theSurf->UDegree(), theSurf->VDegree());
 510       if (aMinRes < aLimRes && aMaxDeg > 3)
 511       {
 512         aLocalStep = 0.0001;
 513       }
 514     }
 515   }
 516   if (theSurf->UContinuity() == GeomAbs_C0)
 517   {
 518     Standard_Integer aNbInt = theSurf->NbUIntervals(GeomAbs_C1);
 519     if (aNbInt > 1)
 520     {
 521       TColStd_Array1OfReal anInts(1, aNbInt + 1);
 522       theSurf->UIntervals(anInts, GeomAbs_C1);
 523       Standard_Integer i;
 524       Standard_Real aMinInt = Precision::Infinite();
 525       for (i = 1; i <= aNbInt; ++i)
 526       {
 527         aMinInt = Min(aMinInt, anInts(i + 1) - anInts(i));
 528       }
 529
 530       aMinInt /= theSurf->LastUParameter() - theSurf->FirstUParameter();
 531       if (aMinInt < 0.002)
 532       {
 533         aLocalStep = 0.0001;
 534       }
 535     }
 536
 537   }
 538
 539   if (theSurf->VContinuity() == GeomAbs_C0)
 540   {
 541     Standard_Integer aNbInt = theSurf->NbVIntervals(GeomAbs_C1);
 542     if (aNbInt > 1)
 543     {
 544       TColStd_Array1OfReal anInts(1, aNbInt + 1);
 545       theSurf->VIntervals(anInts, GeomAbs_C1);
 546       Standard_Integer i;
 547       Standard_Real aMinInt = Precision::Infinite();
 548       for (i = 1; i <= aNbInt; ++i)
 549       {
 550         aMinInt = Min(aMinInt, anInts(i + 1) - anInts(i));
 551       }
 552
 553       aMinInt /= theSurf->LastVParameter() - theSurf->FirstVParameter();
 554       if (aMinInt < 0.002)
 555       {
 556         aLocalStep = 0.0001;
 557       }
 558     }
 559   }
 560
 561   aLocalStep = Min(theStep, aLocalStep);
 562   return aLocalStep;
 563 }
 564 //=======================================================================
 565 //function : Perform
 566 //purpose  :
 567 //=======================================================================
 568 void IntPatch_ImpPrmIntersection::Perform (const Handle(Adaptor3d_Surface)& Surf1,
 569   const Handle(Adaptor3d_TopolTool)& D1,
 570   const Handle(Adaptor3d_Surface)& Surf2,
 571   const Handle(Adaptor3d_TopolTool)& D2,
 572   const Standard_Real TolArc,
 573   const Standard_Real TolTang,
 574   const Standard_Real Fleche,
 575   const Standard_Real Pas)
 576 {
 577   Standard_Boolean reversed, procf, procl, dofirst, dolast;
 578   Standard_Integer indfirst = 0, indlast = 0, ind2, NbSegm;
 579   Standard_Integer NbPointIns, NbPointRst, Nblines, Nbpts, NbPointDep;
 580   Standard_Real U1,V1,U2,V2,paramf,paraml,currentparam;
 581
 582   IntPatch_TheSegmentOfTheSOnBounds thesegm;
 583   IntSurf_PathPoint PPoint;
 584
 585   Handle(IntPatch_RLine) rline;
 586   Handle(IntPatch_WLine) wline;
 587   IntPatch_ThePathPointOfTheSOnBounds PStart,PStartf,PStartl;
 588   IntPatch_Point ptdeb,ptfin,ptbis;
 589
 590   IntPatch_IType typ;
 591   IntSurf_Transition TLine,TArc;
 592   IntSurf_TypeTrans trans1,trans2;
 593   gp_Pnt valpt,ptbid;
 594   gp_Vec tgline,tgrst,norm1,norm2,d1u,d1v;
 595   gp_Dir DirNormale;
 596   gp_Vec VecNormale;
 597
 598   gp_Pnt2d p2d;
 599   gp_Vec2d d2d;
 600
 601   Handle(Adaptor2d_Curve2d) currentarc;
 602   GeomAbs_SurfaceType typeS1, typeS2;
 603   IntSurf_Quadric Quad;
 604   IntPatch_TheSurfFunction Func;
 605   IntPatch_ArcFunction AFunc;
 606   //
 607   typeS1 = Surf1->GetType();
 608   typeS2 = Surf2->GetType();
 609
 610   paramf =0.;
 611   paraml =0.;
 612   trans1 = IntSurf_Undecided;
 613   trans2 = IntSurf_Undecided;
 614   //
 615   done = Standard_False;
 616   empt = Standard_True;
 617   slin.Clear();
 618   spnt.Clear();
 619   //
 620   reversed = Standard_False;
 621   switch (typeS1)
 622   {
 623   case GeomAbs_Plane:
 624     Quad.SetValue(Surf1->Plane());
 625     break;
 626
 627   case GeomAbs_Cylinder:
 628     Quad.SetValue(Surf1->Cylinder());
 629     break;
 630
 631   case GeomAbs_Sphere:
 632     Quad.SetValue(Surf1->Sphere());
 633     break;
 634
 635   case GeomAbs_Cone:
 636     Quad.SetValue(Surf1->Cone());
 637     break;
 638
 639   default:
 640     {
 641       reversed = Standard_True;
 642       switch (typeS2)
 643       {
 644       case GeomAbs_Plane:
 645         Quad.SetValue(Surf2->Plane());
 646         break;
 647
 648       case GeomAbs_Cylinder:
 649         Quad.SetValue(Surf2->Cylinder());
 650         break;
 651
 652       case GeomAbs_Sphere:
 653         Quad.SetValue(Surf2->Sphere());
 654         break;
 655
 656       case GeomAbs_Cone:
 657         Quad.SetValue(Surf2->Cone());
 658         break;
 659       default:
 660         {
 661           throw Standard_ConstructionError();
 662           break;
 663         }
 664       }
 665     }
 666     break;
 667   }
 668   //
 669   Standard_Real aLocalPas = Pas;
 670   if (reversed)
 671     aLocalPas = GetLocalStep(Surf1, Pas);
 672   else
 673     aLocalPas = GetLocalStep(Surf2, Pas);
 674
 675   Func.SetImplicitSurface(Quad);
 676   Func.Set(IntSurf_QuadricTool::Tolerance(Quad));
 677   AFunc.SetQuadric(Quad);
 678   //
 679   if (!reversed) {
 680     Func.Set(Surf2);
 681     AFunc.Set(Surf2);
 682   }
 683   else {
 684     Func.Set(Surf1);
 685     AFunc.Set(Surf1);
 686   }
 687   //
 688   if (!reversed) {
 689     solrst.Perform(AFunc,D2,TolArc,TolTang);
 690   }
 691   else {
 692     solrst.Perform(AFunc,D1,TolArc,TolTang);
 693   }
 694   if (!solrst.IsDone()) {
 695     return;
 696   }
 697   //
 698   IntSurf_SequenceOfPathPoint seqpdep;
 699   IntSurf_SequenceOfInteriorPoint seqpins;
 700   //
 701   NbPointRst = solrst.NbPoints();
 702   TColStd_Array1OfInteger Destination(1,NbPointRst+1); Destination.Init(0);
 703   if (NbPointRst) {
 704     if (!reversed) {
 705       ComputeTangency(solrst,seqpdep,D2,Func,Surf2,Destination);
 706     }
 707     else {
 708       ComputeTangency(solrst,seqpdep,D1,Func,Surf1,Destination);
 709     }
 710   }
 711   //
 712   Standard_Boolean SearchIns = Standard_True;
 713   if(Quad.TypeQuadric() == GeomAbs_Plane && solrst.NbSegments() > 0)
 714   {
 715     //For such kind of cases it is possible that whole surface is on one side of plane,
 716     //plane only touches surface and does not cross it,
 717     //so no inner points exist.
 718     SearchIns = Standard_False;
 719     Handle(Adaptor3d_TopolTool) T;
 720     if(reversed)
 721     {
 722       T = D1;
 723     }
 724     else
 725     {
 726       T = D2;
 727     }
 728     Standard_Integer aNbSamples = 0;
 729     aNbSamples = T->NbSamples();
 730     gp_Pnt2d s2d;
 731     gp_Pnt s3d;
 732     Standard_Real aValf[1], aUVap[2];
 733     math_Vector Valf(aValf,1,1), UVap(aUVap,1,2);
 734     T->SamplePoint(1,s2d, s3d);
 735     UVap(1)=s2d.X();
 736     UVap(2)=s2d.Y();
 737     Func.Value(UVap,Valf);
 738     Standard_Real rvalf = Sign(1.,Valf(1));
 739     for(Standard_Integer i = 2; i <= aNbSamples; ++i)
 740     {
 741       T->SamplePoint(i,s2d, s3d);
 742       UVap(1)=s2d.X();
 743       UVap(2)=s2d.Y();
 744       Func.Value(UVap,Valf);
 745       if(rvalf * Valf(1) < 0.)
 746       {
 747         SearchIns = Standard_True;
 748         break;
 749       }
 750     }
 751   }
 752   // Recherche des points interieurs
 753   NbPointIns = 0;
 754   if(SearchIns) {
 755     if (!reversed) {
 756       if (myIsStartPnt)
 757         solins.Perform(Func,Surf2,myUStart,myVStart);
 758       else
 759         solins.Perform(Func,Surf2,D2,TolTang);
 760     }
 761     else {
 762       if (myIsStartPnt)
 763         solins.Perform(Func,Surf1,myUStart,myVStart);
 764       else
 765         solins.Perform(Func,Surf1,D1,TolTang);
 766     }
 767     NbPointIns = solins.NbPoints();
 768     for (Standard_Integer i=1; i <= NbPointIns; i++) {
 769       seqpins.Append(solins.Value(i));
 770     }
 771   }
 772   //
 773   NbPointDep=seqpdep.Length();
 774   //
 775   if (NbPointDep || NbPointIns) {
 776     IntPatch_TheIWalking iwalk(TolTang, Fleche, aLocalPas);
 777     iwalk.Perform(seqpdep, seqpins, Func, reversed ? Surf1 : Surf2, reversed);
 778
 779     if(!iwalk.IsDone()) {
 780       return;
 781     }
 782
 783     Standard_Real Vmin, Vmax, TolV = 1.e-14;
 784     if (!reversed) { //Surf1 is quadric
 785       Vmin = Surf1->FirstVParameter();
 786       Vmax = Surf1->LastVParameter();
 787     }
 788     else { //Surf2 is quadric
 789       Vmin = Surf2->FirstVParameter();
 790       Vmax = Surf2->LastVParameter();
 791     }
 792     //
 793     Nblines = iwalk.NbLines();
 794     for (Standard_Integer j=1; j<=Nblines; j++) {
 795       const Handle(IntPatch_TheIWLineOfTheIWalking)&  iwline  = iwalk.Value(j);
 796       const Handle(IntSurf_LineOn2S)&                 thelin  = iwline->Line();
 797
 798       Nbpts = thelin->NbPoints();
 799       if(Nbpts>=2) {
 800         Standard_Integer k = 0;
 801         tgline = iwline->TangentVector(k);
 802         if(k>=1 && k<=Nbpts) { } else { k=Nbpts>>1; }
 803         valpt = thelin->Value(k).Value();
 804
 805         if (!reversed) {
 806           thelin->Value(k).ParametersOnS2(U2,V2);
 807           norm1 = Quad.Normale(valpt);
 808           Surf2->D1(U2,V2,ptbid,d1u,d1v);
 809           norm2 = d1u.Crossed(d1v);
 810         }
 811         else {
 812           thelin->Value(k).ParametersOnS1(U2,V2);
 813           norm2 = Quad.Normale(valpt);
 814           Surf1->D1(U2,V2,ptbid,d1u,d1v);
 815           norm1 = d1u.Crossed(d1v);
 816         }
 817         if (tgline.DotCross(norm2,norm1) > 0.) {
 818           trans1 = IntSurf_Out;
 819           trans2 = IntSurf_In;
 820         }
 821         else {
 822           trans1 = IntSurf_In;
 823           trans2 = IntSurf_Out;
 824         }
 825
 826         //
 827         Standard_Real AnU1,AnU2,AnV2;
 828
 829         GeomAbs_SurfaceType typQuad = Quad.TypeQuadric();
 830         Standard_Boolean arecadr=Standard_False;
 831         valpt = thelin->Value(1).Value();
 832         Quad.Parameters(valpt,AnU1,V1);
 833
 834         if((V1 < Vmin) && (Vmin-V1 < TolV)) V1 = Vmin;
 835         if((V1 > Vmax) && (V1-Vmax < TolV)) V1 = Vmax;
 836
 837         if(reversed) {
 838           thelin->SetUV(1,Standard_False,AnU1,V1); //-- on va lire u2,v2
 839           thelin->Value(1).ParametersOnS1(AnU2,AnV2);
 840         }
 841         else {
 842           thelin->SetUV(1,Standard_True,AnU1,V1);  //-- on va lire u1,v1
 843           thelin->Value(1).ParametersOnS2(AnU2,AnV2);
 844         }
 845
 846         if(typQuad==GeomAbs_Cylinder ||
 847           typQuad==GeomAbs_Cone ||
 848           typQuad==GeomAbs_Sphere) {
 849             arecadr=Standard_True;
 850         }
 851         //
 852         for (k=2; k<=Nbpts; ++k) {
 853           valpt = thelin->Value(k).Value();
 854           Quad.Parameters(valpt,U1,V1);
 855           //
 856           if((V1 < Vmin) && (Vmin-V1 < TolV)) {
 857             V1 = Vmin;
 858           }
 859           if((V1 > Vmax) && (V1-Vmax < TolV)) {
 860             V1 = Vmax;
 861           }
 862           //
 863           if(arecadr) {
 864             //modified by NIZNHY-PKV Fri Mar 28 15:06:01 2008f
 865             Standard_Real aCf, aTwoPI;
 866             //
 867             aCf=0.;
 868             aTwoPI=M_PI+M_PI;
 869             if ((U1-AnU1) >  1.5*M_PI) {
 870               while ((U1-AnU1) > (1.5*M_PI+aCf*aTwoPI)) {
 871                 aCf=aCf+1.;
 872               }
 873               U1=U1-aCf*aTwoPI;
 874             }
 875             //
 876             else {
 877               while ((U1-AnU1) < (-1.5*M_PI-aCf*aTwoPI)) {
 878                 aCf=aCf+1.;
 879               }
 880               U1=U1+aCf*aTwoPI;
 881             }
 882             // was:
 883             //if ((U1-AnU1) >  1.5*M_PI) {
 884             //  U1-=M_PI+M_PI;
 885             //}
 886             //else if ((U1-AnU1) < -1.5*M_PI) {
 887             //  U1+=M_PI+M_PI;
 888             //}
 889             //modified by NIZNHY-PKV Fri Mar 28 15:06:11 2008t
 890           }
 891           //
 892           if(reversed) {
 893             thelin->SetUV(k,Standard_False,U1,V1);
 894
 895             thelin->Value(k).ParametersOnS1(U2,V2);
 896             switch(typeS1) {
 897             case GeomAbs_Cylinder:
 898             case GeomAbs_Cone:
 899             case GeomAbs_Sphere:
 900             case GeomAbs_Torus:
 901               while(U2<(AnU2-1.5*M_PI)) U2+=M_PI+M_PI;
 902               while(U2>(AnU2+1.5*M_PI)) U2-=M_PI+M_PI;
 903               break;
 904             default:
 905               break;
 906             }
 907             if(typeS2==GeomAbs_Torus) {
 908               while(V2<(AnV2-1.5*M_PI)) V2+=M_PI+M_PI;
 909               while(V2>(AnV2+1.5*M_PI)) V2-=M_PI+M_PI;
 910             }
 911             thelin->SetUV(k,Standard_True,U2,V2);
 912           }
 913           else {
 914             thelin->SetUV(k,Standard_True,U1,V1);
 915
 916             thelin->Value(k).ParametersOnS2(U2,V2);
 917             switch(typeS2) {
 918             case GeomAbs_Cylinder:
 919             case GeomAbs_Cone:
 920             case GeomAbs_Sphere:
 921             case GeomAbs_Torus:
 922               while(U2<(AnU2-1.5*M_PI)) U2+=M_PI+M_PI;
 923               while(U2>(AnU2+1.5*M_PI)) U2-=M_PI+M_PI;
 924               break;
 925             default:
 926               break;
 927             }
 928             if(typeS2==GeomAbs_Torus) {
 929               while(V2<(AnV2-1.5*M_PI)) V2+=M_PI+M_PI;
 930               while(V2>(AnV2+1.5*M_PI)) V2-=M_PI+M_PI;
 931             }
 932             thelin->SetUV(k,Standard_False,U2,V2);
 933
 934           }
 935
 936           AnU1=U1;
 937           AnU2=U2;
 938           AnV2=V2;
 939         }
 940         // <-A
 941         wline = new IntPatch_WLine(thelin,Standard_False,trans1,trans2);
 942         wline->SetCreatingWayInfo(IntPatch_WLine::IntPatch_WLImpPrm);
 943
 944 #ifdef INTPATCH_IMPPRMINTERSECTION_DEBUG
 945         wline->Dump(0);
 946 #endif
 947
 948         if (   iwline->HasFirstPoint()
 949           && iwline->IsTangentAtBegining() == Standard_False)
 950         {
 951           indfirst = iwline->FirstPointIndex();
 952           PPoint = seqpdep(indfirst);
 953           tgline = PPoint.Direction3d();
 954           Standard_Integer themult = PPoint.Multiplicity();
 955           for (Standard_Integer i=NbPointRst; i>=1; i--) {
 956             if (Destination(i) == indfirst) {
 957               if (!reversed) { //-- typeS1 = Pln || Cyl || Sph || Cone
 958                 Quad.Parameters(PPoint.Value(),U1,V1);
 959
 960                 if((V1 < Vmin) && (Vmin-V1 < TolV)) V1 = Vmin;
 961                 if((V1 > Vmax) && (V1-Vmax < TolV)) V1 = Vmax;
 962
 963                 PPoint.Parameters(themult,U2,V2);
 964                 Surf2->D1(U2,V2,ptbid,d1u,d1v); //-- @@@@
 965               }
 966               else {  //-- typeS1 != Pln && Cyl && Sph && Cone
 967                 Quad.Parameters(PPoint.Value(),U2,V2);
 968
 969                 if((V2 < Vmin) && (Vmin-V2 < TolV)) V2 = Vmin;
 970                 if((V2 > Vmax) && (V2-Vmax < TolV)) V2 = Vmax;
 971
 972                 PPoint.Parameters(themult,U1,V1);
 973                 Surf1->D1(U1,V1,ptbid,d1u,d1v); //-- @@@@
 974               }
 975
 976               VecNormale = d1u.Crossed(d1v);
 977               //-- Modif du 27 Septembre 94 (Recadrage des pts U,V)
 978               ptdeb.SetValue(PPoint.Value(),TolArc,Standard_False);
 979               ptdeb.SetParameters(U1,V1,U2,V2);
 980               ptdeb.SetParameter(1.);
 981
 982               Recadre(reversed,typeS1,typeS2,ptdeb,iwline,1,U1,V1,U2,V2);
 983
 984               currentarc = solrst.Point(i).Arc();
 985               currentparam = solrst.Point(i).Parameter();
 986               currentarc->D1(currentparam,p2d,d2d);
 987               tgrst.SetLinearForm(d2d.X(),d1u,d2d.Y(),d1v);
 988
 989               Standard_Real squaremagnitudeVecNormale = VecNormale.SquareMagnitude();
 990               if(squaremagnitudeVecNormale > 1e-13) {
 991                 DirNormale=VecNormale;
 992                 IntSurf::MakeTransition(tgline,tgrst,DirNormale,TLine,TArc);
 993               }
 994               else {
 995                 TLine.SetValue(Standard_True,IntSurf_Undecided);
 996                 TArc.SetValue(Standard_True,IntSurf_Undecided);
 997               }
 998
 999               ptdeb.SetArc(reversed,currentarc,currentparam,TLine,TArc);
1000               if (!solrst.Point(i).IsNew()) {
1001                 ptdeb.SetVertex(reversed,solrst.Point(i).Vertex());
1002               }
1003               wline->AddVertex(ptdeb);
1004               if (themult == 0) {
1005                 wline->SetFirstPoint(wline->NbVertex());
1006               }
1007
1008               themult--;
1009             }
1010           }
1011         }
1012         else if (iwline->IsTangentAtBegining())
1013         {
1014           gp_Pnt psol = thelin->Value(1).Value();
1015           thelin->Value(1).ParametersOnS1(U1,V1);
1016           thelin->Value(1).ParametersOnS2(U2,V2);
1017           ptdeb.SetValue(psol,TolArc,Standard_True);
1018           ptdeb.SetParameters(U1,V1,U2,V2);
1019           ptdeb.SetParameter(1.);
1020           wline->AddVertex(ptdeb);
1021           wline->SetFirstPoint(wline->NbVertex());
1022         }
1023         else
1024         {
1025           gp_Pnt psol = thelin->Value(1).Value();
1026           thelin->Value(1).ParametersOnS1(U1,V1);
1027           thelin->Value(1).ParametersOnS2(U2,V2);
1028           ptdeb.SetValue(psol,TolArc,Standard_False);
1029           ptdeb.SetParameters(U1,V1,U2,V2);
1030           ptdeb.SetParameter(1.);
1031           wline->AddVertex(ptdeb);
1032           wline->SetFirstPoint(wline->NbVertex());
1033         }
1034
1035
1036         if (   iwline->HasLastPoint()
1037           && iwline->IsTangentAtEnd() == Standard_False)
1038         {
1039           indlast = iwline->LastPointIndex();
1040           PPoint = seqpdep(indlast);
1041           tgline = PPoint.Direction3d().Reversed();
1042           Standard_Integer themult = PPoint.Multiplicity();
1043           for (Standard_Integer i=NbPointRst; i >=1; i--) {
1044             if (Destination(i) == indlast) {
1045               if (!reversed) {
1046                 Quad.Parameters(PPoint.Value(),U1,V1);
1047
1048                 if((V1 < Vmin) && (Vmin-V1 < TolV)) V1 = Vmin;
1049                 if((V1 > Vmax) && (V1-Vmax < TolV)) V1 = Vmax;
1050
1051                 PPoint.Parameters(themult,U2,V2);
1052                 Surf2->D1(U2,V2,ptbid,d1u,d1v); //-- @@@@
1053                 VecNormale = d1u.Crossed(d1v);                    //-- @@@@
1054               }
1055               else {
1056                 Quad.Parameters(PPoint.Value(),U2,V2);
1057
1058                 if((V2 < Vmin) && (Vmin-V2 < TolV)) V2 = Vmin;
1059                 if((V2 > Vmax) && (V2-Vmax < TolV)) V2 = Vmax;
1060
1061                 PPoint.Parameters(themult,U1,V1);
1062                 Surf1->D1(U1,V1,ptbid,d1u,d1v); //-- @@@@
1063                 VecNormale = d1u.Crossed(d1v);                    //-- @@@@
1064               }
1065
1066               ptfin.SetValue(PPoint.Value(),TolArc,Standard_False);
1067               ptfin.SetParameters(U1,V1,U2,V2);
1068               ptfin.SetParameter(Nbpts);
1069
1070               Recadre(reversed,typeS1,typeS2,ptfin,iwline,Nbpts-1,U1,V1,U2,V2);
1071
1072               currentarc = solrst.Point(i).Arc();
1073               currentparam = solrst.Point(i).Parameter();
1074               currentarc->D1(currentparam,p2d,d2d);
1075               tgrst.SetLinearForm(d2d.X(),d1u,d2d.Y(),d1v);
1076
1077
1078               Standard_Real squaremagnitudeVecNormale = VecNormale.SquareMagnitude();
1079               if(squaremagnitudeVecNormale > 1e-13) {
1080                 DirNormale=VecNormale;
1081                 IntSurf::MakeTransition(tgline,tgrst,DirNormale,TLine,TArc);
1082               }
1083               else {
1084                 TLine.SetValue(Standard_True,IntSurf_Undecided);
1085                 TArc.SetValue(Standard_True,IntSurf_Undecided);
1086               }
1087
1088
1089               ptfin.SetArc(reversed,currentarc,currentparam,TLine,TArc);
1090               if (!solrst.Point(i).IsNew()) {
1091                 ptfin.SetVertex(reversed,solrst.Point(i).Vertex());
1092               }
1093               wline->AddVertex(ptfin);
1094               if (themult == 0) {
1095                 wline->SetLastPoint(wline->NbVertex());
1096               }
1097
1098               themult--;
1099             }
1100           }
1101         }
1102         else if (iwline->IsTangentAtEnd())
1103         {
1104           gp_Pnt psol = thelin->Value(Nbpts).Value();
1105           thelin->Value(Nbpts).ParametersOnS1(U1,V1);
1106           thelin->Value(Nbpts).ParametersOnS2(U2,V2);
1107           ptfin.SetValue(psol,TolArc,Standard_True);
1108           ptfin.SetParameters(U1,V1,U2,V2);
1109           ptfin.SetParameter(Nbpts);
1110           wline->AddVertex(ptfin);
1111           wline->SetLastPoint(wline->NbVertex());
1112         }
1113         else
1114         {
1115           gp_Pnt psol = thelin->Value(Nbpts).Value();
1116           thelin->Value(Nbpts).ParametersOnS1(U1,V1);
1117           thelin->Value(Nbpts).ParametersOnS2(U2,V2);
1118           ptfin.SetValue(psol,TolArc,Standard_False);
1119           ptfin.SetParameters(U1,V1,U2,V2);
1120           ptfin.SetParameter(Nbpts);
1121           wline->AddVertex(ptfin);
1122           wline->SetLastPoint(wline->NbVertex());
1123         }
1124         //
1125         // Il faut traiter les points de passage.
1126         slin.Append(wline);
1127       }// if(Nbpts>=2) {
1128     }// for (j=1; j<=Nblines; j++) {
1129
1130     // ON GERE LES RACCORDS ENTRE LIGNES. ELLE NE PEUVENT SE RACCORDER
1131     // QUE SUR DES POINTS DE TANGENCE
1132
1133
1134     Nblines = slin.Length();
1135     for (Standard_Integer j=1; j<=Nblines-1; j++) {
1136       dofirst = dolast = Standard_False;
1137       const  Handle(IntPatch_Line)& slinj = slin(j);
1138       Handle(IntPatch_WLine) wlin1 (Handle(IntPatch_WLine)::DownCast (slinj));
1139       if (wlin1->HasFirstPoint()) {
1140         ptdeb = wlin1->FirstPoint(indfirst);
1141         if (ptdeb.IsTangencyPoint()) {
1142           dofirst = Standard_True;
1143         }
1144       }
1145       if (wlin1->HasLastPoint()) {
1146         ptfin = wlin1->LastPoint(indlast);
1147         if (ptfin.IsTangencyPoint()) {
1148           dolast = Standard_True;
1149         }
1150       }
1151
1152       if (dofirst || dolast) {
1153         for (Standard_Integer k=j+1; k<=Nblines;k++) {
1154           const  Handle(IntPatch_Line)& slink = slin(k);
1155           Handle(IntPatch_WLine) wlin2 (Handle(IntPatch_WLine)::DownCast (slink));
1156           if (wlin2->HasFirstPoint()) {
1157             ptbis = wlin2->FirstPoint(ind2);
1158             if (ptbis.IsTangencyPoint()) {
1159               if (dofirst ) {
1160                 if (ptdeb.Value().Distance(ptbis.Value()) <= TolArc) {
1161                   ptdeb.SetMultiple(Standard_True);
1162                   if (!ptbis.IsMultiple()) {
1163                     ptbis.SetMultiple(Standard_True);
1164                     wlin2->Replace(ind2,ptbis);
1165                   }
1166                 }
1167               }
1168               if (dolast ) {
1169                 if (ptfin.Value().Distance(ptbis.Value()) <= TolArc) {
1170                   ptfin.SetMultiple(Standard_True);
1171                   if (!ptbis.IsMultiple()) {
1172                     ptbis.SetMultiple(Standard_True);
1173                     wlin2->Replace(ind2,ptbis);
1174                   }
1175                 }
1176               }
1177             }
1178           }
1179           if (wlin2->HasLastPoint()) {
1180             ptbis = wlin2->LastPoint(ind2);
1181             if (ptbis.IsTangencyPoint()) {
1182               if (dofirst ) {
1183                 if (ptdeb.Value().Distance(ptbis.Value()) <= TolArc) {
1184                   ptdeb.SetMultiple(Standard_True);
1185                   if (!ptbis.IsMultiple()) {
1186                     ptbis.SetMultiple(Standard_True);
1187                     wlin2->Replace(ind2,ptbis);
1188                   }
1189                 }
1190               }
1191               if (dolast ) {
1192                 if (ptfin.Value().Distance(ptbis.Value()) <= TolArc) {
1193                   ptfin.SetMultiple(Standard_True);
1194                   if (!ptbis.IsMultiple()) {
1195                     ptbis.SetMultiple(Standard_True);
1196                     wlin2->Replace(ind2,ptbis);
1197                   }
1198                 }
1199               }
1200             }
1201           }
1202         }
1203         if(dofirst)
1204           wlin1->Replace(indfirst,ptdeb);
1205         if(dolast)
1206           wlin1->Replace(indlast,ptfin);
1207       }
1208     }
1209   }// if (seqpdep.Length() != 0 || seqpins.Length() != 0) {
1210   //
1211   // Treatment the segments
1212   NbSegm = solrst.NbSegments();
1213   if (NbSegm) {
1214     for(Standard_Integer i=1; i<=NbSegm; i++) {
1215       thesegm = solrst.Segment(i);
1216       //Check if segment is degenerated
1217       if(thesegm.HasFirstPoint() && thesegm.HasLastPoint())
1218       {
1219         Standard_Real tol2 = Precision::Confusion();
1220         tol2 *= tol2;
1221         const gp_Pnt& aPf = thesegm.FirstPoint().Value();
1222         const gp_Pnt& aPl = thesegm.LastPoint().Value();
1223         if(aPf.SquareDistance(aPl) <= tol2)
1224         {
1225           //segment can be degenerated - check inner point
1226           paramf = thesegm.FirstPoint().Parameter();
1227           paraml = thesegm.LastPoint().Parameter();
1228           gp_Pnt2d _p2d =
1229             thesegm.Curve()->Value(.57735 * paramf + 0.42265 * paraml);
1230           gp_Pnt aPm;
1231           if(reversed)
1232           {
1233             Surf1->D0(_p2d.X(), _p2d.Y(), aPm);
1234           }
1235           else
1236           {
1237             Surf2->D0(_p2d.X(), _p2d.Y(), aPm);
1238           }
1239           if(aPm.SquareDistance(aPf) <= tol2)
1240           {
1241             //Degenerated
1242             continue;
1243           }
1244         }
1245       }
1246
1247
1248       //----------------------------------------------------------------------
1249       // on cree une ligne d intersection contenant uniquement le segment.
1250       // VOIR POUR LA TRANSITION DE LA LIGNE
1251       // On ajoute aussi un polygone pour le traitement des intersections
1252       // entre ligne et restrictions de la surface implicite (PutVertexOnLine)
1253       //----------------------------------------------------------------------
1254       //-- Calcul de la transition sur la rline (12 fev 97)
1255       //-- reversed a le sens de OnFirst
1256       //--
1257       dofirst = dolast  = Standard_False;
1258       procf = Standard_False;
1259       procl = Standard_False;
1260       IntSurf_Transition TLineUnk,TArcUnk;
1261
1262       IntPatch_Point _thepointAtBeg;
1263       IntPatch_Point _thepointAtEnd;
1264
1265       Standard_Boolean TransitionOK=Standard_False;
1266
1267       if(thesegm.HasFirstPoint()) {
1268         Standard_Real _u1,_v1,_u2,_v2;
1269
1270         dofirst = Standard_True;
1271         PStartf = thesegm.FirstPoint();
1272         paramf = PStartf.Parameter();
1273
1274         gp_Pnt2d _p2d     = thesegm.Curve()->Value(paramf);
1275         Handle(Adaptor3d_HVertex) _vtx;
1276         if(PStartf.IsNew()==Standard_False)
1277           _vtx= PStartf.Vertex();
1278         const gp_Pnt& _Pp = PStartf.Value();
1279         _thepointAtBeg.SetValue(_Pp,PStartf.Tolerance(),Standard_False);
1280         if (!reversed) { //-- typeS1 = Pln || Cyl || Sph || Cone
1281           Quad.Parameters(_Pp,_u1,_v1);
1282           _u2=_p2d.X(); _v2=_p2d.Y();
1283         }
1284         else {  //-- typeS1 != Pln && Cyl && Sph && Cone
1285           Quad.Parameters(_Pp,_u2,_v2);
1286           _u1=_p2d.X(); _v1=_p2d.Y();
1287         }
1288         _thepointAtBeg.SetParameters(_u1,_v1,_u2,_v2);
1289         _thepointAtBeg.SetParameter(paramf);
1290         if(PStartf.IsNew()==Standard_False)
1291           _thepointAtBeg.SetVertex(reversed,_vtx);
1292         _thepointAtBeg.SetArc(reversed,thesegm.Curve(),paramf,TLineUnk,TArcUnk);
1293
1294
1295         gp_Vec d1u1,d1v1,d1u2,d1v2; gp_Vec2d _d2d;
1296         Surf1->D1(_u1,_v1,ptbid,d1u1,d1v1);
1297         norm1 = d1u1.Crossed(d1v1);
1298         Surf2->D1(_u2,_v2,ptbid,d1u2,d1v2);
1299         norm2 = d1u2.Crossed(d1v2);
1300
1301         thesegm.Curve()->D1(paramf,_p2d,_d2d);
1302         if(reversed) {
1303           tgline.SetLinearForm(_d2d.X(),d1u1,_d2d.Y(),d1v1);
1304         }
1305         else {
1306           tgline.SetLinearForm(_d2d.X(),d1u2,_d2d.Y(),d1v2);
1307         }
1308         _u1=tgline.DotCross(norm2,norm1);
1309         TransitionOK=Standard_True;
1310         if (_u1 > 0.00000001) {
1311           trans1 = IntSurf_Out;
1312           trans2 = IntSurf_In;
1313         }
1314         else if(_u1 < -0.00000001) {
1315           trans1 = IntSurf_In;
1316           trans2 = IntSurf_Out;
1317         }
1318         else {
1319           TransitionOK=Standard_False;
1320         }
1321       }
1322       if(thesegm.HasLastPoint()) {
1323         Standard_Real _u1,_v1,_u2,_v2;
1324
1325         dolast = Standard_True;
1326         PStartl = thesegm.LastPoint();
1327         paraml = PStartl.Parameter();
1328
1329         gp_Pnt2d _p2d = thesegm.Curve()->Value(paraml);
1330         Handle(Adaptor3d_HVertex) _vtx;
1331         if(PStartl.IsNew()==Standard_False)
1332           _vtx = PStartl.Vertex();
1333         const gp_Pnt& _Pp = PStartl.Value();
1334         IntPatch_Point _thepoint;
1335         _thepointAtEnd.SetValue(_Pp,PStartl.Tolerance(),Standard_False);
1336         if (!reversed) { //-- typeS1 = Pln || Cyl || Sph || Cone
1337           Quad.Parameters(_Pp,_u1,_v1);
1338           _u2=_p2d.X(); _v2=_p2d.Y();
1339         }
1340         else {  //-- typeS1 != Pln && Cyl && Sph && Cone
1341           Quad.Parameters(_Pp,_u2,_v2);
1342           _u1=_p2d.X(); _v1=_p2d.Y();
1343         }
1344         _thepointAtEnd.SetParameters(_u1,_v1,_u2,_v2);
1345         _thepointAtEnd.SetParameter(paraml);
1346         if(PStartl.IsNew()==Standard_False)
1347           _thepointAtEnd.SetVertex(reversed,_vtx);
1348         _thepointAtEnd.SetArc(reversed,thesegm.Curve(),paraml,TLineUnk,TArcUnk);
1349
1350
1351
1352         gp_Vec d1u1,d1v1,d1u2,d1v2; gp_Vec2d _d2d;
1353         Surf1->D1(_u1,_v1,ptbid,d1u1,d1v1);
1354         norm1 = d1u1.Crossed(d1v1);
1355         Surf2->D1(_u2,_v2,ptbid,d1u2,d1v2);
1356         norm2 = d1u2.Crossed(d1v2);
1357
1358         thesegm.Curve()->D1(paraml,_p2d,_d2d);
1359         if(reversed) {
1360           tgline.SetLinearForm(_d2d.X(),d1u1,_d2d.Y(),d1v1);
1361         }
1362         else {
1363           tgline.SetLinearForm(_d2d.X(),d1u2,_d2d.Y(),d1v2);
1364         }
1365         _u1=tgline.DotCross(norm2,norm1);
1366         TransitionOK=Standard_True;
1367         if (_u1 > 0.00000001) {
1368           trans1 = IntSurf_Out;
1369           trans2 = IntSurf_In;
1370         }
1371         else if(_u1 < -0.00000001) {
1372           trans1 = IntSurf_In;
1373           trans2 = IntSurf_Out;
1374         }
1375         else {
1376           TransitionOK=Standard_False;
1377         }
1378       }
1379       if(TransitionOK==Standard_False) {
1380         //-- rline = new IntPatch_RLine (thesegm.Curve(),reversed,Standard_False);
1381         rline =  new IntPatch_RLine (Standard_False);
1382         if(reversed) {
1383           rline->SetArcOnS1(thesegm.Curve());
1384         }
1385         else {
1386           rline->SetArcOnS2(thesegm.Curve());
1387         }
1388       }
1389       else {
1390         //-- rline = new IntPatch_RLine (thesegm.Curve(),reversed,Standard_False,trans1,trans2);
1391         rline =  new IntPatch_RLine (Standard_False,trans1,trans2);
1392         if(reversed) {
1393           rline->SetArcOnS1(thesegm.Curve());
1394         }
1395         else {
1396           rline->SetArcOnS2(thesegm.Curve());
1397         }
1398       }
1399
1400       //------------------------------
1401       //-- Ajout des points
1402       //--
1403       if (thesegm.HasFirstPoint()) {
1404         rline->AddVertex(_thepointAtBeg);
1405         rline->SetFirstPoint(rline->NbVertex());
1406       }
1407
1408       if (thesegm.HasLastPoint()) {
1409         rline->AddVertex(_thepointAtEnd);
1410         rline->SetLastPoint(rline->NbVertex());
1411       }
1412
1413       // Polygone sur restriction solution
1414       if (dofirst && dolast) {
1415         Standard_Real prm;
1416         gp_Pnt ptpoly;
1417         IntSurf_PntOn2S p2s;
1418         Handle(IntSurf_LineOn2S) Thelin = new IntSurf_LineOn2S ();
1419         Handle(Adaptor2d_Curve2d) arcsegm = thesegm.Curve();
1420         Standard_Integer nbsample = 100;
1421
1422         if (!reversed) {
1423           for (Standard_Integer j=1; j<=nbsample; j++) {
1424             prm = paramf + (j-1)*(paraml-paramf)/(nbsample-1);
1425             arcsegm->D0(prm,p2d);
1426             Surf2->D0(p2d.X(),p2d.Y(),ptpoly);
1427
1428             Quad.Parameters(ptpoly,U1,V1);
1429             p2s.SetValue(ptpoly,U1,V1,p2d.X(),p2d.Y());
1430             Thelin->Add(p2s);
1431           }
1432         }
1433         else {
1434           for (Standard_Integer j=1; j<=nbsample; j++) {
1435             prm = paramf + (j-1)*(paraml-paramf)/(nbsample-1);
1436             arcsegm->D0(prm,p2d);
1437             Surf1->D0(p2d.X(),p2d.Y(),ptpoly);
1438
1439             Quad.Parameters(ptpoly,U2,V2);
1440             p2s.SetValue(ptpoly,p2d.X(),p2d.Y(),U2,V2);
1441             Thelin->Add(p2s);
1442           }
1443         }
1444         rline->Add(Thelin);
1445       }
1446
1447       if (dofirst || dolast) {
1448         Nblines = slin.Length();
1449         for (Standard_Integer j=1; j<=Nblines; j++) {
1450           const Handle(IntPatch_Line)& slinj = slin(j);
1451           typ = slinj->ArcType();
1452           if (typ == IntPatch_Walking) {
1453             Nbpts = Handle(IntPatch_WLine)::DownCast (slinj)->NbVertex();
1454           }
1455           else {
1456             Nbpts = Handle(IntPatch_RLine)::DownCast (slinj)->NbVertex();
1457           }
1458           for (Standard_Integer k=1; k<=Nbpts;k++) {
1459             if (typ == IntPatch_Walking) {
1460               ptdeb = Handle(IntPatch_WLine)::DownCast (slinj)->Vertex(k);
1461             }
1462             else {
1463               ptdeb = Handle(IntPatch_RLine)::DownCast (slinj)->Vertex(k);
1464             }
1465             if (dofirst) {
1466
1467               if (ptdeb.Value().Distance(PStartf.Value()) <=TolArc) {
1468                 ptdeb.SetMultiple(Standard_True);
1469                 if (typ == IntPatch_Walking) {
1470                   Handle(IntPatch_WLine)::DownCast (slinj)->Replace(k,ptdeb);
1471                 }
1472                 else {
1473                   Handle(IntPatch_RLine)::DownCast (slinj)->Replace(k,ptdeb);
1474                 }
1475                 ptdeb.SetParameter(paramf);
1476                 rline->AddVertex(ptdeb);
1477                 if (!procf){
1478                   procf=Standard_True;
1479                   rline->SetFirstPoint(rline->NbVertex());
1480                 }
1481               }
1482             }
1483             if (dolast) {
1484               if(dofirst) { //-- on recharge le ptdeb
1485                 if (typ == IntPatch_Walking) {
1486                   ptdeb = Handle(IntPatch_WLine)::DownCast (slinj)->Vertex(k);
1487                 }
1488                 else {
1489                   ptdeb = Handle(IntPatch_RLine)::DownCast (slinj)->Vertex(k);
1490                 }
1491               }
1492               if (ptdeb.Value().Distance(PStartl.Value()) <=TolArc) {
1493                 ptdeb.SetMultiple(Standard_True);
1494                 if (typ == IntPatch_Walking) {
1495                   Handle(IntPatch_WLine)::DownCast (slinj)->Replace(k,ptdeb);
1496                 }
1497                 else {
1498                   Handle(IntPatch_RLine)::DownCast (slinj)->Replace(k,ptdeb);
1499                 }
1500                 ptdeb.SetParameter(paraml);
1501                 rline->AddVertex(ptdeb);
1502                 if (!procl){
1503                   procl=Standard_True;
1504                   rline->SetLastPoint(rline->NbVertex());
1505                 }
1506               }
1507             }
1508           }
1509         }
1510       }
1511       slin.Append(rline);
1512     }
1513   }// if (NbSegm)
1514   //
1515   // on traite les restrictions de la surface implicite
1516
1517   for (Standard_Integer i=1, aNbLin = slin.Length(); i<=aNbLin; i++)
1518   {
1519     Handle(IntPatch_PointLine) aL = Handle(IntPatch_PointLine)::DownCast(slin(i));
1520
1521     if (!reversed)
1522       IntPatch_RstInt::PutVertexOnLine(aL,Surf1,D1,Surf2,Standard_True,TolTang);
1523     else
1524       IntPatch_RstInt::PutVertexOnLine(aL,Surf2,D2,Surf1,Standard_False,TolTang);
1525
1526     if (aL->NbPnts() <= 2)
1527     {
1528       Standard_Boolean aCond = aL->NbPnts() < 2;
1529       if (!aCond)
1530         aCond = (aL->Point(1).IsSame(aL->Point(2), Precision::Confusion()));
1531
1532       if (aCond)
1533       {
1534         slin.Remove(i);
1535         i--;
1536         aNbLin--;
1537         continue;
1538       }
1539     }
1540
1541     if(aL->ArcType() == IntPatch_Walking)
1542     {
1543       const Handle(IntPatch_WLine) aWL = Handle(IntPatch_WLine)::DownCast(aL);
1544       slin.Append(aWL);
1545       slin.Remove(i);
1546       i--;
1547       aNbLin--;
1548     }
1549   }
1550
1551   // Now slin is filled as follows: lower indices correspond to Restriction line,
1552   // after (higher indices) - only Walking-line.
1553
1554   const Standard_Real aTol3d = Max(Func.Tolerance(), TolTang);
1555   const Handle(Adaptor3d_Surface)& aQSurf = (reversed) ? Surf2 : Surf1;
1556   const Handle(Adaptor3d_Surface)& anOtherSurf = (reversed) ? Surf1 : Surf2;
1557
1558   for (Standard_Integer i = 1; i <= slin.Length(); i++)
1559   {
1560     const Handle(IntPatch_PointLine)& aL1 = Handle(IntPatch_PointLine)::DownCast(slin(i));
1561     const Handle(IntPatch_RLine)& aRL1 = Handle(IntPatch_RLine)::DownCast(aL1);
1562
1563     if(aRL1.IsNull())
1564     {
1565       //Walking-Walking cases are not supported
1566       break;
1567     }
1568
1569     const Handle(Adaptor2d_Curve2d)& anArc = aRL1->IsArcOnS1() ?
1570                                               aRL1->ArcOnS1() :
1571                                               aRL1->ArcOnS2();
1572     if(anArc->GetType() != GeomAbs_Line)
1573     {
1574       //Restriction line must be isoline.
1575       //Other cases are not supported by
1576       //existing algorithms.
1577
1578       break;
1579     }
1580
1581     Standard_Boolean isFirstDeleted = Standard_False;
1582
1583     for(Standard_Integer j = i + 1; j <= slin.Length(); j++)
1584     {
1585       Handle(IntPatch_PointLine) aL2 = Handle(IntPatch_PointLine)::DownCast(slin(j));
1586       Handle(IntPatch_RLine) aRL2 = Handle(IntPatch_RLine)::DownCast(aL2);
1587
1588       //Here aL1 (i-th line) is Restriction-line and aL2 (j-th line) is
1589       //Restriction or Walking
1590
1591       if(!aRL2.IsNull())
1592       {
1593         const Handle(Adaptor2d_Curve2d)& anArc2 = aRL2->IsArcOnS1() ?
1594                                                    aRL2->ArcOnS1() :
1595                                                    aRL2->ArcOnS2();
1596         if(anArc2->GetType() != GeomAbs_Line)
1597         {
1598           //Restriction line must be isoline.
1599           //Other cases are not supported by
1600           //existing algorithms.
1601
1602           continue;
1603         }
1604       }
1605
1606       //aDir can be equal to one of following four values only
1607       //(because Reastriction line is boundary of rectangular surface):
1608       //either {0, 1} or {0, -1} or {1, 0} or {-1, 0}.
1609       const gp_Dir2d aDir = anArc->Line().Direction();
1610
1611       Standard_Real aTol2d = anOtherSurf->UResolution(aTol3d),
1612                     aPeriod = anOtherSurf->IsVPeriodic() ? anOtherSurf->VPeriod() : 0.0;
1613
1614       if(Abs(aDir.X()) < 0.5)
1615       {//Restriction directs along V-direction
1616         aTol2d = anOtherSurf->VResolution(aTol3d);
1617         aPeriod = anOtherSurf->IsUPeriodic() ? anOtherSurf->UPeriod() : 0.0;
1618       }
1619
1620       const Standard_Boolean isCoincide = IsCoincide(Func, aL2, anArc, aRL1->IsArcOnS1(),
1621                                                       aTol3d, aTol2d, aPeriod);
1622
1623       if(isCoincide)
1624       {
1625         if(aRL2.IsNull())
1626         {//Delete Walking-line
1627           slin.Remove(j);
1628           j--;
1629         }
1630         else
1631         {//Restriction-Restriction
1632           const Handle(Adaptor2d_Curve2d)& anArc2 = aRL2->IsArcOnS1() ?
1633                                                      aRL2->ArcOnS1() :
1634                                                      aRL2->ArcOnS2();
1635
1636           const Standard_Real aRange2 = anArc2->LastParameter() -
1637                                         anArc2->FirstParameter();
1638           const Standard_Real aRange1 = anArc->LastParameter() -
1639                                         anArc->FirstParameter();
1640
1641           if(aRange2 > aRange1)
1642           {
1643             isFirstDeleted = Standard_True;
1644             break;
1645           }
1646           else
1647           {//Delete j-th line
1648             slin.Remove(j);
1649             j--;
1650           }
1651         }
1652       }
1653     } //for(Standard_Integer j = i + 1; j <= slin.Length(); j++)
1654
1655     if(isFirstDeleted)
1656     {//Delete i-th line
1657       slin.Remove(i--);
1658     }
1659   }//for (Standard_Integer i = 1; i <= slin.Length(); i++)
1660
1661   empt = (slin.Length() == 0 && spnt.Length() == 0);
1662   done = Standard_True;
1663
1664
1665   if(slin.Length() == 0)
1666     return;
1667
1668   Standard_Boolean isDecomposeRequired =  (Quad.TypeQuadric() == GeomAbs_Cone) ||
1669                                           (Quad.TypeQuadric() == GeomAbs_Sphere) ||
1670                                           (Quad.TypeQuadric() == GeomAbs_Cylinder) ||
1671                                           (Quad.TypeQuadric() == GeomAbs_Torus);
1672
1673   if(!isDecomposeRequired)
1674     return;
1675
1676   // post processing for cones and spheres
1677
1678   const Handle(Adaptor3d_TopolTool)& PDomain = (reversed) ? D1 : D2;
1679
1680   IntPatch_SequenceOfLine dslin;
1681   Standard_Boolean isDecompose = Standard_False;
1682   for(Standard_Integer i = 1; i <= slin.Length(); i++ )
1683   {
1684     if(DecomposeResult( Handle(IntPatch_PointLine)::DownCast(slin(i)),
1685                                         reversed, Quad, PDomain, aQSurf,
1686                                         anOtherSurf, TolArc, aTol3d, dslin))
1687     {
1688       isDecompose = Standard_True;
1689     }
1690   }
1691
1692   if(!isDecompose)
1693     return;
1694
1695   slin.Clear();
1696   for(Standard_Integer i = 1; i <= dslin.Length(); i++ )
1697     slin.Append(dslin(i));
1698 }
1699
1700 // correct U parameter of the start point of line on Quadric
1701 // (change 0->2PI or vs, if necessary)
1702 static Standard_Real AdjustUFirst(Standard_Real U1,Standard_Real U2)
1703 {
1704   Standard_Real u = U1;
1705
1706   // case: no adjustment
1707   if( U1 > 0. && U1 < (2.*M_PI) )
1708     return u;
1709
1710   // case: near '0'
1711   if( U1 == 0. || fabs(U1) <= 1.e-9 ) {
1712     if( U2 > 0. && U2 < (2.*M_PI) )
1713       u = ( U2 < ((2.*M_PI)-U2) ) ? 0. : (2.*M_PI);
1714     else {
1715       Standard_Real uu = U2;
1716       if( U2 > (2.*M_PI) )
1717         while( uu > (2.*M_PI) )
1718           uu -= (2.*M_PI);
1719       else
1720         while( uu < 0.)
1721           uu += (2.*M_PI);
1722
1723       u = ( uu < ((2.*M_PI)-uu) ) ? 0. : (2.*M_PI);
1724     }
1725   }
1726   // case: near '2PI'
1727   else if( U1 == (2.*M_PI) || fabs((2.*M_PI)-fabs(U1)) <= 1.e-9 ) {
1728     if( U2 > 0. && U2 < (2.*M_PI) )
1729       u = ( U2 < ((2.*M_PI)-U2) ) ? 0. : (2.*M_PI);
1730     else {
1731       Standard_Real uu = U2;
1732       if( U2 > (2.*M_PI) )
1733         while( uu > (2.*M_PI) )
1734           uu -= (2.*M_PI);
1735       else
1736         while( uu < 0.)
1737           uu += (2.*M_PI);
1738
1739       u = ( uu < ((2.*M_PI)-uu) ) ? 0. : (2.*M_PI);
1740     }
1741   }
1742   // case: '<0. || >2PI'
1743   else {
1744     if(U1 < 0.)
1745       while(u < 0.)
1746         u += 2.*M_PI;
1747     if(U1 > (2.*M_PI))
1748       while(u > (2.*M_PI))
1749         u -= (2.*M_PI);
1750   }
1751
1752   return u;
1753 }
1754
1755 // collect vertices, reject equals
1756 static Handle(IntSurf_LineOn2S) GetVertices(const Handle(IntPatch_PointLine)& thePLine,
1757                                             const Standard_Real           TOL3D,
1758                                             const Standard_Real           TOL2D)
1759 {
1760   //  Standard_Real TOL3D = 1.e-12, TOL2D = 1.e-8;
1761
1762   Handle(IntSurf_LineOn2S) vertices = new IntSurf_LineOn2S();
1763
1764   Standard_Real U1 = 0., U2 = 0., V1 = 0., V2 = 0.;
1765   Standard_Integer i = 0, k = 0;
1766   Standard_Integer NbVrt = thePLine->NbVertex();
1767
1768   TColStd_Array1OfInteger anVrts(1,NbVrt);
1769   anVrts.Init(0);
1770
1771   // check equal vertices
1772   for(i = 1; i <= NbVrt; i++) {
1773
1774     if( anVrts(i) == -1 ) continue;
1775
1776     const IntPatch_Point& Pi = thePLine->Vertex(i);
1777
1778     for(k = (i+1); k <= NbVrt; k++) {
1779
1780       if( anVrts(k) == -1 ) continue;
1781
1782       const IntPatch_Point& Pk = thePLine->Vertex(k);
1783
1784       if(Pi.Value().Distance(Pk.Value()) <= TOL3D) {
1785         // suggest the points are equal;
1786         // test 2d parameters on surface
1787         Standard_Boolean sameU1 = Standard_False;
1788         Standard_Boolean sameV1 = Standard_False;
1789         Standard_Boolean sameU2 = Standard_False;
1790         Standard_Boolean sameV2 = Standard_False;
1791
1792         Pi.ParametersOnS1(U1,V1);
1793         Pk.ParametersOnS1(U2,V2);
1794         if(fabs(U1-U2) <= TOL2D) sameU1 = Standard_True;
1795         if(fabs(V1-V2) <= TOL2D) sameV1 = Standard_True;
1796
1797         Pi.ParametersOnS2(U1,V1);
1798         Pk.ParametersOnS2(U2,V2);
1799         if(fabs(U1-U2) <= TOL2D) sameU2 = Standard_True;
1800         if(fabs(V1-V2) <= TOL2D) sameV2 = Standard_True;
1801
1802         if((sameU1 && sameV1) && (sameU2 && sameV2))
1803           anVrts(k) = -1;
1804       }
1805     }
1806   }
1807
1808   // copy further processed vertices
1809   for(i = 1; i <= NbVrt; i++) {
1810     if( anVrts(i) == -1 ) continue;
1811     vertices->Add(thePLine->Vertex(i).PntOn2S());
1812   }
1813   return vertices;
1814 }
1815
1816 static void SearchVertices(const Handle(IntSurf_LineOn2S)& Line,
1817   const Handle(IntSurf_LineOn2S)& Vertices,
1818   TColStd_Array1OfInteger&        PTypes)
1819 {
1820   Standard_Integer nbp = Line->NbPoints(), nbv = Vertices->NbPoints();
1821   Standard_Integer ip = 0, iv = 0;
1822   for(ip = 1; ip <= nbp; ip++) {
1823     const IntSurf_PntOn2S& aP = Line->Value(ip);
1824     Standard_Integer type = 0;
1825     for(iv = 1; iv <= nbv; iv++) {
1826       const IntSurf_PntOn2S& aV = Vertices->Value(iv);
1827       if(aP.IsSame(aV, Precision::Confusion(), Precision::PConfusion())) {
1828         type = iv;
1829         break;
1830       }
1831     }
1832     PTypes(ip) = type;
1833   }
1834 }
1835
1836 static inline Standard_Boolean IsSeamParameter(const Standard_Real U,
1837   const Standard_Real TOL2D)
1838 {
1839   return (fabs(U) <= TOL2D || fabs(2.*M_PI - U) <= TOL2D);
1840 }
1841
1842 static inline Standard_Real AdjustU(const Standard_Real U)
1843 {
1844   Standard_Real u = U, DBLPI = 2.*M_PI;
1845   if(u < 0. || u > DBLPI) {
1846     if(u < 0.)
1847       while(u < 0.)
1848         u += DBLPI;
1849     else
1850       while(u > DBLPI)
1851         u -= DBLPI;
1852   }
1853   return u;
1854 }
1855
1856 static inline void Correct2DBounds(const Standard_Real UF,
1857   const Standard_Real UL,
1858   const Standard_Real VF,
1859   const Standard_Real VL,
1860   const Standard_Real TOL2D,
1861   Standard_Real&      U,
1862   Standard_Real&      V)
1863 {
1864   Standard_Real Eps = 1.e-16;
1865   Standard_Real dUF = fabs(U - UF);
1866   Standard_Real dUL = fabs(U - UL);
1867   Standard_Real dVF = fabs(V - VF);
1868   Standard_Real dVL = fabs(V - VL);
1869   if(dUF <= TOL2D && dUF > Eps) U = UF;
1870   if(dUL <= TOL2D && dUL > Eps) U = UL;
1871   if(dVF <= TOL2D && dVF > Eps) V = VF;
1872   if(dVL <= TOL2D && dVL > Eps) V = VL;
1873 }
1874
1875 static void AdjustLine(Handle(IntSurf_LineOn2S)& Line,
1876   const Standard_Boolean    IsReversed,
1877   const Handle(Adaptor3d_Surface)&         QSurf,
1878   const Standard_Real       TOL2D)
1879 {
1880   Standard_Real VF = QSurf->FirstVParameter();
1881   Standard_Real VL = QSurf->LastVParameter();
1882   Standard_Real UF = QSurf->FirstUParameter();
1883   Standard_Real UL = QSurf->LastUParameter();
1884
1885   Standard_Integer nbp = Line->NbPoints(), ip = 0;
1886   Standard_Real U = 0., V = 0.;
1887   for(ip = 1; ip <= nbp; ip++) {
1888     if(IsReversed) {
1889       Line->Value(ip).ParametersOnS2(U,V);
1890       U = AdjustU(U);
1891       Correct2DBounds(UF,UL,VF,VL,TOL2D,U,V);
1892       Line->SetUV(ip,Standard_False,U,V);
1893     }
1894     else {
1895       Line->Value(ip).ParametersOnS1(U,V);
1896       U = AdjustU(U);
1897       Correct2DBounds(UF,UL,VF,VL,TOL2D,U,V);
1898       Line->SetUV(ip,Standard_True,U,V);
1899     }
1900   }
1901 }
1902
1903 static Standard_Boolean InsertSeamVertices(Handle(IntSurf_LineOn2S)&       Line,
1904   const Standard_Boolean          IsReversed,
1905   Handle(IntSurf_LineOn2S)&       Vertices,
1906   const TColStd_Array1OfInteger&  PTypes,
1907   const Standard_Real             TOL2D)
1908 {
1909   Standard_Boolean result = Standard_False;
1910   Standard_Integer ip = 0, nbp = Line->NbPoints();
1911   Standard_Real U = 0., V = 0.;
1912   for(ip = 1; ip <= nbp; ip++) {
1913     Standard_Integer ipt = PTypes(ip);
1914     if(ipt != 0) {
1915       const IntSurf_PntOn2S& aP = Line->Value(ip);
1916       if(IsReversed)
1917         aP.ParametersOnS2(U,V); // S2 - quadric
1918       else
1919         aP.ParametersOnS1(U,V); // S1 - quadric
1920       U = AdjustU(U);
1921       if(IsSeamParameter(U,TOL2D)) {
1922         if(ip == 1 || ip == nbp) {
1923           Standard_Real U1 = 0., V1 = 0.;
1924           Standard_Integer ipp = (ip == 1) ? (ip+1) : (ip-1);
1925           if(IsReversed)
1926             Line->Value(ipp).ParametersOnS2(U1,V1); // S2 - quadric
1927           else
1928             Line->Value(ipp).ParametersOnS1(U1,V1); // S1 - quadric
1929           Standard_Real u = AdjustUFirst(U,U1);
1930           if(fabs(u-U) >= 1.5*M_PI) {
1931             Standard_Real U2 = 0., V2 = 0.;
1932             if(IsReversed) {
1933               Line->Value(ip).ParametersOnS1(U2,V2); // prm
1934               Line->SetUV(ip,Standard_False,u,V);
1935               Line->SetUV(ip,Standard_True,U2,V2);
1936             }
1937             else {
1938               Line->Value(ip).ParametersOnS2(U2,V2); // prm
1939               Line->SetUV(ip,Standard_True,u,V);
1940               Line->SetUV(ip,Standard_False,U2,V2);
1941             }
1942           }
1943         }
1944         else {
1945           Standard_Integer ipp = ip - 1;
1946           Standard_Integer ipn = ip + 1;
1947           Standard_Real U1 = 0., V1 = 0., U2 = 0., V2 = 0.;
1948           if(IsReversed) {
1949             Line->Value(ipp).ParametersOnS2(U1,V1); // quad
1950             Line->Value(ipn).ParametersOnS2(U2,V2); // quad
1951           }
1952           else {
1953             Line->Value(ipp).ParametersOnS1(U1,V1); // quad
1954             Line->Value(ipn).ParametersOnS1(U2,V2); // quad
1955           }
1956           U1 = AdjustU(U1);
1957           U2 = AdjustU(U2);
1958           Standard_Boolean pnearZero = (fabs(U1) < fabs(2.*M_PI-U1)) ? Standard_True : Standard_False;
1959           Standard_Boolean cnearZero = (fabs(U) < fabs(2.*M_PI-U)) ? Standard_True : Standard_False;
1960           if(pnearZero == cnearZero) {
1961             if(!IsSeamParameter(U2,TOL2D) && !IsSeamParameter(U1,TOL2D)) {
1962               Standard_Real nU = (cnearZero) ? (2.*M_PI) : 0.;
1963               IntSurf_PntOn2S nP;
1964               nP.SetValue(aP.Value());
1965               Standard_Real U3 = 0., V3 = 0.;
1966               if(IsReversed) {
1967                 Line->Value(ip).ParametersOnS1(U3,V3); // prm
1968                 nP.SetValue(Standard_False,nU,V);
1969                 nP.SetValue(Standard_True,U3,V3);
1970               }
1971               else {
1972                 Line->Value(ip).ParametersOnS2(U3,V3); // prm
1973                 nP.SetValue(Standard_True,nU,V);
1974                 nP.SetValue(Standard_False,U3,V3);
1975               }
1976               Line->InsertBefore(ipn,nP);
1977               Vertices->Add(nP);
1978               result = Standard_True;
1979               break;
1980             }
1981           }
1982           else {
1983             if(!IsSeamParameter(U2,TOL2D) && !IsSeamParameter(U1,TOL2D)) {
1984               Standard_Real nU = (cnearZero) ? (2.*M_PI) : 0.;
1985               IntSurf_PntOn2S nP;
1986               nP.SetValue(aP.Value());
1987               Standard_Real U3 = 0., V3 = 0.;
1988               if(IsReversed) {
1989                 Line->Value(ip).ParametersOnS1(U3,V3); // prm
1990                 nP.SetValue(Standard_False,nU,V);
1991                 nP.SetValue(Standard_True,U3,V3);
1992               }
1993               else {
1994                 Line->Value(ip).ParametersOnS2(U3,V3); // prm
1995                 nP.SetValue(Standard_True,nU,V);
1996                 nP.SetValue(Standard_False,U3,V3);
1997               }
1998               Line->InsertBefore(ip,nP);
1999               Vertices->Add(nP);
2000               result = Standard_True;
2001               break;
2002             }
2003             else {
2004               // Line->InsertBefore(ip,Line->Value(ipn));
2005               // Line->RemovePoint(ip+2);
2006               // result = Standard_True;
2007               // std::cout << "swap vertex " << std::endl;
2008               // break;
2009             }
2010           }
2011         }
2012       }
2013     }
2014   }
2015   return result;
2016 }
2017
2018 static void ToSmooth( const Handle(IntSurf_LineOn2S)& Line,
2019                       const Standard_Boolean          IsReversed,
2020                       const IntSurf_Quadric&          Quad,
2021                       const Standard_Boolean          IsFirst,
2022                       Standard_Real&                  D3D)
2023 {
2024   if(Line->NbPoints() <= 10)
2025     return;
2026
2027   D3D = 0.;
2028   Standard_Integer NbTestPnts = Line->NbPoints() / 5;
2029   if(NbTestPnts < 5) NbTestPnts = 5;
2030
2031   Standard_Integer startp = (IsFirst) ? 2 : (Line->NbPoints() - NbTestPnts - 2);
2032   Standard_Integer ip = 0;
2033   Standard_Real Uc = 0., Vc = 0., Un = 0., Vn = 0., DDU = 0.;
2034   //Standard_Real DDV = 0.;
2035
2036   for(ip = startp; ip <= NbTestPnts; ip++) {
2037     if(IsReversed) {
2038       Line->Value(ip).ParametersOnS2(Uc,Vc); // S2 - quadric
2039       Line->Value(ip+1).ParametersOnS2(Un,Vn);
2040     }
2041     else {
2042       Line->Value(ip).ParametersOnS1(Uc,Vc); // S1 - quadric
2043       Line->Value(ip+1).ParametersOnS1(Un,Vn);
2044     }
2045     DDU += fabs(fabs(Uc)-fabs(Un));
2046     //DDV += fabs(fabs(Vc)-fabs(Vn));
2047
2048     if(ip > startp) {
2049       Standard_Real DP = Line->Value(ip).Value().Distance(Line->Value(ip-1).Value());
2050       D3D += DP;
2051     }
2052   }
2053
2054   DDU /= (Standard_Real) NbTestPnts + 1;
2055   //DDV /= (Standard_Real) NbTestPnts + 1;
2056
2057   D3D /= (Standard_Real) NbTestPnts + 1;
2058
2059
2060   Standard_Integer Index1 = (IsFirst) ? 1 : (Line->NbPoints());
2061   Standard_Integer Index2 = (IsFirst) ? 2 : (Line->NbPoints()-1);
2062   Standard_Integer Index3 = (IsFirst) ? 3 : (Line->NbPoints()-2);
2063
2064   Standard_Boolean doU = Standard_False;
2065
2066   Standard_Real U1 = 0., U2 = 0., V1 = 0., V2 = 0., U3 = 0., V3 = 0.;
2067
2068   if(IsReversed) {
2069     Line->Value(Index1).ParametersOnS2(U1,V1); // S2 - quadric
2070     Line->Value(Index2).ParametersOnS2(U2,V2);
2071     Line->Value(Index3).ParametersOnS2(U3,V3);
2072   }
2073   else {
2074     Line->Value(Index1).ParametersOnS1(U1,V1); // S1 - quadric
2075     Line->Value(Index2).ParametersOnS1(U2,V2);
2076     Line->Value(Index3).ParametersOnS1(U3,V3);
2077   }
2078
2079   if(!doU && Quad.TypeQuadric() == GeomAbs_Sphere) {
2080     if(fabs(fabs(U1)-fabs(U2)) > (M_PI/16.)) doU = Standard_True;
2081
2082     if(doU && (fabs(U1) <= 1.e-9 || fabs(U1-2.*M_PI) <= 1.e-9)) {
2083       if(fabs(V1-M_PI/2.) <= 1.e-9 || fabs(V1+M_PI/2.) <= 1.e-9) {}
2084       else {
2085         doU = Standard_False;
2086       }
2087     }
2088   }
2089
2090   if(Quad.TypeQuadric() == GeomAbs_Cone) {
2091     Standard_Real Uapx = 0., Vapx = 0.;
2092     Quad.Parameters(Quad.Cone().Apex(),Uapx,Vapx);
2093
2094     if(fabs(fabs(U1)-fabs(U2)) > M_PI/32.) doU = Standard_True;
2095
2096     if(doU && (fabs(U1) <= 1.e-9 || fabs(U1-2.*M_PI) <= 1.e-9)) {
2097       if(fabs(V1-Vapx) <= 1.e-9) {}
2098       else {
2099         doU = Standard_False;
2100       }
2101     }
2102   }
2103
2104   if(doU) {
2105     Standard_Real dU = Min((DDU/10.),5.e-8);
2106     Standard_Real U = (U2 > U3) ? (U2 + dU) : (U2 - dU);
2107     if(IsReversed)
2108       Line->SetUV(Index1,Standard_False,U,V1);
2109     else
2110       Line->SetUV(Index1,Standard_True,U,V1);
2111     U1 = U;
2112   }
2113 }
2114
2115 static Standard_Boolean TestMiddleOnPrm(const IntSurf_PntOn2S& aP,
2116                                         const IntSurf_PntOn2S& aV,
2117                                         const Standard_Boolean IsReversed,
2118                                         const Standard_Real    ArcTol,
2119                                         const Handle(Adaptor3d_TopolTool)&  PDomain)
2120
2121 {
2122   Standard_Boolean result = Standard_False;
2123   Standard_Real Up = 0., Vp = 0., Uv = 0., Vv = 0.;
2124   if(IsReversed) {
2125     aP.ParametersOnS1(Up,Vp); //S1 - parametric
2126     aV.ParametersOnS1(Uv,Vv);
2127   }
2128   else {
2129     aP.ParametersOnS2(Up,Vp); // S2 - parametric
2130     aV.ParametersOnS2(Uv,Vv);
2131   }
2132   Standard_Real Um = (Up + Uv)*0.5, Vm = (Vp + Vv)*0.5;
2133   gp_Pnt2d a2DPntM(Um,Vm);
2134   TopAbs_State PosM = PDomain->Classify(a2DPntM,ArcTol);
2135   if(PosM == TopAbs_ON || PosM == TopAbs_IN )
2136     result = Standard_True;
2137   return result;
2138 }
2139
2140 static void VerifyVertices( const Handle(IntSurf_LineOn2S)&    Line,
2141                             const Standard_Boolean             IsReversed,
2142                             const Handle(IntSurf_LineOn2S)&    Vertices,
2143                             const Standard_Real                TOL2D,
2144                             const Standard_Real       ArcTol,
2145                             const Handle(Adaptor3d_TopolTool)& PDomain,
2146                             IntSurf_PntOn2S&          VrtF,
2147                             Standard_Boolean&         AddFirst,
2148                             IntSurf_PntOn2S&          VrtL,
2149                             Standard_Boolean&         AddLast)
2150 {
2151   Standard_Integer nbp = Line->NbPoints(), nbv = Vertices->NbPoints();
2152   Standard_Integer FIndexSame = 0, FIndexNear = 0, LIndexSame = 0, LIndexNear = 0;
2153   const IntSurf_PntOn2S& aPF = Line->Value(1);
2154   const IntSurf_PntOn2S& aPL = Line->Value(nbp);
2155   Standard_Real UF = 0., VF = 0., UL = 0., VL = 0.;
2156   if(IsReversed) {
2157     aPF.ParametersOnS2(UF,VF);
2158     aPL.ParametersOnS2(UL,VL);
2159   }
2160   else {
2161     aPF.ParametersOnS1(UF,VF);
2162     aPL.ParametersOnS1(UL,VL);
2163   }
2164   gp_Pnt2d a2DPF(UF,VF);
2165   gp_Pnt2d a2DPL(UL,VL);
2166   Standard_Real DistMinF = 1.e+100, DistMinL = 1.e+100;
2167   Standard_Integer FConjugated = 0, LConjugated = 0;
2168
2169   Standard_Integer iv = 0;
2170
2171   for(iv = 1; iv <= nbv; iv++) {
2172     Standard_Real Uv = 0., Vv = 0.;
2173     if(IsReversed) {
2174       Vertices->Value(iv).ParametersOnS2(Uv,Vv);
2175       Uv = AdjustU(Uv);
2176       Vertices->SetUV(iv,Standard_False,Uv,Vv);
2177     }
2178     else {
2179       Vertices->Value(iv).ParametersOnS1(Uv,Vv);
2180       Uv = AdjustU(Uv);
2181       Vertices->SetUV(iv,Standard_True,Uv,Vv);
2182     }
2183   }
2184
2185   for(iv = 1; iv <= nbv; iv++) {
2186     const IntSurf_PntOn2S& aV = Vertices->Value(iv);
2187     if(aPF.IsSame(aV, Precision::Confusion(), Precision::PConfusion())) {
2188       FIndexSame = iv;
2189       break;
2190     }
2191     else {
2192       Standard_Real Uv = 0., Vv = 0.;
2193       if(IsReversed)
2194         aV.ParametersOnS2(Uv,Vv);
2195       else
2196         aV.ParametersOnS1(Uv,Vv);
2197       gp_Pnt2d a2DV(Uv,Vv);
2198       Standard_Real Dist = a2DV.Distance(a2DPF);
2199       if(Dist < DistMinF) {
2200         DistMinF = Dist;
2201         FIndexNear = iv;
2202         if(FConjugated != 0)
2203           FConjugated = 0;
2204       }
2205       if(IsSeamParameter(Uv,TOL2D)) {
2206         Standard_Real Ucv = (fabs(Uv) < fabs(2.*M_PI-Uv)) ? (2.*M_PI) : 0.;
2207         gp_Pnt2d a2DCV(Ucv,Vv);
2208         Standard_Real CDist = a2DCV.Distance(a2DPF);
2209         if(CDist < DistMinF) {
2210           DistMinF = CDist;
2211           FConjugated = iv;
2212           FIndexNear = iv;
2213         }
2214       }
2215     }
2216   }
2217
2218   for(iv = 1; iv <= nbv; iv++) {
2219     const IntSurf_PntOn2S& aV = Vertices->Value(iv);
2220     if(aPL.IsSame(aV, Precision::Confusion(), Precision::PConfusion())) {
2221       LIndexSame = iv;
2222       break;
2223     }
2224     else {
2225       Standard_Real Uv = 0., Vv = 0.;
2226       if(IsReversed)
2227         aV.ParametersOnS2(Uv,Vv);
2228       else
2229         aV.ParametersOnS1(Uv,Vv);
2230       gp_Pnt2d a2DV(Uv,Vv);
2231       Standard_Real Dist = a2DV.Distance(a2DPL);
2232       if(Dist < DistMinL) {
2233         DistMinL = Dist;
2234         LIndexNear = iv;
2235         if(LConjugated != 0)
2236           LConjugated = 0;
2237       }
2238       if(IsSeamParameter(Uv,TOL2D)) {
2239         Standard_Real Ucv = (fabs(Uv) < fabs(2.*M_PI-Uv)) ? (2.*M_PI) : 0.;
2240         gp_Pnt2d a2DCV(Ucv,Vv);
2241         Standard_Real CDist = a2DCV.Distance(a2DPL);
2242         if(CDist < DistMinL) {
2243           DistMinL = CDist;
2244           LConjugated = iv;
2245           LIndexNear = iv;
2246         }
2247       }
2248     }
2249   }
2250
2251   AddFirst = Standard_False;
2252   AddLast  = Standard_False;
2253
2254   if(FIndexSame == 0) {
2255     if(FIndexNear != 0) {
2256       const IntSurf_PntOn2S& aV = Vertices->Value(FIndexNear);
2257       Standard_Real Uv = 0., Vv = 0.;
2258       if(IsReversed)
2259         aV.ParametersOnS2(Uv,Vv);
2260       else
2261         aV.ParametersOnS1(Uv,Vv);
2262       if(IsSeamParameter(Uv,TOL2D)) {
2263         Standard_Real Ucv = (fabs(Uv) < fabs(2.*M_PI-Uv)) ? (2.*M_PI) : 0.;
2264         Standard_Boolean test = TestMiddleOnPrm(aPF,aV,IsReversed,ArcTol,PDomain);
2265         if(test) {
2266           VrtF.SetValue(aV.Value());
2267           if(IsReversed) {
2268             Standard_Real U2 = 0., V2 = 0.;
2269             aV.ParametersOnS1(U2,V2); // S1 - prm
2270             VrtF.SetValue(Standard_True,U2,V2);
2271             if(FConjugated == 0)
2272               VrtF.SetValue(Standard_False,Uv,Vv);
2273             else
2274               VrtF.SetValue(Standard_False,Ucv,Vv);
2275           }
2276           else {
2277             Standard_Real U2 = 0., V2 = 0.;
2278             aV.ParametersOnS2(U2,V2); // S2 - prm
2279             VrtF.SetValue(Standard_False,U2,V2);
2280             if(FConjugated == 0)
2281               VrtF.SetValue(Standard_True,Uv,Vv);
2282             else
2283               VrtF.SetValue(Standard_True,Ucv,Vv);
2284           }
2285           Standard_Real Dist3D = VrtF.Value().Distance(aPF.Value());
2286           if(Dist3D > 1.5e-7 && DistMinF > TOL2D) {
2287             AddFirst = Standard_True;
2288           }
2289         }
2290       }
2291       else {
2292         // to do: analyze internal vertex
2293       }
2294     }
2295   }
2296
2297   if(LIndexSame == 0) {
2298     if(LIndexNear != 0) {
2299       const IntSurf_PntOn2S& aV = Vertices->Value(LIndexNear);
2300       Standard_Real Uv = 0., Vv = 0.;
2301       if(IsReversed)
2302         aV.ParametersOnS2(Uv,Vv);
2303       else
2304         aV.ParametersOnS1(Uv,Vv);
2305       if(IsSeamParameter(Uv,TOL2D)) {
2306         Standard_Real Ucv = (fabs(Uv) < fabs(2.*M_PI-Uv)) ? (2.*M_PI) : 0.;
2307         Standard_Boolean test = TestMiddleOnPrm(aPL,aV,IsReversed,ArcTol,PDomain);
2308         if(test) {
2309           VrtL.SetValue(aV.Value());
2310           if(IsReversed) {
2311             Standard_Real U2 = 0., V2 = 0.;
2312             aV.ParametersOnS1(U2,V2); // S1 - prm
2313             VrtL.SetValue(Standard_True,U2,V2);
2314             if(LConjugated == 0)
2315               VrtL.SetValue(Standard_False,Uv,Vv);
2316             else
2317               VrtL.SetValue(Standard_False,Ucv,Vv);
2318           }
2319           else {
2320             Standard_Real U2 = 0., V2 = 0.;
2321             aV.ParametersOnS2(U2,V2); // S2 - prm
2322             VrtL.SetValue(Standard_False,U2,V2);
2323             if(LConjugated == 0)
2324               VrtL.SetValue(Standard_True,Uv,Vv);
2325             else
2326               VrtL.SetValue(Standard_True,Ucv,Vv);
2327           }
2328           Standard_Real Dist3D = VrtL.Value().Distance(aPL.Value());
2329           if(Dist3D > 1.5e-7 && DistMinL > TOL2D) {
2330             AddLast = Standard_True;
2331           }
2332         }
2333       }
2334       else {
2335         // to do: analyze internal vertex
2336       }
2337     }
2338   }
2339 }
2340
2341 static Standard_Boolean AddVertices(Handle(IntSurf_LineOn2S)& Line,
2342   const IntSurf_PntOn2S&    VrtF,
2343   const Standard_Boolean    AddFirst,
2344   const IntSurf_PntOn2S&    VrtL,
2345   const Standard_Boolean    AddLast,
2346   const Standard_Real       D3DF,
2347   const Standard_Real       D3DL)
2348 {
2349   Standard_Boolean result = Standard_False;
2350   if(AddFirst) {
2351     Standard_Real DF = Line->Value(1).Value().Distance(VrtF.Value());
2352     if((D3DF*2.) > DF && DF > 1.5e-7) {
2353       Line->InsertBefore(1,VrtF);
2354       result = Standard_True;
2355     }
2356   }
2357   if(AddLast) {
2358     Standard_Real DL = Line->Value(Line->NbPoints()).Value().Distance(VrtL.Value());
2359     if((D3DL*2.) > DL && DL > 1.5e-7) {
2360       Line->Add(VrtL);
2361       result = Standard_True;
2362     }
2363   }
2364   return result;
2365 }
2366
2367
2368 static void PutIntVertices(const Handle(IntPatch_PointLine)&    Line,
2369   Handle(IntSurf_LineOn2S)& Result,
2370   Standard_Boolean          theIsReversed,
2371   Handle(IntSurf_LineOn2S)& Vertices,
2372   const Standard_Real       ArcTol)
2373 {
2374   Standard_Integer nbp = Result->NbPoints(), nbv = Vertices->NbPoints();
2375
2376   if(nbp < 3)
2377     return;
2378
2379   const Handle(IntPatch_RLine) aRLine = Handle(IntPatch_RLine)::DownCast(Line);
2380
2381   Standard_Integer ip = 0, iv = 0;
2382   gp_Pnt aPnt;
2383   IntPatch_Point thePnt;
2384   Standard_Real U1 = 0., V1 = 0., U2 = 0., V2 = 0.;
2385
2386   for(ip = 2; ip <= (nbp-1); ip++) {
2387     const IntSurf_PntOn2S& aP = Result->Value(ip);
2388     for(iv = 1; iv <= nbv; iv++) {
2389       const IntSurf_PntOn2S& aV = Vertices->Value(iv);
2390       if(aP.IsSame(aV, Precision::Confusion(), Precision::PConfusion())) {
2391         aPnt = Result->Value(ip).Value();
2392         Result->Value(ip).ParametersOnS1(U1,V1);
2393         Result->Value(ip).ParametersOnS2(U2,V2);
2394         thePnt.SetValue(aPnt,ArcTol,Standard_False);
2395         thePnt.SetParameters(U1,V1,U2,V2);
2396
2397         Standard_Real aParam = (Standard_Real)ip;
2398
2399         if(!aRLine.IsNull())
2400         {
2401           //In fact, aRLine is always on the parametric surface.
2402           //If (theIsReversed == TRUE) then (U1, V1) - point on
2403           //parametric surface, otherwise - point on quadric.
2404           const Handle(Adaptor2d_Curve2d)& anArc = aRLine->IsArcOnS1() ?
2405                                                     aRLine->ArcOnS1() :
2406                                                     aRLine->ArcOnS2();
2407
2408           const gp_Lin2d aLin(anArc->Line());
2409           gp_Pnt2d aPSurf;
2410
2411           if(theIsReversed)
2412           {
2413             aPSurf.SetCoord(U1, V1);
2414           }
2415           else
2416           {
2417             aPSurf.SetCoord(U2, V2);
2418           }
2419
2420           aParam = ElCLib::Parameter(aLin, aPSurf);
2421         }
2422
2423         thePnt.SetParameter(aParam);
2424         Line->AddVertex(thePnt);
2425       }
2426     }
2427   }
2428 }
2429
2430 static Standard_Boolean HasInternals(Handle(IntSurf_LineOn2S)& Line,
2431   Handle(IntSurf_LineOn2S)& Vertices)
2432 {
2433   Standard_Integer nbp = Line->NbPoints(), nbv = Vertices->NbPoints();
2434   Standard_Integer ip = 0, iv = 0;
2435   Standard_Boolean result = Standard_False;
2436
2437   if(nbp < 3)
2438     return result;
2439
2440   for(ip = 2; ip <= (nbp-1); ip++) {
2441     const IntSurf_PntOn2S& aP = Line->Value(ip);
2442     for(iv = 1; iv <= nbv; iv++) {
2443       const IntSurf_PntOn2S& aV = Vertices->Value(iv);
2444       if(aP.IsSame(aV, Precision::Confusion(), Precision::PConfusion())) {
2445         result = Standard_True;
2446         break;
2447       }
2448     }
2449     if(result)
2450       break;
2451   }
2452
2453   return result;
2454 }
2455 static Handle(IntPatch_WLine) MakeSplitWLine (Handle(IntPatch_WLine)&        WLine,
2456   Standard_Boolean         Tang,
2457   IntSurf_TypeTrans        Trans1,
2458   IntSurf_TypeTrans        Trans2,
2459   Standard_Real            ArcTol,
2460   Standard_Integer         ParFirst,
2461   Standard_Integer         ParLast)
2462 {
2463   Handle(IntSurf_LineOn2S) SLine = WLine->Curve();
2464   Handle(IntSurf_LineOn2S) sline = new IntSurf_LineOn2S();
2465
2466   Standard_Integer ip = 0;
2467   for(ip = ParFirst; ip <= ParLast; ip++)
2468     sline->Add(SLine->Value(ip));
2469
2470   Handle(IntPatch_WLine) wline = new IntPatch_WLine(sline,Tang,Trans1,Trans2);
2471   wline->SetCreatingWayInfo(IntPatch_WLine::IntPatch_WLImpPrm);
2472
2473   gp_Pnt aSPnt;
2474   IntPatch_Point TPntF,TPntL;
2475   Standard_Real uu1 = 0., vv1 = 0., uu2 = 0., vv2 = 0.;
2476
2477   aSPnt = sline->Value(1).Value();
2478   sline->Value(1).ParametersOnS1(uu1,vv1);
2479   sline->Value(1).ParametersOnS2(uu2,vv2);
2480   TPntF.SetValue(aSPnt,ArcTol,Standard_False);
2481   TPntF.SetParameters(uu1,vv1,uu2,vv2);
2482   TPntF.SetParameter(1.);
2483   wline->AddVertex(TPntF);
2484   wline->SetFirstPoint(1);
2485
2486   aSPnt =  sline->Value(sline->NbPoints()).Value();
2487   sline->Value(sline->NbPoints()).ParametersOnS1(uu1,vv1);
2488   sline->Value(sline->NbPoints()).ParametersOnS2(uu2,vv2);
2489   TPntL.SetValue(aSPnt,ArcTol,Standard_False);
2490   TPntL.SetParameters(uu1,vv1,uu2,vv2);
2491   TPntL.SetParameter((Standard_Real)sline->NbPoints());
2492   wline->AddVertex(TPntL);
2493   wline->SetLastPoint(wline->NbVertex());
2494
2495   return wline;
2496 }
2497
2498 static Standard_Boolean SplitOnSegments(Handle(IntPatch_WLine)&        WLine,
2499   Standard_Boolean         Tang,
2500   IntSurf_TypeTrans        Trans1,
2501   IntSurf_TypeTrans        Trans2,
2502   Standard_Real            ArcTol,
2503   IntPatch_SequenceOfLine& Segments)
2504 {
2505   Standard_Boolean result = Standard_False;
2506   Segments.Clear();
2507
2508   Standard_Integer nbv = WLine->NbVertex();
2509   if(nbv > 3) {
2510     Standard_Integer iv = 0;
2511     for(iv = 1; iv < nbv; iv++) {
2512       Standard_Integer firstPar =
2513                     (Standard_Integer) WLine->Vertex(iv).ParameterOnLine();
2514       Standard_Integer lastPar  =
2515                     (Standard_Integer) WLine->Vertex(iv+1).ParameterOnLine();
2516       if((lastPar - firstPar) <= 1)
2517         continue;
2518       else {
2519         Handle(IntPatch_WLine) splitwline = MakeSplitWLine(WLine,Tang,Trans1,Trans2,
2520                                                                 ArcTol,firstPar,lastPar);
2521         Segments.Append(splitwline);
2522         if(!result)
2523           result = Standard_True;
2524       }
2525     }
2526   }
2527   return result;
2528 }
2529
2530 //=======================================================================
2531 //function : IsPointOnBoundary
2532 //purpose  : Returns TRUE if point <theParam> matches <theBoundary +/- thePeriod>
2533 //            with given tolerance criterion.
2534 //            For not-periodic case, thePeriod must be equal to 0.0.
2535 //=======================================================================
2536 static Standard_Boolean IsPointOnBoundary(const Standard_Real theToler2D,
2537                                           const Standard_Real theBoundary,
2538                                           const Standard_Real thePeriod,
2539                                           const Standard_Real theParam)
2540 {
2541   Standard_Real aDelta = Abs(theParam - theBoundary);
2542   if (thePeriod != 0.0)
2543   {
2544     aDelta = fmod(aDelta, thePeriod);
2545
2546     // 0 <= aDelta < thePeriod
2547     return ((aDelta < theToler2D) || ((thePeriod - aDelta) < theToler2D));
2548   }
2549
2550   // Here, thePeriod == 0.0, aDelta > 0.0
2551
2552   return (aDelta < theToler2D);
2553 }
2554
2555 //=======================================================================
2556 //function : DetectOfBoundaryAchievement
2557 //purpose  : Can change values of theNewLine (by adding the computed point on boundary,
2558 //            which parameter will be adjusted) and theIsOnBoundary variables.
2559 //=======================================================================
2560 static void DetectOfBoundaryAchievement(const Handle(Adaptor3d_Surface)& theQSurf, // quadric
2561                                         const Standard_Boolean theIsReversed,
2562                                         const Handle(IntSurf_LineOn2S)& theSourceLine,
2563                                         const Standard_Integer thePointIndex,
2564                                         const Standard_Real theToler2D,
2565                                         Handle(IntSurf_LineOn2S)& theNewLine,
2566                                         Standard_Boolean& theIsOnBoundary)
2567 {
2568   const Standard_Real aUPeriod = theQSurf->IsUPeriodic() ? theQSurf->UPeriod() : 0.0,
2569                       aVPeriod = theQSurf->IsVPeriodic() ? theQSurf->VPeriod() : 0.0;
2570   const Standard_Real aUf = theQSurf->FirstUParameter(),
2571                       aUl = theQSurf->LastUParameter(),
2572                       aVf = theQSurf->FirstVParameter(),
2573                       aVl = theQSurf->LastVParameter();
2574
2575   const IntSurf_PntOn2S &aPPrev = theSourceLine->Value(thePointIndex - 1),
2576                         &aPCurr = theSourceLine->Value(thePointIndex);
2577   Standard_Real aUPrev, aVPrev, aUCurr, aVCurr;
2578   if (theIsReversed)
2579   {
2580     aPPrev.ParametersOnS2(aUPrev, aVPrev);   // S2 - quadric, set U,V by Pnt3D
2581     aPCurr.ParametersOnS2(aUCurr, aVCurr);   // S2 - quadric, set U,V by Pnt3D
2582   }
2583   else
2584   {
2585     aPPrev.ParametersOnS1(aUPrev, aVPrev);    // S1 - quadric, set U,V by Pnt3D
2586     aPCurr.ParametersOnS1(aUCurr, aVCurr);    // S1 - quadric, set U,V by Pnt3D
2587   }
2588
2589   // Ignore cases when the WLine goes along the surface boundary completely.
2590
2591   if (IsPointOnBoundary(theToler2D, aUf, aUPeriod, aUCurr) &&
2592       !IsPointOnBoundary(theToler2D, aUf, aUPeriod, aUPrev))
2593   {
2594     theIsOnBoundary = Standard_True;
2595   }
2596   else if (IsPointOnBoundary(theToler2D, aUl, aUPeriod, aUCurr) &&
2597            !IsPointOnBoundary(theToler2D, aUl, aUPeriod, aUPrev))
2598   {
2599     theIsOnBoundary = Standard_True;
2600   }
2601   else if (IsPointOnBoundary(theToler2D, aVf, aVPeriod, aVCurr) &&
2602            !IsPointOnBoundary(theToler2D, aVf, aVPeriod, aVPrev))
2603   {
2604     theIsOnBoundary = Standard_True;
2605   }
2606   else if (IsPointOnBoundary(theToler2D, aVl, aVPeriod, aVCurr) &&
2607            !IsPointOnBoundary(theToler2D, aVl, aVPeriod, aVPrev))
2608   {
2609     theIsOnBoundary = Standard_True;
2610   }
2611
2612   if (theIsOnBoundary)
2613   {
2614     // Adjust, to avoid bad jumping of the WLine.
2615
2616     const Standard_Real aDu = (aUPrev - aUCurr);
2617     const Standard_Real aDv = (aVPrev - aVCurr);
2618     if (aUPeriod > 0.0 && (2.0*Abs(aDu) > aUPeriod))
2619     {
2620       aUCurr += Sign(aUPeriod, aDu);
2621     }
2622
2623     if (aVPeriod > 0.0 && (2.0*Abs(aDv) > aVPeriod))
2624     {
2625       aVCurr += Sign(aVPeriod, aDv);
2626     }
2627
2628     IntSurf_PntOn2S aPoint = aPCurr;
2629     aPoint.SetValue(!theIsReversed, aUCurr, aVCurr);
2630     theNewLine->Add(aPoint);
2631   }
2632 }
2633 //=======================================================================
2634 //function : DecomposeResult
2635 //purpose  : Split <theLine> in the places where it passes through seam edge
2636 //            or singularity (apex of cone or pole of sphere).
2637 //            This passage is detected by jump of U-parameter
2638 //            from point to point.
2639 //=======================================================================
2640 static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLine,
2641                                         const Standard_Boolean       IsReversed,
2642                                         const IntSurf_Quadric&       theQuad,
2643                                         const Handle(Adaptor3d_TopolTool)& thePDomain,
2644                                         const Handle(Adaptor3d_Surface)&  theQSurf, //quadric
2645                                         const Handle(Adaptor3d_Surface)&  thePSurf, //parametric
2646                                         const Standard_Real                theArcTol,
2647                                         const Standard_Real                theTolTang,
2648                                         IntPatch_SequenceOfLine&           theLines)
2649 {
2650   if(theLine->ArcType() == IntPatch_Restriction)
2651   {
2652     const Handle(IntPatch_RLine)& aRL = Handle(IntPatch_RLine)::DownCast(theLine);
2653     if(!aRL.IsNull())
2654     {
2655       const Handle(Adaptor2d_Curve2d)& anArc = aRL->IsArcOnS1() ?
2656                                         aRL->ArcOnS1() :
2657                                         aRL->ArcOnS2();
2658       if(anArc->GetType() != GeomAbs_Line)
2659       {
2660         //Restriction line must be isoline.
2661         //Other cases are not supported by
2662         //existing algorithms.
2663
2664         return Standard_False;
2665       }
2666     }
2667   }
2668
2669   const Standard_Real aDeltaUmax = M_PI_2;
2670   const Standard_Real aTOL3D = 1.e-10,
2671                       aTOL2D = Precision::PConfusion(),
2672                       aTOL2DS = Precision::PConfusion();
2673
2674   const Handle(IntSurf_LineOn2S)& aSLine = theLine->Curve();
2675
2676   if(aSLine->NbPoints() <= 2)
2677   {
2678     return Standard_False;
2679   }
2680
2681   //Deletes repeated vertices
2682   Handle(IntSurf_LineOn2S) aVLine = GetVertices(theLine,aTOL3D,aTOL2D);
2683
2684   Handle(IntSurf_LineOn2S) aSSLine(aSLine);
2685
2686   if(aSSLine->NbPoints() <= 1)
2687     return Standard_False;
2688
2689   AdjustLine(aSSLine,IsReversed,theQSurf,aTOL2D);
2690
2691   if(theLine->ArcType() == IntPatch_Walking)
2692   {
2693     Standard_Boolean isInserted = Standard_True;
2694     while(isInserted)
2695     {
2696       const Standard_Integer aNbPnts = aSSLine->NbPoints();
2697       TColStd_Array1OfInteger aPTypes(1,aNbPnts);
2698       SearchVertices(aSSLine,aVLine,aPTypes);
2699       isInserted = InsertSeamVertices(aSSLine,IsReversed,aVLine,aPTypes,aTOL2D);
2700     }
2701   }
2702
2703   const Standard_Integer aLindex = aSSLine->NbPoints();
2704   Standard_Integer aFindex = 1, aBindex = 0;
2705
2706   // build WLine parts (if any)
2707   Standard_Boolean flNextLine = Standard_True;
2708   Standard_Boolean hasBeenDecomposed = Standard_False;
2709   IntPatch_SpecPntType aPrePointExist = IntPatch_SPntNone;
2710
2711   IntSurf_PntOn2S PrePoint;
2712   while(flNextLine)
2713   {
2714     // reset variables
2715     flNextLine = Standard_False;
2716     Standard_Boolean isDecomposited = Standard_False;
2717
2718     Handle(IntSurf_LineOn2S) sline = new IntSurf_LineOn2S();
2719
2720     //if((Lindex-Findex+1) <= 2 )
2721     if((aLindex <= aFindex) && !aPrePointExist)
2722     {
2723       //break of "while(flNextLine)" cycle
2724       break;
2725     }
2726
2727     if(aPrePointExist)
2728     {
2729       const IntSurf_PntOn2S& aRefPt = aSSLine->Value(aFindex);
2730
2731       const Standard_Real aURes = theQSurf->UResolution(theArcTol),
2732                           aVRes = theQSurf->VResolution(theArcTol);
2733
2734       const Standard_Real aTol2d = (aPrePointExist == IntPatch_SPntPole) ? -1.0 :
2735               (aPrePointExist == IntPatch_SPntSeamV)? aVRes :
2736               (aPrePointExist == IntPatch_SPntSeamUV)? Max(aURes, aVRes) : aURes;
2737
2738       if(IntPatch_SpecialPoints::ContinueAfterSpecialPoint(theQSurf, thePSurf, aRefPt,
2739                                                               aPrePointExist, aTol2d,
2740                                                               PrePoint, IsReversed))
2741       {
2742         sline->Add(PrePoint);
2743
2744         //Avoid adding duplicate points.
2745         for (;aFindex <= aLindex; aFindex++)
2746         {
2747           if (!PrePoint.IsSame(aSSLine->Value(aFindex), theTolTang))
2748           {
2749             break;
2750           }
2751         }
2752       }
2753       else
2754       {
2755         //break of "while(flNextLine)" cycle
2756         break;
2757       }
2758     }
2759
2760     aPrePointExist = IntPatch_SPntNone;
2761
2762     // analyze other points
2763     for(Standard_Integer k = aFindex; k <= aLindex; k++)
2764     {
2765       if( k == aFindex )
2766       {
2767         PrePoint = aSSLine->Value(k);
2768         sline->Add(PrePoint);
2769         continue;
2770       }
2771
2772       //Check whether the current point is on the boundary of theQSurf.
2773       //If that is TRUE then the Walking-line will be decomposed in this point.
2774       //However, this boundary is not singular-point (like seam or pole of sphere).
2775       //Therefore, its processing will be simplified.
2776       Standard_Boolean isOnBoundary = Standard_False;
2777
2778       // Values of sline and isOnBoundary can be changed by this function
2779       DetectOfBoundaryAchievement(theQSurf, IsReversed, aSSLine,
2780                                   k, aTOL2D, sline, isOnBoundary);
2781
2782       aPrePointExist = IsSeamOrPole(theQSurf, aSSLine, IsReversed,
2783                                     k - 1, theTolTang, aDeltaUmax);
2784
2785       if (isOnBoundary && (aPrePointExist != IntPatch_SPntPoleSeamU))
2786       {
2787         // If the considered point is on seam then its UV-parameters
2788         // are defined to within the surface period. Therefore, we can
2789         // trust already computed parameters of this point.
2790         // However, if this point (which is on the surface boundary) is
2791         // a sphere pole or cone apex then its (point's) parameters
2792         // have to be recomputed in the code below
2793         // (see IntPatch_SpecialPoints::AddSingularPole() method).
2794         // E.g. see "bugs modalg_6 bug26684_2" test case.
2795
2796         aPrePointExist = IntPatch_SPntNone;
2797       }
2798
2799       if (aPrePointExist != IntPatch_SPntNone)
2800       {
2801         aBindex = k;
2802         isDecomposited = Standard_True;
2803         ////
2804         const IntSurf_PntOn2S& aRefPt = aSSLine->Value(aBindex-1);
2805
2806         Standard_Real aCompareTol3D = Precision::Confusion();
2807         Standard_Real aCompareTol2D = Precision::PConfusion();
2808
2809         IntSurf_PntOn2S aNewPoint = aRefPt;
2810         IntPatch_SpecPntType aLastType = IntPatch_SPntNone;
2811
2812         if(aPrePointExist == IntPatch_SPntSeamUV)
2813         {
2814           aPrePointExist = IntPatch_SPntNone;
2815           aLastType = IntPatch_SPntSeamUV;
2816           IntPatch_SpecialPoints::AddCrossUVIsoPoint(theQSurf, thePSurf,
2817                                                         aRefPt, theTolTang,
2818                                                         aNewPoint, IsReversed);
2819         }
2820         else if(aPrePointExist == IntPatch_SPntSeamV)
2821         {//WLine goes through seam
2822           aPrePointExist = IntPatch_SPntNone;
2823           aLastType = IntPatch_SPntSeamV;
2824
2825           //Not quadric point
2826           Standard_Real aU0 = 0.0, aV0 = 0.0;
2827           //Quadric point
2828           Standard_Real aUQuadRef = 0.0, aVQuadRef = 0.0;
2829
2830           if(IsReversed)
2831           {
2832             aRefPt.Parameters(aU0, aV0, aUQuadRef, aVQuadRef);
2833           }
2834           else
2835           {
2836             aRefPt.Parameters(aUQuadRef, aVQuadRef, aU0, aV0);
2837           }
2838
2839           math_Vector aTol(1, 3), aStartPoint(1,3),
2840             anInfBound(1, 3), aSupBound(1, 3);
2841
2842           //Parameters on parametric surface
2843           Standard_Real aUp = 0.0, aVp = 0.0, aUq = 0.0, aVq = 0.0;
2844           if(IsReversed)
2845           {
2846             aSSLine->Value(k).Parameters(aUp, aVp, aUq, aVq);
2847           }
2848           else
2849           {
2850             aSSLine->Value(k).Parameters(aUq, aVq, aUp, aVp);
2851           }
2852
2853           aTol(1) = thePSurf->UResolution(theArcTol);
2854           aTol(2) = thePSurf->VResolution(theArcTol);
2855           aTol(3) = theQSurf->UResolution(theArcTol);
2856           aStartPoint(1) = 0.5*(aU0 + aUp);
2857           aStartPoint(2) = 0.5*(aV0 + aVp);
2858           aStartPoint(3) = 0.5*(aUQuadRef + aUq);
2859           anInfBound(1) = thePSurf->FirstUParameter();
2860           anInfBound(2) = thePSurf->FirstVParameter();
2861           anInfBound(3) = theQSurf->FirstUParameter();
2862           aSupBound(1) = thePSurf->LastUParameter();
2863           aSupBound(2) = thePSurf->LastVParameter();
2864           aSupBound(3) = theQSurf->LastUParameter();
2865
2866           IntPatch_SpecialPoints::
2867                       AddPointOnUorVIso(theQSurf, thePSurf, aRefPt, Standard_False, 0.0,
2868                                         aTol, aStartPoint, anInfBound, aSupBound,
2869                                         aNewPoint, IsReversed);
2870         }
2871         else if(aPrePointExist == IntPatch_SPntPoleSeamU)
2872         {
2873           aPrePointExist = IntPatch_SPntNone;
2874
2875           IntPatch_Point aVert;
2876           aVert.SetValue(aRefPt);
2877           aVert.SetTolerance(theTolTang);
2878
2879           if(IntPatch_SpecialPoints::
2880                       AddSingularPole(theQSurf, thePSurf, aRefPt,
2881                                       aVert, aNewPoint, IsReversed))
2882           {
2883             aPrePointExist = IntPatch_SPntPole;
2884             aLastType = IntPatch_SPntPole;
2885             if (isOnBoundary)
2886             {
2887               // It is necessary to replace earlier added point on
2888               // the surface boundary with the pole. For that,
2889               // here we delete excess point. New point will be added later.
2890               isOnBoundary = Standard_False;
2891               sline->RemovePoint(sline->NbPoints());
2892             }
2893
2894             aCompareTol2D = -1.0;
2895           } //if(IntPatch_AddSpecialPoints::AddSingularPole(...))
2896           else
2897           {//Pole is not an intersection point
2898             aPrePointExist = IntPatch_SPntSeamU;
2899           }
2900         }
2901
2902         if(aPrePointExist == IntPatch_SPntSeamU)
2903         {//WLine goes through seam
2904           aPrePointExist = IntPatch_SPntNone;
2905           aLastType = IntPatch_SPntSeamU;
2906
2907           //Not quadric point
2908           Standard_Real aU0 = 0.0, aV0 = 0.0;
2909           //Quadric point
2910           Standard_Real aUQuadRef = 0.0, aVQuadRef = 0.0;
2911
2912           if(IsReversed)
2913           {
2914             aRefPt.Parameters(aU0, aV0, aUQuadRef, aVQuadRef);
2915           }
2916           else
2917           {
2918             aRefPt.Parameters(aUQuadRef, aVQuadRef, aU0, aV0);
2919           }
2920
2921           math_Vector aTol(1, 3), aStartPoint(1,3),
2922                       anInfBound(1, 3), aSupBound(1, 3);
2923
2924           //Parameters on parametric surface
2925           Standard_Real aUp = 0.0, aVp = 0.0, aUq = 0.0, aVq = 0.0;
2926           if (IsReversed)
2927           {
2928             aSSLine->Value(k).Parameters(aUp, aVp, aUq, aVq);
2929           }
2930           else
2931           {
2932             aSSLine->Value(k).Parameters(aUq, aVq, aUp, aVp);
2933           }
2934
2935           aTol(1) = thePSurf->UResolution(theArcTol);
2936           aTol(2) = thePSurf->VResolution(theArcTol);
2937           aTol(3) = theQSurf->VResolution(theArcTol);
2938           aStartPoint(1) = 0.5*(aU0 + aUp);
2939           aStartPoint(2) = 0.5*(aV0 + aVp);
2940           aStartPoint(3) = 0.5*(aVQuadRef + aVq);
2941           anInfBound(1) = thePSurf->FirstUParameter();
2942           anInfBound(2) = thePSurf->FirstVParameter();
2943           anInfBound(3) = theQSurf->FirstVParameter();
2944           aSupBound(1) = thePSurf->LastUParameter();
2945           aSupBound(2) = thePSurf->LastVParameter();
2946           aSupBound(3) = theQSurf->LastVParameter();
2947
2948           IntPatch_SpecialPoints::
2949                 AddPointOnUorVIso(theQSurf, thePSurf, aRefPt, Standard_True, 0.0, aTol,
2950                                   aStartPoint, anInfBound, aSupBound, aNewPoint,
2951                                   IsReversed);
2952         }
2953
2954         if(!aNewPoint.IsSame(aRefPt, aCompareTol3D, aCompareTol2D))
2955         {
2956           if (isOnBoundary)
2957             break;
2958
2959           sline->Add(aNewPoint);
2960           aPrePointExist = aLastType;
2961           PrePoint = aNewPoint;
2962         }
2963         else
2964         {
2965           if (isOnBoundary || (sline->NbPoints() == 1))
2966           {
2967             //FIRST point of the sline is the pole of the quadric.
2968             //Therefore, there is no point in decomposition.
2969
2970             // If the considered point is on surface boundary then
2971             // it is already marked as vertex. So, decomposition is
2972             // not required, too.
2973
2974             PrePoint = aRefPt;
2975             aPrePointExist = isOnBoundary ? IntPatch_SPntNone : aLastType;
2976           }
2977         }
2978
2979         ////
2980         break;
2981       } //if (aPrePointExist != IntPatch_SPntNone) cond.
2982
2983       PrePoint = aSSLine->Value(k);
2984
2985       if (isOnBoundary)
2986       {
2987         aBindex = k;
2988         isDecomposited = Standard_True;
2989         aPrePointExist = IntPatch_SPntNone;
2990         break;
2991       }
2992       else
2993       {
2994         sline->Add(aSSLine->Value(k));
2995       }
2996     } //for(Standard_Integer k = aFindex; k <= aLindex; k++)
2997
2998     //Creation of new line as part of existing theLine.
2999     //This part is defined by sline.
3000
3001     if(sline->NbPoints() == 1)
3002     {
3003       flNextLine = Standard_True;
3004
3005       if (aFindex < aBindex)
3006         aFindex = aBindex;
3007
3008       //Go to the next part of aSSLine
3009       //because we cannot create the line
3010       //with single point.
3011
3012       continue;
3013     }
3014
3015     IntSurf_PntOn2S aVF, aVL;
3016     Standard_Boolean addVF = Standard_False, addVL = Standard_False;
3017     VerifyVertices(sline,IsReversed,aVLine,aTOL2DS,theArcTol,
3018                                     thePDomain,aVF,addVF,aVL,addVL);
3019
3020     Standard_Boolean hasInternals = HasInternals(sline,aVLine);
3021
3022     Standard_Real D3F = 0., D3L = 0.;
3023     ToSmooth(sline,IsReversed,theQuad,Standard_True,D3F);
3024     ToSmooth(sline,IsReversed,theQuad,Standard_False,D3L);
3025
3026     //if(D3F <= 1.5e-7 && sline->NbPoints() >=3) {
3027     //  D3F = sline->Value(2).Value().Distance(sline->Value(3).Value());
3028     //}
3029     //if(D3L <= 1.5e-7 && sline->NbPoints() >=3) {
3030     //  D3L = sline->Value(sline->NbPoints()-1).Value().Distance(sline->
3031     //                                  Value(sline->NbPoints()-2).Value());
3032     //}
3033
3034     if(addVF || addVL)
3035     {
3036       Standard_Boolean isAdded = AddVertices(sline,aVF,addVF,aVL,addVL,D3F,D3L);
3037       if(isAdded)
3038       {
3039         ToSmooth(sline,IsReversed,theQuad,Standard_True,D3F);
3040         ToSmooth(sline,IsReversed,theQuad,Standard_False,D3L);
3041       }
3042     }
3043
3044     if(theLine->ArcType() == IntPatch_Walking)
3045     {
3046       IntPatch_Point aTPntF, aTPntL;
3047
3048       Handle(IntPatch_WLine) wline =
3049                           new IntPatch_WLine(sline,Standard_False,
3050                           theLine->TransitionOnS1(),theLine->TransitionOnS2());
3051       wline->SetCreatingWayInfo(IntPatch_WLine::IntPatch_WLImpPrm);
3052
3053       Standard_Real aU1 = 0.0, aV1 = 0.0, aU2 = 0.0, aV2 = 0.0;
3054       gp_Pnt aSPnt(sline->Value(1).Value());
3055       sline->Value(1).Parameters(aU1, aV1, aU2, aV2);
3056       aTPntF.SetValue(aSPnt,theArcTol,Standard_False);
3057       aTPntF.SetParameters(aU1, aV1, aU2, aV2);
3058       aTPntF.SetParameter(1.0);
3059       wline->AddVertex(aTPntF);
3060       wline->SetFirstPoint(1);
3061
3062       if(hasInternals)
3063       {
3064         PutIntVertices(wline,sline,IsReversed,aVLine,theArcTol);
3065       }
3066
3067       aSPnt =  sline->Value(sline->NbPoints()).Value();
3068       sline->Value(sline->NbPoints()).Parameters(aU1, aV1, aU2, aV2);
3069       aTPntL.SetValue(aSPnt,theArcTol,Standard_False);
3070       aTPntL.SetParameters(aU1, aV1, aU2, aV2);
3071       aTPntL.SetParameter(sline->NbPoints());
3072       wline->AddVertex(aTPntL);
3073       wline->SetLastPoint(wline->NbVertex());
3074
3075       IntPatch_SequenceOfLine segm;
3076       Standard_Boolean isSplited = SplitOnSegments(wline,Standard_False,
3077                       theLine->TransitionOnS1(),theLine->TransitionOnS2(),theArcTol,segm);
3078
3079       if(!isSplited)
3080       {
3081         theLines.Append(wline);
3082       }
3083       else
3084       {
3085         Standard_Integer nbsegms = segm.Length();
3086         Standard_Integer iseg = 0;
3087         for(iseg = 1; iseg <= nbsegms; iseg++)
3088           theLines.Append(segm(iseg));
3089       }
3090     }
3091     else
3092     {//theLine->ArcType() == IntPatch_Restriction
3093       if(!isDecomposited && !hasBeenDecomposed)
3094       {
3095         //The line has not been changed
3096         theLines.Append(Handle(IntPatch_RLine)::DownCast(theLine));
3097         return hasBeenDecomposed;
3098       }
3099
3100       IntPatch_Point aTPnt;
3101       gp_Pnt2d aPSurf;
3102       gp_Pnt aSPnt;
3103
3104       Handle(IntPatch_RLine) aRLine = new IntPatch_RLine(*Handle(IntPatch_RLine)::DownCast(theLine));
3105
3106       aRLine->ClearVertexes();
3107       aRLine->SetCurve(sline);
3108
3109       if(hasInternals)
3110       {
3111         PutIntVertices(aRLine,sline,IsReversed,aVLine,theArcTol);
3112       }
3113
3114       const Handle(Adaptor2d_Curve2d)& anArc = aRLine->IsArcOnS1() ?
3115                                                 aRLine->ArcOnS1() :
3116                                                 aRLine->ArcOnS2();
3117
3118       Standard_Real aFPar = anArc->FirstParameter(),
3119                     aLPar = anArc->LastParameter();
3120
3121       const IntSurf_PntOn2S &aRFirst = sline->Value(1),
3122                             &aRLast = sline->Value(sline->NbPoints());
3123
3124       const gp_Lin2d aLin(anArc->Line());
3125
3126       for(Standard_Integer aFLIndex = 0; aFLIndex < 2; aFLIndex++)
3127       {
3128         Standard_Real aU1 = 0.0, aV1 = 0.0, aU2 = 0.0, aV2 = 0.0;
3129         if(aFLIndex == 0)
3130         {
3131           aRFirst.Parameters(aU1, aV1, aU2, aV2);
3132           aSPnt.SetXYZ(aRFirst.Value().XYZ());
3133         }
3134         else
3135         {
3136           aRLast.Parameters(aU1, aV1, aU2, aV2);
3137           aSPnt.SetXYZ(aRLast.Value().XYZ());
3138         }
3139
3140         if(IsReversed)
3141         {
3142           aPSurf.SetCoord(aU1, aV1);
3143         }
3144         else
3145         {
3146           aPSurf.SetCoord(aU2, aV2);
3147         }
3148
3149         Standard_Real aPar = ElCLib::Parameter(aLin, aPSurf);
3150
3151         if(aFLIndex == 0)
3152         {
3153           aFPar = Max(aFPar, aPar);
3154           aPar = aFPar;
3155         }
3156         else
3157         {
3158           aLPar = Min(aLPar, aPar);
3159           aPar = aLPar;
3160         }
3161
3162         aTPnt.SetParameter(aPar);
3163         aTPnt.SetValue(aSPnt,theArcTol,Standard_False);
3164         aTPnt.SetParameters(aU1, aV1, aU2, aV2);
3165
3166         aRLine->AddVertex(aTPnt);
3167       }
3168
3169       if(aLPar - aFPar > Precision::PConfusion())
3170       {
3171         aRLine->SetFirstPoint(1);
3172         aRLine->SetLastPoint(aRLine->NbVertex());
3173
3174         anArc->Trim(aFPar, aLPar, theArcTol);
3175
3176         theLines.Append(aRLine);
3177       }
3178     }
3179
3180     if(isDecomposited)
3181     {
3182       aFindex = aBindex;
3183       flNextLine = hasBeenDecomposed = Standard_True;
3184     }
3185   }
3186
3187   return hasBeenDecomposed;
3188 }
3189
3190 //=======================================================================
3191 //function : CheckSegmSegm
3192 //purpose  : Returns TRUE if the segment [theParF, theParL] is included
3193 //            in the segment [theRefParF, theRefParL] segment.
3194 //=======================================================================
3195 static Standard_Boolean CheckSegmSegm(const Standard_Real theRefParF,
3196                                       const Standard_Real theRefParL,
3197                                       const Standard_Real theParF,
3198                                       const Standard_Real theParL)
3199 {
3200   if((theParF < theRefParF) || (theParF > theRefParL))
3201   {
3202     return Standard_False;
3203   }
3204
3205   if((theParL < theRefParF) || (theParL > theRefParL))
3206   {
3207     return Standard_False;
3208   }
3209
3210   return Standard_True;
3211 }
3212
3213 //=======================================================================
3214 //function : IsCoincide
3215 //purpose  : Check, if theLine is coincided with theArc (in 2d-space).
3216 //
3217 // Attention!!!
3218 //            Cases when theArc is not 2d-line adaptor are supported by
3219 //          TopOpeBRep classes only (i.e. are archaic).
3220 //=======================================================================
3221 Standard_Boolean IsCoincide(IntPatch_TheSurfFunction& theFunc,
3222                             const Handle(IntPatch_PointLine)& theLine,
3223                             const Handle(Adaptor2d_Curve2d)& theArc,
3224                             const Standard_Boolean isTheSurface1Using, //Surf1 is parametric?
3225                             const Standard_Real theToler3D,
3226                             const Standard_Real theToler2D,
3227                             const Standard_Real thePeriod) // Period of parametric surface in direction which is perpendicular to theArc direction.
3228 {
3229   const Standard_Real aCoeffs[] = { 0.02447174185,  0.09549150281, 0.20610737385, 0.34549150281, /*Sin(x)*Sin(x)*/
3230                                     0.5, 0.65450849719, 0.79389262615 };
3231   if(theLine->ArcType() == IntPatch_Restriction)
3232   {//Restriction-restriction processing
3233     const Handle(IntPatch_RLine)& aRL2 = Handle(IntPatch_RLine)::DownCast(theLine);
3234     const Handle(Adaptor2d_Curve2d)& anArc = aRL2->IsArcOnS1() ? aRL2->ArcOnS1() : aRL2->ArcOnS2();
3235
3236     if(anArc->GetType() != GeomAbs_Line)
3237     {
3238       //Restriction line must be isoline.
3239       //Other cases are not supported by
3240       //existing algorithms.
3241
3242       return Standard_False;
3243     }
3244
3245     const gp_Lin2d aLin1(theArc->Line()),
3246                    aLin2(anArc->Line());
3247
3248     if(!aLin1.Direction().IsParallel(aLin2.Direction(), Precision::Angular()))
3249     {
3250       return Standard_False;
3251     }
3252
3253     const Standard_Real aDist =
3254             theArc->Line().Distance(anArc->Line());
3255     if((aDist < theToler2D) || (Abs(aDist - thePeriod) < theToler2D))
3256     {
3257       const Standard_Real aRf = theArc->FirstParameter(),
3258                           aRl = theArc->LastParameter();
3259       const Standard_Real aParf = anArc->FirstParameter(),
3260                           aParl = anArc->LastParameter();
3261       const gp_Pnt2d aP1(ElCLib::Value(aParf, aLin2)),
3262                      aP2(ElCLib::Value(aParl, aLin2));
3263
3264       Standard_Real aParam1 = ElCLib::Parameter(aLin1, aP1),
3265                     aParam2 = ElCLib::Parameter(aLin1, aP2);
3266
3267       if(CheckSegmSegm(aRf, aRl, aParam1, aParam2))
3268         return Standard_True;
3269
3270       //Lines are parallel. Therefore, there is no point in
3271       //projecting points to another line in order to check
3272       //if segment second line is included in segment of first one.
3273
3274       return CheckSegmSegm(aParam1, aParam2, aRf, aRl);
3275     }
3276
3277     return Standard_False;
3278   }
3279
3280   const Standard_Integer aNbPnts = theLine->NbPnts();
3281   const Standard_Real aUAf = theArc->FirstParameter(),
3282                       aUAl = theArc->LastParameter();
3283   const gp_Lin2d anArcLin(theArc->Line());
3284
3285   math_Vector aX(1, 2), aVal(1, 1);
3286
3287   for(Standard_Integer aPtID = 1; aPtID <= aNbPnts; aPtID++)
3288   {
3289     Standard_Real aUf = 0.0, aVf = 0.0;
3290     if(isTheSurface1Using)
3291       theLine->Point(aPtID).ParametersOnS1(aUf, aVf);
3292     else
3293       theLine->Point(aPtID).ParametersOnS2(aUf, aVf);
3294
3295     //Take 2d-point in parametric surface (because theArc is
3296     //2d-line in parametric surface).
3297     const gp_Pnt2d aPloc(aUf, aVf);
3298
3299     const Standard_Real aRParam = ElCLib::Parameter(anArcLin, aPloc);
3300
3301     if((aRParam < aUAf) || (aRParam > aUAl))
3302       return Standard_False;
3303
3304     const gp_Pnt2d aPmin(ElCLib::Value(aRParam, anArcLin));
3305
3306     const Standard_Real aDist = aPloc.Distance(aPmin);
3307     if((aDist < theToler2D) || (Abs(aDist - thePeriod) < theToler2D))
3308     {//Considered point is in Restriction line.
3309      //Go to the next point.
3310       continue;
3311     }
3312
3313     //Check if intermediate points between aPloc and theArc are
3314     //intersection point (i.e. if aPloc is in tangent zone between
3315     //two intersected surfaces).
3316
3317     const Standard_Real aUl = aPmin.X(), aVl = aPmin.Y();
3318
3319     Standard_Real aU, aV;
3320     Standard_Real dU = aUl - aUf, dV = aVl - aVf;
3321     for(Standard_Integer i = 0; i < 7; i++)
3322     {
3323       aU = aUf + aCoeffs[i] * dU;
3324       aV = aVf + aCoeffs[i] * dV;
3325
3326       aX.Value(1) = aU;
3327       aX.Value(2) = aV;
3328
3329       if(!theFunc.Value(aX, aVal))
3330       {
3331         return Standard_False;
3332       }
3333
3334       if(Abs(theFunc.Root()) > theToler3D)
3335       {
3336         return Standard_False;
3337       }
3338     }
3339   }
3340
3341   return Standard_True;
3342 }