src/ProjLib/ProjLib_CompProjectedCurve.cxx

   1 // Created on: 1997-09-23
   2 // Created by: Roman BORISOV
   3 // Copyright (c) 1997-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17
  18 #include <algorithm>
  19
  20 #include <Adaptor2d_HCurve2d.hxx>
  21 #include <Adaptor3d_HCurve.hxx>
  22 #include <Adaptor3d_HSurface.hxx>
  23 #include <Extrema_ExtCS.hxx>
  24 #include <Extrema_ExtPS.hxx>
  25 #include <Extrema_GenLocateExtPS.hxx>
  26 #include <Extrema_POnCurv.hxx>
  27 #include <Extrema_POnSurf.hxx>
  28 #include <GeomAbs_CurveType.hxx>
  29 #include <GeomLib.hxx>
  30 #include <gp_Mat2d.hxx>
  31 #include <gp_Pnt2d.hxx>
  32 #include <gp_Vec2d.hxx>
  33 #include <gp_XY.hxx>
  34 #include <Precision.hxx>
  35 #include <ProjLib_CompProjectedCurve.hxx>
  36 #include <ProjLib_HCompProjectedCurve.hxx>
  37 #include <ProjLib_PrjResolve.hxx>
  38 #include <Standard_DomainError.hxx>
  39 #include <Standard_NoSuchObject.hxx>
  40 #include <Standard_NotImplemented.hxx>
  41 #include <Standard_OutOfRange.hxx>
  42 #include <TColgp_HSequenceOfPnt.hxx>
  43 #include <Adaptor3d_CurveOnSurface.hxx>
  44 #include <Geom2d_Line.hxx>
  45 #include <Geom2dAdaptor_HCurve.hxx>
  46 #include <Extrema_ExtCC.hxx>
  47 #include <NCollection_Vector.hxx>
  48
  49 #define FuncTol 1.e-10
  50
  51 #ifdef OCCT_DEBUG_CHRONO
  52 #include <OSD_Timer.hxx>
  53
  54 static OSD_Chronometer chr_init_point, chr_dicho_bound;
  55
  56 Standard_EXPORT Standard_Real t_init_point, t_dicho_bound;
  57 Standard_EXPORT Standard_Integer init_point_count, dicho_bound_count;
  58
  59 static void InitChron(OSD_Chronometer& ch)
  60 {
  61   ch.Reset();
  62   ch.Start();
  63 }
  64
  65 static void ResultChron( OSD_Chronometer & ch, Standard_Real & time)
  66 {
  67   Standard_Real tch ;
  68   ch.Stop();
  69   ch.Show(tch);
  70   time=time +tch;
  71 }
  72 #endif
  73
  74 // Structure to perform splits computation.
  75 // This structure is not thread-safe since operations under mySplits should be performed in a critical section.
  76 // myPeriodicDir - 0 for U periodicity and 1 for V periodicity.
  77 struct SplitDS
  78 {
  79   SplitDS(const Handle(Adaptor3d_HCurve)   &theCurve,
  80           const Handle(Adaptor3d_HSurface) &theSurface,
  81           NCollection_Vector<Standard_Real> &theSplits)
  82   : myCurve(theCurve),
  83     mySurface(theSurface),
  84     mySplits(theSplits)
  85   { }
  86
  87   // Assignment operator is forbidden.
  88   void operator=(const SplitDS &theSplitDS);
  89
  90   const Handle(Adaptor3d_HCurve) myCurve;
  91   const Handle(Adaptor3d_HSurface) mySurface;
  92   NCollection_Vector<Standard_Real> &mySplits;
  93
  94   Standard_Real myPerMinParam;
  95   Standard_Real myPerMaxParam;
  96   Standard_Integer myPeriodicDir;
  97
  98   Extrema_ExtCC *myExtCC;
  99   Extrema_ExtPS *myExtPS;
 100 };
 101
 102   //! Compute split points in the parameter space of the curve.
 103   static void BuildCurveSplits(const Handle(Adaptor3d_HCurve)   &theCurve,
 104                                const Handle(Adaptor3d_HSurface) &theSurface,
 105                                const Standard_Real theTolU,
 106                                const Standard_Real theTolV,
 107                                NCollection_Vector<Standard_Real> &theSplits);
 108
 109   //! Perform splitting on a specified direction. Sub-method in BuildCurveSplits.
 110   static void SplitOnDirection(SplitDS & theSplitDS);
 111
 112   //! Perform recursive search of the split points.
 113   static void FindSplitPoint(SplitDS & theSplitDS,
 114                              const Standard_Real theMinParam,
 115                              const Standard_Real theMaxParam);
 116
 117
 118 //=======================================================================
 119 //function : Comparator
 120 //purpose  : used in sort algorithm
 121 //=======================================================================
 122 inline Standard_Boolean Comparator(const Standard_Real theA,
 123                                    const Standard_Real theB)
 124 {
 125   return theA < theB;
 126 }
 127
 128 //=======================================================================
 129 //function : d1
 130 //purpose  : computes first derivative of the projected curve
 131 //=======================================================================
 132
 133 static void d1(const Standard_Real t,
 134   const Standard_Real u,
 135   const Standard_Real v,
 136   gp_Vec2d& V,
 137   const Handle(Adaptor3d_HCurve)& Curve,
 138   const Handle(Adaptor3d_HSurface)& Surface)
 139 {
 140   gp_Pnt S, C;
 141   gp_Vec DS1_u, DS1_v, DS2_u, DS2_uv, DS2_v, DC1_t;
 142   Surface->D2(u, v, S, DS1_u, DS1_v, DS2_u, DS2_v, DS2_uv);
 143   Curve->D1(t, C, DC1_t);
 144   gp_Vec Ort(C, S);// Ort = S - C
 145
 146   gp_Vec2d dE_dt(-DC1_t*DS1_u, -DC1_t*DS1_v);
 147   gp_XY dE_du(DS1_u*DS1_u + Ort*DS2_u,
 148     DS1_u*DS1_v + Ort*DS2_uv);
 149   gp_XY dE_dv(DS1_v*DS1_u + Ort*DS2_uv,
 150     DS1_v*DS1_v + Ort*DS2_v);
 151
 152   Standard_Real det = dE_du.X()*dE_dv.Y() - dE_du.Y()*dE_dv.X();
 153   if (fabs(det) < gp::Resolution()) throw Standard_ConstructionError();
 154
 155   gp_Mat2d M(gp_XY(dE_dv.Y()/det, -dE_du.Y()/det),
 156     gp_XY(-dE_dv.X()/det, dE_du.X()/det));
 157
 158   V = - gp_Vec2d(gp_Vec2d(M.Row(1))*dE_dt, gp_Vec2d(M.Row(2))*dE_dt);
 159 }
 160
 161 //=======================================================================
 162 //function : d2
 163 //purpose  : computes second derivative of the projected curve
 164 //=======================================================================
 165
 166 static void d2(const Standard_Real t,
 167   const Standard_Real u,
 168   const Standard_Real v,
 169   gp_Vec2d& V1, gp_Vec2d& V2,
 170   const Handle(Adaptor3d_HCurve)& Curve,
 171   const Handle(Adaptor3d_HSurface)& Surface)
 172 {
 173   gp_Pnt S, C;
 174   gp_Vec DS1_u, DS1_v, DS2_u, DS2_uv, DS2_v,
 175     DS3_u, DS3_v, DS3_uuv, DS3_uvv,
 176     DC1_t, DC2_t;
 177   Surface->D3(u, v, S, DS1_u, DS1_v, DS2_u, DS2_v, DS2_uv,
 178     DS3_u, DS3_v, DS3_uuv, DS3_uvv);
 179   Curve->D2(t, C, DC1_t, DC2_t);
 180   gp_Vec Ort(C, S);
 181
 182   gp_Vec2d dE_dt(-DC1_t*DS1_u, -DC1_t*DS1_v);
 183   gp_XY dE_du(DS1_u*DS1_u + Ort*DS2_u,
 184     DS1_u*DS1_v + Ort*DS2_uv);
 185   gp_XY dE_dv(DS1_v*DS1_u + Ort*DS2_uv,
 186     DS1_v*DS1_v + Ort*DS2_v);
 187
 188   Standard_Real det = dE_du.X()*dE_dv.Y() - dE_du.Y()*dE_dv.X();
 189   if (fabs(det) < gp::Resolution()) throw Standard_ConstructionError();
 190
 191   gp_Mat2d M(gp_XY(dE_dv.Y()/det, -dE_du.Y()/det),
 192     gp_XY(-dE_dv.X()/det, dE_du.X()/det));
 193
 194   // First derivative
 195   V1 = - gp_Vec2d(gp_Vec2d(M.Row(1))*dE_dt, gp_Vec2d(M.Row(2))*dE_dt);
 196
 197   /* Second derivative */
 198
 199   // Computation of d2E_dt2 = S1
 200   gp_Vec2d d2E_dt(-DC2_t*DS1_u, -DC2_t*DS1_v);
 201
 202   // Computation of 2*(d2E/dtdX)(dX/dt) = S2
 203   gp_Vec2d d2E1_dtdX(-DC1_t*DS2_u,
 204     -DC1_t*DS2_uv);
 205   gp_Vec2d d2E2_dtdX(-DC1_t*DS2_uv,
 206     -DC1_t*DS2_v);
 207   gp_Vec2d S2 = 2*gp_Vec2d(d2E1_dtdX*V1, d2E2_dtdX*V1);
 208
 209   // Computation of (d2E/dX2)*(dX/dt)2 = S3
 210
 211   // Row11 = (d2E1/du2, d2E1/dudv)
 212   Standard_Real tmp;
 213   gp_Vec2d Row11(3*DS1_u*DS2_u + Ort*DS3_u,
 214     tmp = 2*DS1_u*DS2_uv +
 215     DS1_v*DS2_u + Ort*DS3_uuv);
 216
 217   // Row12 = (d2E1/dudv, d2E1/dv2)
 218   gp_Vec2d Row12(tmp, DS2_v*DS1_u + 2*DS1_v*DS2_uv +
 219     Ort*DS3_uvv);
 220
 221   // Row21 = (d2E2/du2, d2E2/dudv)
 222   gp_Vec2d Row21(DS2_u*DS1_v + 2*DS1_u*DS2_uv + Ort*DS3_uuv,
 223     tmp = 2*DS2_uv*DS1_v + DS1_u*DS2_v + Ort*DS3_uvv);
 224
 225   // Row22 = (d2E2/duv, d2E2/dvdv)
 226   gp_Vec2d Row22(tmp, 3*DS1_v*DS2_v + Ort*DS3_v);
 227
 228   gp_Vec2d S3(V1*gp_Vec2d(Row11*V1, Row12*V1),
 229     V1*gp_Vec2d(Row21*V1, Row22*V1));
 230
 231   gp_Vec2d Sum = d2E_dt + S2 + S3;
 232
 233   V2 = - gp_Vec2d(gp_Vec2d(M.Row(1))*Sum, gp_Vec2d(M.Row(2))*Sum);
 234 }
 235 //=======================================================================
 236 //function : d1CurveOnSurf
 237 //purpose  : computes first derivative of the 3d projected curve
 238 //=======================================================================
 239
 240 #if 0
 241 static void d1CurvOnSurf(const Standard_Real t,
 242   const Standard_Real u,
 243   const Standard_Real v,
 244   gp_Vec& V,
 245   const Handle(Adaptor3d_HCurve)& Curve,
 246   const Handle(Adaptor3d_HSurface)& Surface)
 247 {
 248   gp_Pnt S, C;
 249   gp_Vec2d V2d;
 250   gp_Vec DS1_u, DS1_v, DS2_u, DS2_uv, DS2_v, DC1_t;
 251   Surface->D2(u, v, S, DS1_u, DS1_v, DS2_u, DS2_v, DS2_uv);
 252   Curve->D1(t, C, DC1_t);
 253   gp_Vec Ort(C, S);// Ort = S - C
 254
 255   gp_Vec2d dE_dt(-DC1_t*DS1_u, -DC1_t*DS1_v);
 256   gp_XY dE_du(DS1_u*DS1_u + Ort*DS2_u,
 257     DS1_u*DS1_v + Ort*DS2_uv);
 258   gp_XY dE_dv(DS1_v*DS1_u + Ort*DS2_uv,
 259     DS1_v*DS1_v + Ort*DS2_v);
 260
 261   Standard_Real det = dE_du.X()*dE_dv.Y() - dE_du.Y()*dE_dv.X();
 262   if (fabs(det) < gp::Resolution()) throw Standard_ConstructionError();
 263
 264   gp_Mat2d M(gp_XY(dE_dv.Y()/det, -dE_du.Y()/det),
 265     gp_XY(-dE_dv.X()/det, dE_du.X()/det));
 266
 267   V2d = - gp_Vec2d(gp_Vec2d(M.Row(1))*dE_dt, gp_Vec2d(M.Row(2))*dE_dt);
 268
 269   V = DS1_u * V2d.X() + DS1_v * V2d.Y();
 270
 271 }
 272 #endif
 273
 274 //=======================================================================
 275 //function : d2CurveOnSurf
 276 //purpose  : computes second derivative of the 3D projected curve
 277 //=======================================================================
 278
 279 static void d2CurvOnSurf(const Standard_Real t,
 280   const Standard_Real u,
 281   const Standard_Real v,
 282   gp_Vec& V1 , gp_Vec& V2 ,
 283   const Handle(Adaptor3d_HCurve)& Curve,
 284   const Handle(Adaptor3d_HSurface)& Surface)
 285 {
 286   gp_Pnt S, C;
 287   gp_Vec2d V12d,V22d;
 288   gp_Vec DS1_u, DS1_v, DS2_u, DS2_uv, DS2_v,
 289     DS3_u, DS3_v, DS3_uuv, DS3_uvv,
 290     DC1_t, DC2_t;
 291   Surface->D3(u, v, S, DS1_u, DS1_v, DS2_u, DS2_v, DS2_uv,
 292     DS3_u, DS3_v, DS3_uuv, DS3_uvv);
 293   Curve->D2(t, C, DC1_t, DC2_t);
 294   gp_Vec Ort(C, S);
 295
 296   gp_Vec2d dE_dt(-DC1_t*DS1_u, -DC1_t*DS1_v);
 297   gp_XY dE_du(DS1_u*DS1_u + Ort*DS2_u,
 298     DS1_u*DS1_v + Ort*DS2_uv);
 299   gp_XY dE_dv(DS1_v*DS1_u + Ort*DS2_uv,
 300     DS1_v*DS1_v + Ort*DS2_v);
 301
 302   Standard_Real det = dE_du.X()*dE_dv.Y() - dE_du.Y()*dE_dv.X();
 303   if (fabs(det) < gp::Resolution()) throw Standard_ConstructionError();
 304
 305   gp_Mat2d M(gp_XY(dE_dv.Y()/det, -dE_du.Y()/det),
 306     gp_XY(-dE_dv.X()/det, dE_du.X()/det));
 307
 308   // First derivative
 309   V12d = - gp_Vec2d(gp_Vec2d(M.Row(1))*dE_dt, gp_Vec2d(M.Row(2))*dE_dt);
 310
 311   /* Second derivative */
 312
 313   // Computation of d2E_dt2 = S1
 314   gp_Vec2d d2E_dt(-DC2_t*DS1_u, -DC2_t*DS1_v);
 315
 316   // Computation of 2*(d2E/dtdX)(dX/dt) = S2
 317   gp_Vec2d d2E1_dtdX(-DC1_t*DS2_u,
 318     -DC1_t*DS2_uv);
 319   gp_Vec2d d2E2_dtdX(-DC1_t*DS2_uv,
 320     -DC1_t*DS2_v);
 321   gp_Vec2d S2 = 2*gp_Vec2d(d2E1_dtdX*V12d, d2E2_dtdX*V12d);
 322
 323   // Computation of (d2E/dX2)*(dX/dt)2 = S3
 324
 325   // Row11 = (d2E1/du2, d2E1/dudv)
 326   Standard_Real tmp;
 327   gp_Vec2d Row11(3*DS1_u*DS2_u + Ort*DS3_u,
 328     tmp = 2*DS1_u*DS2_uv +
 329     DS1_v*DS2_u + Ort*DS3_uuv);
 330
 331   // Row12 = (d2E1/dudv, d2E1/dv2)
 332   gp_Vec2d Row12(tmp, DS2_v*DS1_u + 2*DS1_v*DS2_uv +
 333     Ort*DS3_uvv);
 334
 335   // Row21 = (d2E2/du2, d2E2/dudv)
 336   gp_Vec2d Row21(DS2_u*DS1_v + 2*DS1_u*DS2_uv + Ort*DS3_uuv,
 337     tmp = 2*DS2_uv*DS1_v + DS1_u*DS2_v + Ort*DS3_uvv);
 338
 339   // Row22 = (d2E2/duv, d2E2/dvdv)
 340   gp_Vec2d Row22(tmp, 3*DS1_v*DS2_v + Ort*DS3_v);
 341
 342   gp_Vec2d S3(V12d*gp_Vec2d(Row11*V12d, Row12*V12d),
 343     V12d*gp_Vec2d(Row21*V12d, Row22*V12d));
 344
 345   gp_Vec2d Sum = d2E_dt + S2 + S3;
 346
 347   V22d = - gp_Vec2d(gp_Vec2d(M.Row(1))*Sum, gp_Vec2d(M.Row(2))*Sum);
 348
 349   V1 = DS1_u * V12d.X() + DS1_v * V12d.Y();
 350   V2 =     DS2_u * V12d.X() *V12d.X()
 351     +     DS1_u * V22d.X()
 352     + 2 * DS2_uv * V12d.X() *V12d.Y()
 353     +     DS2_v * V12d.Y() * V12d.Y()
 354     +     DS1_v * V22d.Y();
 355 }
 356
 357 //=======================================================================
 358 //function : ExactBound
 359 //purpose  : computes exact boundary point
 360 //=======================================================================
 361
 362 static Standard_Boolean ExactBound(gp_Pnt& Sol,
 363   const Standard_Real NotSol,
 364   const Standard_Real Tol,
 365   const Standard_Real TolU,
 366   const Standard_Real TolV,
 367   const Handle(Adaptor3d_HCurve)& Curve,
 368   const Handle(Adaptor3d_HSurface)& Surface)
 369 {
 370   Standard_Real U0, V0, t, t1, t2, FirstU, LastU, FirstV, LastV;
 371   gp_Pnt2d POnS;
 372   U0 = Sol.Y();
 373   V0 = Sol.Z();
 374   FirstU = Surface->FirstUParameter();
 375   LastU = Surface->LastUParameter();
 376   FirstV = Surface->FirstVParameter();
 377   LastV = Surface->LastVParameter();
 378   // Here we have to compute the boundary that projection is going to intersect
 379   gp_Vec2d D2d;
 380   //these variables are to estimate which boundary has more apportunity
 381   //to be intersected
 382   Standard_Real RU1, RU2, RV1, RV2;
 383   d1(Sol.X(), U0, V0, D2d, Curve, Surface);
 384   // Here we assume that D2d != (0, 0)
 385   if(Abs(D2d.X()) < gp::Resolution())
 386   {
 387     RU1 = Precision::Infinite();
 388     RU2 = Precision::Infinite();
 389     RV1 = V0 - FirstV;
 390     RV2 = LastV - V0;
 391   }
 392   else if(Abs(D2d.Y()) < gp::Resolution())
 393   {
 394     RU1 = U0 - FirstU;
 395     RU2 = LastU - U0;
 396     RV1 = Precision::Infinite();
 397     RV2 = Precision::Infinite();
 398   }
 399   else
 400   {
 401     RU1 = gp_Pnt2d(U0, V0).
 402       Distance(gp_Pnt2d(FirstU, V0 + (FirstU - U0)*D2d.Y()/D2d.X()));
 403     RU2 = gp_Pnt2d(U0, V0).
 404       Distance(gp_Pnt2d(LastU, V0 + (LastU - U0)*D2d.Y()/D2d.X()));
 405     RV1 = gp_Pnt2d(U0, V0).
 406       Distance(gp_Pnt2d(U0 + (FirstV - V0)*D2d.X()/D2d.Y(), FirstV));
 407     RV2 = gp_Pnt2d(U0, V0).
 408       Distance(gp_Pnt2d(U0 + (LastV - V0)*D2d.X()/D2d.Y(), LastV));
 409   }
 410   TColgp_SequenceOfPnt Seq;
 411   Seq.Append(gp_Pnt(FirstU, RU1, 2));
 412   Seq.Append(gp_Pnt(LastU, RU2, 2));
 413   Seq.Append(gp_Pnt(FirstV, RV1, 3));
 414   Seq.Append(gp_Pnt(LastV, RV2, 3));
 415   Standard_Integer i, j;
 416   for(i = 1; i <= 3; i++)
 417   {
 418     for(j = 1; j <= 4-i; j++)
 419     {
 420       if(Seq(j).Y() < Seq(j+1).Y())
 421       {
 422         gp_Pnt swp;
 423         swp = Seq.Value(j+1);
 424         Seq.ChangeValue(j+1) = Seq.Value(j);
 425         Seq.ChangeValue(j) = swp;
 426       }
 427     }
 428   }
 429
 430   t = Sol.X ();
 431   t1 = Min (Sol.X (), NotSol);
 432   t2 = Max (Sol.X (), NotSol);
 433
 434   Standard_Boolean isDone = Standard_False;
 435   while (!Seq.IsEmpty ())
 436   {
 437     gp_Pnt P;
 438     P = Seq.Last ();
 439     Seq.Remove (Seq.Length ());
 440     ProjLib_PrjResolve aPrjPS (Curve->Curve (),
 441       Surface->Surface (),
 442       Standard_Integer (P.Z ()));
 443     if (Standard_Integer (P.Z ()) == 2)
 444     {
 445       aPrjPS.Perform (t, P.X (), V0, gp_Pnt2d (Tol, TolV),
 446         gp_Pnt2d (t1, Surface->FirstVParameter ()),
 447         gp_Pnt2d (t2, Surface->LastVParameter ()), FuncTol);
 448       if (!aPrjPS.IsDone ()) continue;
 449       POnS = aPrjPS.Solution ();
 450       Sol = gp_Pnt (POnS.X (), P.X (), POnS.Y ());
 451       isDone = Standard_True;
 452       break;
 453     }
 454     else
 455     {
 456       aPrjPS.Perform (t, U0, P.X (), gp_Pnt2d (Tol, TolU),
 457         gp_Pnt2d (t1, Surface->FirstUParameter ()),
 458         gp_Pnt2d (t2, Surface->LastUParameter ()), FuncTol);
 459       if (!aPrjPS.IsDone ()) continue;
 460       POnS = aPrjPS.Solution ();
 461       Sol = gp_Pnt (POnS.X (), POnS.Y (), P.X ());
 462       isDone = Standard_True;
 463       break;
 464     }
 465   }
 466
 467   return isDone;
 468 }
 469
 470 //=======================================================================
 471 //function : DichExactBound
 472 //purpose  : computes exact boundary point
 473 //=======================================================================
 474
 475 static void DichExactBound(gp_Pnt& Sol,
 476   const Standard_Real NotSol,
 477   const Standard_Real Tol,
 478   const Standard_Real TolU,
 479   const Standard_Real TolV,
 480   const Handle(Adaptor3d_HCurve)& Curve,
 481   const Handle(Adaptor3d_HSurface)& Surface)
 482 {
 483 #ifdef OCCT_DEBUG_CHRONO
 484   InitChron(chr_dicho_bound);
 485 #endif
 486
 487   Standard_Real U0, V0, t;
 488   gp_Pnt2d POnS;
 489   U0 = Sol.Y();
 490   V0 = Sol.Z();
 491   ProjLib_PrjResolve aPrjPS(Curve->Curve(), Surface->Surface(), 1);
 492
 493   Standard_Real aNotSol = NotSol;
 494   while (fabs(Sol.X() - aNotSol) > Tol)
 495   {
 496     t = (Sol.X() + aNotSol)/2;
 497     aPrjPS.Perform(t, U0, V0, gp_Pnt2d(TolU, TolV),
 498       gp_Pnt2d(Surface->FirstUParameter(),Surface->FirstVParameter()),
 499       gp_Pnt2d(Surface->LastUParameter(),Surface->LastVParameter()),
 500       FuncTol, Standard_True);
 501
 502     if (aPrjPS.IsDone())
 503     {
 504       POnS = aPrjPS.Solution();
 505       Sol = gp_Pnt(t, POnS.X(), POnS.Y());
 506       U0=Sol.Y();
 507       V0=Sol.Z();
 508     }
 509     else aNotSol = t;
 510   }
 511 #ifdef OCCT_DEBUG_CHRONO
 512   ResultChron(chr_dicho_bound,t_dicho_bound);
 513   dicho_bound_count++;
 514 #endif
 515 }
 516
 517 //=======================================================================
 518 //function : InitialPoint
 519 //purpose  :
 520 //=======================================================================
 521
 522 static Standard_Boolean InitialPoint(const gp_Pnt& Point,
 523   const Standard_Real t,
 524   const Handle(Adaptor3d_HCurve)& C,
 525   const Handle(Adaptor3d_HSurface)& S,
 526   const Standard_Real TolU,
 527   const Standard_Real TolV,
 528   Standard_Real& U,
 529   Standard_Real& V)
 530 {
 531
 532   ProjLib_PrjResolve aPrjPS(C->Curve(), S->Surface(), 1);
 533   Standard_Real ParU,ParV;
 534   Extrema_ExtPS aExtPS;
 535   aExtPS.Initialize(S->Surface(), S->FirstUParameter(),
 536     S->LastUParameter(), S->FirstVParameter(),
 537     S->LastVParameter(), TolU, TolV);
 538
 539   aExtPS.Perform(Point);
 540   Standard_Integer argmin = 0;
 541   if (aExtPS.IsDone() && aExtPS.NbExt())
 542   {
 543     Standard_Integer i, Nend;
 544     // Search for the nearest solution which is also a normal projection
 545     Nend = aExtPS.NbExt();
 546     for(i = 1; i <= Nend; i++)
 547     {
 548       Extrema_POnSurf POnS = aExtPS.Point(i);
 549       POnS.Parameter(ParU, ParV);
 550       aPrjPS.Perform(t, ParU, ParV, gp_Pnt2d(TolU, TolV),
 551         gp_Pnt2d(S->FirstUParameter(), S->FirstVParameter()),
 552         gp_Pnt2d(S->LastUParameter(), S->LastVParameter()),
 553         FuncTol, Standard_True);
 554       if(aPrjPS.IsDone() )
 555         if  (argmin == 0 || aExtPS.SquareDistance(i) < aExtPS.SquareDistance(argmin)) argmin = i;
 556     }
 557   }
 558   if( argmin == 0 ) return Standard_False;
 559   else
 560   {
 561     Extrema_POnSurf POnS = aExtPS.Point(argmin);
 562     POnS.Parameter(U, V);
 563     return Standard_True;
 564   }
 565 }
 566
 567 //=======================================================================
 568 //function : ProjLib_CompProjectedCurve
 569 //purpose  :
 570 //=======================================================================
 571
 572 ProjLib_CompProjectedCurve::ProjLib_CompProjectedCurve()
 573 : myNbCurves(0),
 574   myTolU    (0.0),
 575   myTolV    (0.0),
 576   myMaxDist (0.0)
 577 {
 578 }
 579
 580 //=======================================================================
 581 //function : ProjLib_CompProjectedCurve
 582 //purpose  :
 583 //=======================================================================
 584
 585 ProjLib_CompProjectedCurve::ProjLib_CompProjectedCurve
 586                            (const Handle(Adaptor3d_HSurface)& theSurface,
 587                             const Handle(Adaptor3d_HCurve)&   theCurve,
 588                             const Standard_Real               theTolU,
 589                             const Standard_Real               theTolV)
 590 : mySurface (theSurface),
 591   myCurve   (theCurve),
 592   myNbCurves(0),
 593   mySequence(new ProjLib_HSequenceOfHSequenceOfPnt()),
 594   myTolU    (theTolU),
 595   myTolV    (theTolV),
 596   myMaxDist (-1.0)
 597 {
 598   Init();
 599 }
 600
 601 //=======================================================================
 602 //function : ProjLib_CompProjectedCurve
 603 //purpose  :
 604 //=======================================================================
 605
 606 ProjLib_CompProjectedCurve::ProjLib_CompProjectedCurve
 607                            (const Handle(Adaptor3d_HSurface)& theSurface,
 608                             const Handle(Adaptor3d_HCurve)&   theCurve,
 609                             const Standard_Real               theTolU,
 610                             const Standard_Real               theTolV,
 611                             const Standard_Real               theMaxDist)
 612 : mySurface (theSurface),
 613   myCurve   (theCurve),
 614   myNbCurves(0),
 615   mySequence(new ProjLib_HSequenceOfHSequenceOfPnt()),
 616   myTolU    (theTolU),
 617   myTolV    (theTolV),
 618   myMaxDist (theMaxDist)
 619 {
 620   Init();
 621 }
 622
 623 //=======================================================================
 624 //function : Init
 625 //purpose  :
 626 //=======================================================================
 627
 628 void ProjLib_CompProjectedCurve::Init()
 629 {
 630   myTabInt.Nullify();
 631   NCollection_Vector<Standard_Real> aSplits;
 632   aSplits.Clear();
 633
 634   Standard_Real Tol;// Tolerance for ExactBound
 635   Standard_Integer i, Nend = 0, aSplitIdx = 0;
 636   Standard_Boolean FromLastU = Standard_False,
 637                    isSplitsComputed = Standard_False;
 638
 639   const Standard_Real aTol3D = Precision::Confusion();
 640   Extrema_ExtCS CExt(myCurve->Curve(), mySurface->Surface(), aTol3D, aTol3D);
 641   if (CExt.IsDone() && CExt.NbExt())
 642   {
 643     // Search for the minimum solution.
 644     // Avoid usage of extrema result that can be wrong for extrusion.
 645     if(myMaxDist > 0 &&
 646
 647        mySurface->GetType() != GeomAbs_SurfaceOfExtrusion)
 648     {
 649       Standard_Real min_val2;
 650       min_val2 = CExt.SquareDistance(1);
 651
 652       Nend = CExt.NbExt();
 653       for(i = 2; i <= Nend; i++)
 654       {
 655         if (CExt.SquareDistance(i) < min_val2)
 656           min_val2 = CExt.SquareDistance(i);
 657       }
 658       if (min_val2 > myMaxDist * myMaxDist)
 659         return; // No near solution -> exit.
 660     }
 661   }
 662
 663   Standard_Real FirstU, LastU, Step, SearchStep, WalkStep, t;
 664
 665   FirstU = myCurve->FirstParameter();
 666   LastU  = myCurve->LastParameter();
 667   const Standard_Real GlobalMinStep = 1.e-4;
 668   //<GlobalMinStep> is sufficiently small to provide solving from initial point
 669   //and, on the other hand, it is sufficiently large to avoid too close solutions.
 670   const Standard_Real MinStep = 0.01*(LastU - FirstU),
 671     MaxStep = 0.1*(LastU - FirstU);
 672   SearchStep = 10*MinStep;
 673   Step = SearchStep;
 674
 675   gp_Pnt2d aLowBorder(mySurface->FirstUParameter(),mySurface->FirstVParameter());
 676   gp_Pnt2d aUppBorder(mySurface->LastUParameter(), mySurface->LastVParameter());
 677   gp_Pnt2d aTol(myTolU, myTolV);
 678   ProjLib_PrjResolve aPrjPS(myCurve->Curve(), mySurface->Surface(), 1);
 679
 680   t = FirstU;
 681   Standard_Boolean new_part;
 682   Standard_Real prevDeb=0.;
 683   Standard_Boolean SameDeb=Standard_False;
 684
 685
 686   gp_Pnt Triple, prevTriple;
 687
 688   //Basic loop
 689   while(t <= LastU)
 690   {
 691     // Search for the beginning of a new continuous part
 692     // to avoid infinite computation in some difficult cases.
 693     new_part = Standard_False;
 694     if(t > FirstU && Abs(t-prevDeb) <= Precision::PConfusion()) SameDeb=Standard_True;
 695     while(t <= LastU && !new_part && !FromLastU && !SameDeb)
 696     {
 697       prevDeb=t;
 698       if (t == LastU) FromLastU=Standard_True;
 699       Standard_Boolean initpoint=Standard_False;
 700       Standard_Real U = 0., V = 0.;
 701       gp_Pnt CPoint;
 702       Standard_Real ParT,ParU,ParV;
 703
 704       // Search an initial point in the list of Extrema Curve-Surface
 705       if(Nend != 0 && !CExt.IsParallel())
 706       {
 707         for (i=1;i<=Nend;i++)
 708         {
 709           Extrema_POnCurv P1;
 710           Extrema_POnSurf P2;
 711           CExt.Points(i,P1,P2);
 712           ParT=P1.Parameter();
 713           P2.Parameter(ParU, ParV);
 714
 715           aPrjPS.Perform(ParT, ParU, ParV, aTol, aLowBorder, aUppBorder, FuncTol, Standard_True);
 716
 717           if ( aPrjPS.IsDone() && P1.Parameter() > Max(FirstU,t-Step+Precision::PConfusion())
 718             && P1.Parameter() <= t)
 719           {
 720             t=ParT;
 721             U=ParU;
 722             V=ParV;
 723             CPoint=P1.Value();
 724             initpoint = Standard_True;
 725             break;
 726           }
 727         }
 728       }
 729       if (!initpoint)
 730       {
 731         myCurve->D0(t,CPoint);
 732 #ifdef OCCT_DEBUG_CHRONO
 733         InitChron(chr_init_point);
 734 #endif
 735         // PConfusion - use geometric tolerances in extrema / optimization.
 736         initpoint=InitialPoint(CPoint, t,myCurve,mySurface, Precision::PConfusion(), Precision::PConfusion(), U, V);
 737 #ifdef OCCT_DEBUG_CHRONO
 738         ResultChron(chr_init_point,t_init_point);
 739         init_point_count++;
 740 #endif
 741       }
 742       if(initpoint)
 743       {
 744         // When U or V lie on surface joint in some cases we cannot use them
 745         // as initial point for aPrjPS, so we switch them
 746         gp_Vec2d D;
 747
 748         if ((mySurface->IsUPeriodic() &&
 749             Abs(aUppBorder.X() - aLowBorder.X() - mySurface->UPeriod()) < Precision::Confusion()) ||
 750             (mySurface->IsVPeriodic() &&
 751             Abs(aUppBorder.Y() - aLowBorder.Y() - mySurface->VPeriod()) < Precision::Confusion()))
 752         {
 753           if((Abs(U - aLowBorder.X()) < mySurface->UResolution(Precision::PConfusion())) &&
 754             mySurface->IsUPeriodic())
 755           {
 756             d1(t, U, V, D, myCurve, mySurface);
 757             if (D.X() < 0 ) U = aUppBorder.X();
 758           }
 759           else if((Abs(U - aUppBorder.X()) < mySurface->UResolution(Precision::PConfusion())) &&
 760             mySurface->IsUPeriodic())
 761           {
 762             d1(t, U, V, D, myCurve, mySurface);
 763             if (D.X() > 0) U = aLowBorder.X();
 764           }
 765
 766           if((Abs(V - aLowBorder.Y()) < mySurface->VResolution(Precision::PConfusion())) &&
 767             mySurface->IsVPeriodic())
 768           {
 769             d1(t, U, V, D, myCurve, mySurface);
 770             if (D.Y() < 0) V = aUppBorder.Y();
 771           }
 772           else if((Abs(V - aUppBorder.Y()) <= mySurface->VResolution(Precision::PConfusion())) &&
 773             mySurface->IsVPeriodic())
 774           {
 775             d1(t, U, V, D, myCurve, mySurface);
 776             if (D.Y() > 0) V = aLowBorder.Y();
 777           }
 778         }
 779
 780         if (myMaxDist > 0)
 781         {
 782           // Here we are going to stop if the distance between projection and
 783           // corresponding curve point is greater than myMaxDist
 784           gp_Pnt POnS;
 785           Standard_Real d;
 786           mySurface->D0(U, V, POnS);
 787           d = CPoint.Distance(POnS);
 788           if (d > myMaxDist)
 789           {
 790             mySequence->Clear();
 791             myNbCurves = 0;
 792             return;
 793           }
 794         }
 795         Triple = gp_Pnt(t, U, V);
 796         if (t != FirstU)
 797         {
 798           //Search for exact boundary point
 799           Tol = Min(myTolU, myTolV);
 800           gp_Vec2d aD;
 801           d1(Triple.X(), Triple.Y(), Triple.Z(), aD, myCurve, mySurface);
 802           Tol /= Max(Abs(aD.X()), Abs(aD.Y()));
 803
 804           if(!ExactBound(Triple, t - Step, Tol,
 805             myTolU, myTolV, myCurve, mySurface))
 806           {
 807 #ifdef OCCT_DEBUG
 808             cout<<"There is a problem with ExactBound computation"<<endl;
 809 #endif
 810             DichExactBound(Triple, t - Step, Tol, myTolU, myTolV,
 811               myCurve, mySurface);
 812           }
 813         }
 814         new_part = Standard_True;
 815       }
 816       else
 817       {
 818         if(t == LastU) break;
 819         t += Step;
 820         if(t>LastU)
 821         {
 822           Step =Step+LastU-t;
 823           t=LastU;
 824         }
 825       }
 826     }
 827     if (!new_part) break;
 828
 829     //We have found a new continuous part
 830     Handle(TColgp_HSequenceOfPnt) hSeq = new TColgp_HSequenceOfPnt();
 831     mySequence->Append(hSeq);
 832     myNbCurves++;
 833     mySequence->Value(myNbCurves)->Append(Triple);
 834     prevTriple = Triple;
 835
 836     if (Triple.X() == LastU) break;//return;
 837
 838     //Computation of WalkStep
 839     gp_Vec D1, D2;
 840     Standard_Real MagnD1, MagnD2;
 841     d2CurvOnSurf(Triple.X(), Triple.Y(), Triple.Z(), D1, D2, myCurve, mySurface);
 842     MagnD1 = D1.Magnitude();
 843     MagnD2 = D2.Magnitude();
 844     if(MagnD2 < Precision::Confusion()) WalkStep = MaxStep;
 845     else WalkStep = Min(MaxStep, Max(MinStep, 0.1*MagnD1/MagnD2));
 846
 847     Step = WalkStep;
 848
 849     t = Triple.X() + Step;
 850     if (t > LastU) t = LastU;
 851     Standard_Real prevStep = Step;
 852     Standard_Real U0, V0;
 853
 854     //Here we are trying to prolong continuous part
 855     while (t <= LastU && new_part)
 856     {
 857
 858       U0 = Triple.Y() + (Step / prevStep) * (Triple.Y() - prevTriple.Y());
 859       V0 = Triple.Z() + (Step / prevStep) * (Triple.Z() - prevTriple.Z());
 860       // adjust U0 to be in [mySurface->FirstUParameter(),mySurface->LastUParameter()]
 861       U0 = Min(Max(U0, aLowBorder.X()), aUppBorder.X());
 862       // adjust V0 to be in [mySurface->FirstVParameter(),mySurface->LastVParameter()]
 863       V0 = Min(Max(V0, aLowBorder.Y()), aUppBorder.Y());
 864
 865
 866       aPrjPS.Perform(t, U0, V0, aTol,
 867                      aLowBorder, aUppBorder, FuncTol, Standard_True);
 868       if(!aPrjPS.IsDone())
 869       {
 870         if (Step <= GlobalMinStep)
 871         {
 872           //Search for exact boundary point
 873           Tol = Min(myTolU, myTolV);
 874           gp_Vec2d D;
 875           d1(Triple.X(), Triple.Y(), Triple.Z(), D, myCurve, mySurface);
 876           Tol /= Max(Abs(D.X()), Abs(D.Y()));
 877
 878           if(!ExactBound(Triple, t, Tol, myTolU, myTolV,
 879             myCurve, mySurface))
 880           {
 881 #ifdef OCCT_DEBUG
 882             cout<<"There is a problem with ExactBound computation"<<endl;
 883 #endif
 884             DichExactBound(Triple, t, Tol, myTolU, myTolV,
 885               myCurve, mySurface);
 886           }
 887
 888           if((Triple.X() - mySequence->Value(myNbCurves)->Value(mySequence->Value(myNbCurves)->Length()).X()) > 1.e-10)
 889             mySequence->Value(myNbCurves)->Append(Triple);
 890           if((LastU - Triple.X()) < Tol) {t = LastU + 1; break;}//return;
 891
 892           Step = SearchStep;
 893           t = Triple.X() + Step;
 894           if (t > (LastU-MinStep/2) )
 895           {
 896             Step =Step+LastU-t;
 897             t = LastU;
 898           }
 899           new_part = Standard_False;
 900         }
 901         else
 902         {
 903           // decrease step
 904           Standard_Real SaveStep = Step;
 905           Step /= 2.;
 906           t = Triple .X() + Step;
 907           if (t > (LastU-MinStep/4) )
 908           {
 909             Step =Step+LastU-t;
 910             if (Abs(Step - SaveStep) <= Precision::PConfusion())
 911               Step = GlobalMinStep; //to avoid looping
 912             t = LastU;
 913           }
 914         }
 915       }
 916       // Go further
 917       else
 918       {
 919         prevTriple = Triple;
 920         prevStep = Step;
 921         Triple = gp_Pnt(t, aPrjPS.Solution().X(), aPrjPS.Solution().Y());
 922
 923         // Check for possible local traps.
 924         UpdateTripleByTrapCriteria(Triple);
 925
 926         // Protection from case when the whole curve lies on a seam.
 927         if (!isSplitsComputed)
 928         {
 929           Standard_Boolean isUPossible = Standard_False;
 930           if (mySurface->IsUPeriodic() &&
 931              (Abs(Triple.Y() - mySurface->FirstUParameter() ) > Precision::PConfusion() &&
 932               Abs(Triple.Y() - mySurface->LastUParameter()  ) > Precision::PConfusion()))
 933           {
 934             isUPossible = Standard_True;
 935           }
 936
 937           Standard_Boolean isVPossible = Standard_False;
 938           if (mySurface->IsVPeriodic() &&
 939              (Abs(Triple.Z() - mySurface->FirstVParameter() ) > Precision::PConfusion() &&
 940               Abs(Triple.Z() - mySurface->LastVParameter()  ) > Precision::PConfusion()))
 941           {
 942             isVPossible = Standard_True;
 943           }
 944
 945           if (isUPossible || isVPossible)
 946           {
 947             // When point is good conditioned.
 948             BuildCurveSplits(myCurve, mySurface, myTolU, myTolV, aSplits);
 949             isSplitsComputed = Standard_True;
 950           }
 951         }
 952
 953         if((Triple.X() - mySequence->Value(myNbCurves)->Value(mySequence->Value(myNbCurves)->Length()).X()) > 1.e-10)
 954           mySequence->Value(myNbCurves)->Append(Triple);
 955         if (t == LastU) {t = LastU + 1; break;}//return;
 956         //Computation of WalkStep
 957         d2CurvOnSurf(Triple.X(), Triple.Y(), Triple.Z(), D1, D2, myCurve, mySurface);
 958         MagnD1 = D1.Magnitude();
 959         MagnD2 = D2.Magnitude();
 960         if(MagnD2 < Precision::Confusion() ) WalkStep = MaxStep;
 961         else WalkStep = Min(MaxStep, Max(MinStep, 0.1*MagnD1/MagnD2));
 962
 963         Step = WalkStep;
 964         t += Step;
 965         if (t > (LastU-MinStep/2))
 966         {
 967           Step = Step + LastU - t;
 968           t = LastU;
 969         }
 970
 971         // We assume at least one point of cache inside of a split.
 972         const Standard_Integer aSize = aSplits.Size();
 973         for(Standard_Integer anIdx = aSplitIdx; anIdx < aSize; ++anIdx)
 974         {
 975           const Standard_Real aParam = aSplits(anIdx);
 976           if (Abs(aParam - Triple.X() ) < Precision::PConfusion())
 977           {
 978             // The current point is equal to a split point.
 979             new_part = Standard_False;
 980
 981             // Move split index to avoid check of the whole list.
 982             ++aSplitIdx;
 983             break;
 984           }
 985           else if (aParam < t + Precision::PConfusion() )
 986           {
 987             // The next point crosses the split point.
 988             t = aParam;
 989             Step = t - prevTriple.X();
 990           }
 991         } // for(Standard_Integer anIdx = aSplitIdx; anIdx < aSize; ++anIdx)
 992       }
 993     }
 994   }
 995
 996   // Sequence post-proceeding.
 997   Standard_Integer j;
 998
 999   // 1. Removing poor parts
1000   Standard_Integer NbPart=myNbCurves;
1001   Standard_Integer ipart=1;
1002   for(i = 1; i <= NbPart; i++) {
1003     //    Standard_Integer NbPoints = mySequence->Value(i)->Length();
1004     if(mySequence->Value(ipart)->Length() < 2) {
1005       mySequence->Remove(ipart);
1006       myNbCurves--;
1007     }
1008     else ipart++;
1009   }
1010
1011   if(myNbCurves == 0) return;
1012
1013   // 2. Removing common parts of bounds
1014   for(i = 1; i < myNbCurves; i++)
1015   {
1016     if(mySequence->Value(i)->Value(mySequence->Value(i)->Length()).X() >=
1017       mySequence->Value(i+1)->Value(1).X())
1018     {
1019       mySequence->ChangeValue(i+1)->ChangeValue(1).SetX(mySequence->Value(i)->Value(mySequence->Value(i)->Length()).X() + 1.e-12);
1020     }
1021   }
1022
1023   // 3. Computation of the maximum distance from each part of curve to surface
1024
1025   myMaxDistance = new TColStd_HArray1OfReal(1, myNbCurves);
1026   myMaxDistance->Init(0);
1027   for(i = 1; i <= myNbCurves; i++)
1028   {
1029     for(j = 1; j <= mySequence->Value(i)->Length(); j++)
1030     {
1031       gp_Pnt POnC, POnS, aTriple;
1032       Standard_Real Distance;
1033       aTriple = mySequence->Value(i)->Value(j);
1034       myCurve->D0(aTriple.X(), POnC);
1035       mySurface->D0(aTriple.Y(), aTriple.Z(), POnS);
1036       Distance = POnC.Distance(POnS);
1037       if (myMaxDistance->Value(i) < Distance)
1038       {
1039         myMaxDistance->ChangeValue(i) = Distance;
1040       }
1041     }
1042   }
1043
1044   // 4. Check the projection to be a single point
1045
1046   gp_Pnt2d Pmoy, Pcurr, P;
1047   Standard_Real AveU, AveV;
1048   mySnglPnts = new TColStd_HArray1OfBoolean(1, myNbCurves);
1049   mySnglPnts->Init (Standard_True);
1050
1051   for(i = 1; i <= myNbCurves; i++)
1052   {
1053     //compute an average U and V
1054
1055     for(j = 1, AveU = 0., AveV = 0.; j <= mySequence->Value(i)->Length(); j++)
1056     {
1057       AveU += mySequence->Value(i)->Value(j).Y();
1058       AveV += mySequence->Value(i)->Value(j).Z();
1059     }
1060     AveU /= mySequence->Value(i)->Length();
1061     AveV /= mySequence->Value(i)->Length();
1062
1063     Pmoy.SetCoord(AveU,AveV);
1064     for(j = 1; j <= mySequence->Value(i)->Length(); j++)
1065     {
1066       Pcurr =
1067         gp_Pnt2d(mySequence->Value(i)->Value(j).Y(), mySequence->Value(i)->Value(j).Z());
1068       if (Pcurr.Distance(Pmoy) > ((myTolU < myTolV) ? myTolV : myTolU))
1069       {
1070         mySnglPnts->SetValue(i, Standard_False);
1071         break;
1072       }
1073     }
1074   }
1075
1076   // 5. Check the projection to be an isoparametric curve of the surface
1077
1078   myUIso = new TColStd_HArray1OfBoolean(1, myNbCurves);
1079   myUIso->Init (Standard_True);
1080
1081   myVIso = new TColStd_HArray1OfBoolean(1, myNbCurves);
1082   myVIso->Init (Standard_True);
1083
1084   for(i = 1; i <= myNbCurves; i++) {
1085     if (IsSinglePnt(i, P)|| mySequence->Value(i)->Length() <=2) {
1086       myUIso->SetValue(i, Standard_False);
1087       myVIso->SetValue(i, Standard_False);
1088       continue;
1089     }
1090
1091     // new test for isoparametrics
1092
1093     if ( mySequence->Value(i)->Length() > 2) {
1094       //compute an average U and V
1095
1096       for(j = 1, AveU = 0., AveV = 0.; j <= mySequence->Value(i)->Length(); j++) {
1097         AveU += mySequence->Value(i)->Value(j).Y();
1098         AveV += mySequence->Value(i)->Value(j).Z();
1099       }
1100       AveU /= mySequence->Value(i)->Length();
1101       AveV /= mySequence->Value(i)->Length();
1102
1103       // is i-part U-isoparametric ?
1104       for(j = 1; j <= mySequence->Value(i)->Length(); j++)
1105       {
1106         if(Abs(mySequence->Value(i)->Value(j).Y() - AveU) > myTolU)
1107         {
1108           myUIso->SetValue(i, Standard_False);
1109           break;
1110         }
1111       }
1112
1113       // is i-part V-isoparametric ?
1114       for(j = 1; j <= mySequence->Value(i)->Length(); j++)
1115       {
1116         if(Abs(mySequence->Value(i)->Value(j).Z() - AveV) > myTolV)
1117         {
1118           myVIso->SetValue(i, Standard_False);
1119           break;
1120         }
1121       }
1122       //
1123     }
1124   }
1125 }
1126 //=======================================================================
1127 //function : Load
1128 //purpose  :
1129 //=======================================================================
1130
1131 void ProjLib_CompProjectedCurve::Load(const Handle(Adaptor3d_HSurface)& S)
1132 {
1133   mySurface = S;
1134 }
1135
1136 //=======================================================================
1137 //function : Load
1138 //purpose  :
1139 //=======================================================================
1140
1141 void ProjLib_CompProjectedCurve::Load(const Handle(Adaptor3d_HCurve)& C)
1142 {
1143   myCurve = C;
1144 }
1145
1146 //=======================================================================
1147 //function : GetSurface
1148 //purpose  :
1149 //=======================================================================
1150
1151 const Handle(Adaptor3d_HSurface)& ProjLib_CompProjectedCurve::GetSurface() const
1152 {
1153   return mySurface;
1154 }
1155
1156
1157 //=======================================================================
1158 //function : GetCurve
1159 //purpose  :
1160 //=======================================================================
1161
1162 const Handle(Adaptor3d_HCurve)& ProjLib_CompProjectedCurve::GetCurve() const
1163 {
1164   return myCurve;
1165 }
1166
1167 //=======================================================================
1168 //function : GetTolerance
1169 //purpose  :
1170 //=======================================================================
1171
1172 void ProjLib_CompProjectedCurve::GetTolerance(Standard_Real& TolU,
1173   Standard_Real& TolV) const
1174 {
1175   TolU = myTolU;
1176   TolV = myTolV;
1177 }
1178
1179 //=======================================================================
1180 //function : NbCurves
1181 //purpose  :
1182 //=======================================================================
1183
1184 Standard_Integer ProjLib_CompProjectedCurve::NbCurves() const
1185 {
1186   return myNbCurves;
1187 }
1188 //=======================================================================
1189 //function : Bounds
1190 //purpose  :
1191 //=======================================================================
1192
1193 void ProjLib_CompProjectedCurve::Bounds(const Standard_Integer Index,
1194   Standard_Real& Udeb,
1195   Standard_Real& Ufin) const
1196 {
1197   if(Index < 1 || Index > myNbCurves) throw Standard_NoSuchObject();
1198   Udeb = mySequence->Value(Index)->Value(1).X();
1199   Ufin = mySequence->Value(Index)->Value(mySequence->Value(Index)->Length()).X();
1200 }
1201 //=======================================================================
1202 //function : IsSinglePnt
1203 //purpose  :
1204 //=======================================================================
1205
1206 Standard_Boolean ProjLib_CompProjectedCurve::IsSinglePnt(const Standard_Integer Index, gp_Pnt2d& P) const
1207 {
1208   if(Index < 1 || Index > myNbCurves) throw Standard_NoSuchObject();
1209   P = gp_Pnt2d(mySequence->Value(Index)->Value(1).Y(), mySequence->Value(Index)->Value(1).Z());
1210   return mySnglPnts->Value(Index);
1211 }
1212
1213 //=======================================================================
1214 //function : IsUIso
1215 //purpose  :
1216 //=======================================================================
1217
1218 Standard_Boolean ProjLib_CompProjectedCurve::IsUIso(const Standard_Integer Index, Standard_Real& U) const
1219 {
1220   if(Index < 1 || Index > myNbCurves) throw Standard_NoSuchObject();
1221   U = mySequence->Value(Index)->Value(1).Y();
1222   return myUIso->Value(Index);
1223 }
1224 //=======================================================================
1225 //function : IsVIso
1226 //purpose  :
1227 //=======================================================================
1228
1229 Standard_Boolean ProjLib_CompProjectedCurve::IsVIso(const Standard_Integer Index, Standard_Real& V) const
1230 {
1231   if(Index < 1 || Index > myNbCurves) throw Standard_NoSuchObject();
1232   V = mySequence->Value(Index)->Value(1).Z();
1233   return myVIso->Value(Index);
1234 }
1235 //=======================================================================
1236 //function : Value
1237 //purpose  :
1238 //=======================================================================
1239
1240 gp_Pnt2d ProjLib_CompProjectedCurve::Value(const Standard_Real t) const
1241 {
1242   gp_Pnt2d P;
1243   D0(t, P);
1244   return P;
1245 }
1246 //=======================================================================
1247 //function : D0
1248 //purpose  :
1249 //=======================================================================
1250
1251 void ProjLib_CompProjectedCurve::D0(const Standard_Real U,gp_Pnt2d& P) const
1252 {
1253   Standard_Integer i, j;
1254   Standard_Real Udeb, Ufin;
1255   Standard_Boolean found = Standard_False;
1256
1257   for(i = 1; i <= myNbCurves; i++)
1258   {
1259     Bounds(i, Udeb, Ufin);
1260     if (U >= Udeb && U <= Ufin)
1261     {
1262       found = Standard_True;
1263       break;
1264     }
1265   }
1266   if (!found) throw Standard_DomainError("ProjLib_CompProjectedCurve::D0");
1267
1268   Standard_Real U0, V0;
1269
1270   Standard_Integer End = mySequence->Value(i)->Length();
1271   for(j = 1; j < End; j++)
1272     if ((U >= mySequence->Value(i)->Value(j).X()) && (U <= mySequence->Value(i)->Value(j + 1).X())) break;
1273
1274   //  U0 = mySequence->Value(i)->Value(j).Y();
1275   //  V0 = mySequence->Value(i)->Value(j).Z();
1276
1277   //  Cubic Interpolation
1278   if(mySequence->Value(i)->Length() < 4 ||
1279     (Abs(U-mySequence->Value(i)->Value(j).X()) <= Precision::PConfusion()) )
1280   {
1281     U0 = mySequence->Value(i)->Value(j).Y();
1282     V0 = mySequence->Value(i)->Value(j).Z();
1283   }
1284   else if (Abs(U-mySequence->Value(i)->Value(j+1).X())
1285     <= Precision::PConfusion())
1286   {
1287     U0 = mySequence->Value(i)->Value(j+1).Y();
1288     V0 = mySequence->Value(i)->Value(j+1).Z();
1289   }
1290   else
1291   {
1292     if (j == 1) j = 2;
1293     if (j > mySequence->Value(i)->Length() - 2)
1294       j = mySequence->Value(i)->Length() - 2;
1295
1296     gp_Vec2d I1, I2, I3, I21, I22, I31, Y1, Y2, Y3, Y4, Res;
1297     Standard_Real X1, X2, X3, X4;
1298
1299     X1 = mySequence->Value(i)->Value(j - 1).X();
1300     X2 = mySequence->Value(i)->Value(j).X();
1301     X3 = mySequence->Value(i)->Value(j + 1).X();
1302     X4 = mySequence->Value(i)->Value(j + 2).X();
1303
1304     Y1 = gp_Vec2d(mySequence->Value(i)->Value(j - 1).Y(),
1305       mySequence->Value(i)->Value(j - 1).Z());
1306     Y2 = gp_Vec2d(mySequence->Value(i)->Value(j).Y(),
1307       mySequence->Value(i)->Value(j).Z());
1308     Y3 = gp_Vec2d(mySequence->Value(i)->Value(j + 1).Y(),
1309       mySequence->Value(i)->Value(j + 1).Z());
1310     Y4 = gp_Vec2d(mySequence->Value(i)->Value(j + 2).Y(),
1311       mySequence->Value(i)->Value(j + 2).Z());
1312
1313     I1 = (Y1 - Y2)/(X1 - X2);
1314     I2 = (Y2 - Y3)/(X2 - X3);
1315     I3 = (Y3 - Y4)/(X3 - X4);
1316
1317     I21 = (I1 - I2)/(X1 - X3);
1318     I22 = (I2 - I3)/(X2 - X4);
1319
1320     I31 = (I21 - I22)/(X1 - X4);
1321
1322     Res = Y1 + (U - X1)*(I1 + (U - X2)*(I21 + (U - X3)*I31));
1323
1324     U0 = Res.X();
1325     V0 = Res.Y();
1326
1327     if(U0 < mySurface->FirstUParameter()) U0 = mySurface->FirstUParameter();
1328     else if(U0 > mySurface->LastUParameter()) U0 = mySurface->LastUParameter();
1329
1330     if(V0 < mySurface->FirstVParameter()) V0 = mySurface->FirstVParameter();
1331     else if(V0 > mySurface->LastVParameter()) V0 = mySurface->LastVParameter();
1332   }
1333   //End of cubic interpolation
1334
1335   ProjLib_PrjResolve aPrjPS(myCurve->Curve(), mySurface->Surface(), 1);
1336   aPrjPS.Perform(U, U0, V0, gp_Pnt2d(myTolU, myTolV),
1337     gp_Pnt2d(mySurface->FirstUParameter(), mySurface->FirstVParameter()),
1338     gp_Pnt2d(mySurface->LastUParameter(), mySurface->LastVParameter()));
1339   if (aPrjPS.IsDone())
1340     P = aPrjPS.Solution();
1341   else
1342   {
1343     gp_Pnt thePoint = myCurve->Value(U);
1344     Extrema_ExtPS aExtPS(thePoint, mySurface->Surface(), myTolU, myTolV);
1345     if (aExtPS.IsDone() && aExtPS.NbExt())
1346     {
1347       Standard_Integer k, Nend, imin = 1;
1348       // Search for the nearest solution which is also a normal projection
1349       Nend = aExtPS.NbExt();
1350       for(k = 2; k <= Nend; k++)
1351         if (aExtPS.SquareDistance(k) < aExtPS.SquareDistance(imin))
1352           imin = k;
1353       const Extrema_POnSurf& POnS = aExtPS.Point(imin);
1354       Standard_Real ParU,ParV;
1355       POnS.Parameter(ParU, ParV);
1356       P.SetCoord(ParU, ParV);
1357     }
1358     else
1359       P.SetCoord(U0,V0);
1360   }
1361 }
1362 //=======================================================================
1363 //function : D1
1364 //purpose  :
1365 //=======================================================================
1366
1367 void ProjLib_CompProjectedCurve::D1(const Standard_Real t,
1368   gp_Pnt2d& P,
1369   gp_Vec2d& V) const
1370 {
1371   Standard_Real u, v;
1372   D0(t, P);
1373   u = P.X();
1374   v = P.Y();
1375   d1(t, u, v, V, myCurve, mySurface);
1376 }
1377 //=======================================================================
1378 //function : D2
1379 //purpose  :
1380 //=======================================================================
1381
1382 void ProjLib_CompProjectedCurve::D2(const Standard_Real t,
1383   gp_Pnt2d& P,
1384   gp_Vec2d& V1,
1385   gp_Vec2d& V2) const
1386 {
1387   Standard_Real u, v;
1388   D0(t, P);
1389   u = P.X();
1390   v = P.Y();
1391   d2(t, u, v, V1, V2, myCurve, mySurface);
1392 }
1393 //=======================================================================
1394 //function : DN
1395 //purpose  :
1396 //=======================================================================
1397
1398 gp_Vec2d ProjLib_CompProjectedCurve::DN(const Standard_Real t,
1399   const Standard_Integer N) const
1400 {
1401   if (N < 1 ) throw Standard_OutOfRange("ProjLib_CompProjectedCurve : N must be greater than 0");
1402   else if (N ==1)
1403   {
1404     gp_Pnt2d P;
1405     gp_Vec2d V;
1406     D1(t,P,V);
1407     return V;
1408   }
1409   else if ( N==2)
1410   {
1411     gp_Pnt2d P;
1412     gp_Vec2d V1,V2;
1413     D2(t,P,V1,V2);
1414     return V2;
1415   }
1416   else if (N > 2 )
1417     throw Standard_NotImplemented("ProjLib_CompProjectedCurve::DN");
1418   return gp_Vec2d();
1419 }
1420
1421 //=======================================================================
1422 //function : GetSequence
1423 //purpose  :
1424 //=======================================================================
1425
1426 const Handle(ProjLib_HSequenceOfHSequenceOfPnt)& ProjLib_CompProjectedCurve::GetSequence() const
1427 {
1428   return mySequence;
1429 }
1430 //=======================================================================
1431 //function : FirstParameter
1432 //purpose  :
1433 //=======================================================================
1434
1435 Standard_Real ProjLib_CompProjectedCurve::FirstParameter() const
1436 {
1437   return myCurve->FirstParameter();
1438 }
1439
1440 //=======================================================================
1441 //function : LastParameter
1442 //purpose  :
1443 //=======================================================================
1444
1445 Standard_Real ProjLib_CompProjectedCurve::LastParameter() const
1446 {
1447   return myCurve->LastParameter();
1448 }
1449
1450 //=======================================================================
1451 //function : MaxDistance
1452 //purpose  :
1453 //=======================================================================
1454
1455 Standard_Real ProjLib_CompProjectedCurve::MaxDistance(const Standard_Integer Index) const
1456 {
1457   if(Index < 1 || Index > myNbCurves) throw Standard_NoSuchObject();
1458   return myMaxDistance->Value(Index);
1459 }
1460
1461 //=======================================================================
1462 //function : NbIntervals
1463 //purpose  :
1464 //=======================================================================
1465
1466 Standard_Integer ProjLib_CompProjectedCurve::NbIntervals(const GeomAbs_Shape S) const
1467 {
1468   const_cast<ProjLib_CompProjectedCurve*>(this)->myTabInt.Nullify();
1469   BuildIntervals(S);
1470   return myTabInt->Length() - 1;
1471 }
1472
1473 //=======================================================================
1474 //function : Intervals
1475 //purpose  :
1476 //=======================================================================
1477
1478 void ProjLib_CompProjectedCurve::Intervals(TColStd_Array1OfReal& T,const GeomAbs_Shape S) const
1479 {
1480   if (myTabInt.IsNull()) BuildIntervals (S);
1481   T = myTabInt->Array1();
1482 }
1483
1484 //=======================================================================
1485 //function : BuildIntervals
1486 //purpose  :
1487 //=======================================================================
1488
1489 void ProjLib_CompProjectedCurve::BuildIntervals(const GeomAbs_Shape S) const
1490 {
1491   GeomAbs_Shape SforS = GeomAbs_CN;
1492   switch(S) {
1493   case GeomAbs_C0:
1494     SforS = GeomAbs_C1;
1495     break;
1496   case GeomAbs_C1:
1497     SforS = GeomAbs_C2;
1498     break;
1499   case GeomAbs_C2:
1500     SforS = GeomAbs_C3;
1501     break;
1502   case GeomAbs_C3:
1503     SforS = GeomAbs_CN;
1504     break;
1505   case GeomAbs_CN:
1506     SforS = GeomAbs_CN;
1507     break;
1508   default:
1509     throw Standard_OutOfRange();
1510   }
1511   Standard_Integer i, j, k;
1512   Standard_Integer NbIntCur = myCurve->NbIntervals(S);
1513   Standard_Integer NbIntSurU = mySurface->NbUIntervals(SforS);
1514   Standard_Integer NbIntSurV = mySurface->NbVIntervals(SforS);
1515
1516   TColStd_Array1OfReal CutPntsT(1, NbIntCur+1);
1517   TColStd_Array1OfReal CutPntsU(1, NbIntSurU+1);
1518   TColStd_Array1OfReal CutPntsV(1, NbIntSurV+1);
1519
1520   myCurve->Intervals(CutPntsT, S);
1521   mySurface->UIntervals(CutPntsU, SforS);
1522   mySurface->VIntervals(CutPntsV, SforS);
1523
1524   Standard_Real Tl, Tr, Ul, Ur, Vl, Vr, Tol;
1525
1526   Handle(TColStd_HArray1OfReal) BArr = NULL,
1527     CArr = NULL,
1528     UArr = NULL,
1529     VArr = NULL;
1530
1531   // proccessing projection bounds
1532   BArr = new TColStd_HArray1OfReal(1, 2*myNbCurves);
1533   for(i = 1; i <= myNbCurves; i++)
1534   {
1535     Bounds(i, BArr->ChangeValue(2*i - 1), BArr->ChangeValue(2*i));
1536   }
1537
1538   // proccessing curve discontinuities
1539   if(NbIntCur > 1) {
1540     CArr = new TColStd_HArray1OfReal(1, NbIntCur - 1);
1541     for(i = 1; i <= CArr->Length(); i++)
1542     {
1543       CArr->ChangeValue(i) = CutPntsT(i + 1);
1544     }
1545   }
1546
1547   // proccessing U-surface discontinuities
1548   TColStd_SequenceOfReal TUdisc;
1549
1550   for(k = 2; k <= NbIntSurU; k++) {
1551     //    cout<<"CutPntsU("<<k<<") = "<<CutPntsU(k)<<endl;
1552     for(i = 1; i <= myNbCurves; i++)
1553     {
1554       for(j = 1; j < mySequence->Value(i)->Length(); j++)
1555       {
1556         Ul = mySequence->Value(i)->Value(j).Y();
1557         Ur = mySequence->Value(i)->Value(j + 1).Y();
1558
1559         if(Abs(Ul - CutPntsU(k)) <= myTolU)
1560           TUdisc.Append(mySequence->Value(i)->Value(j).X());
1561         else if(Abs(Ur - CutPntsU(k)) <= myTolU)
1562           TUdisc.Append(mySequence->Value(i)->Value(j + 1).X());
1563         else if((Ul < CutPntsU(k) && CutPntsU(k) < Ur) ||
1564           (Ur < CutPntsU(k) && CutPntsU(k) < Ul))
1565         {
1566           Standard_Real V;
1567           V = (mySequence->Value(i)->Value(j).Z()
1568             + mySequence->Value(i)->Value(j +1).Z())/2;
1569           ProjLib_PrjResolve Solver(myCurve->Curve(), mySurface->Surface(), 2);
1570
1571           gp_Vec2d D;
1572           gp_Pnt Triple;
1573           Triple = mySequence->Value(i)->Value(j);
1574           d1(Triple.X(), Triple.Y(), Triple.Z(), D, myCurve, mySurface);
1575           if (Abs(D.X()) < Precision::Confusion())
1576             Tol = myTolU;
1577           else
1578             Tol = Min(myTolU, myTolU / Abs(D.X()));
1579
1580           Tl = mySequence->Value(i)->Value(j).X();
1581           Tr = mySequence->Value(i)->Value(j + 1).X();
1582
1583           Solver.Perform((Tl + Tr)/2, CutPntsU(k), V,
1584             gp_Pnt2d(Tol, myTolV),
1585             gp_Pnt2d(Tl, mySurface->FirstVParameter()),
1586             gp_Pnt2d(Tr, mySurface->LastVParameter()));
1587           //
1588           if(Solver.IsDone())
1589           {
1590             TUdisc.Append(Solver.Solution().X());
1591           }
1592         }
1593       }
1594     }
1595   }
1596   for(i = 2; i <= TUdisc.Length(); i++)
1597   {
1598     if(TUdisc(i) - TUdisc(i-1) < Precision::PConfusion())
1599     {
1600       TUdisc.Remove(i--);
1601     }
1602   }
1603
1604   if(TUdisc.Length())
1605   {
1606     UArr = new TColStd_HArray1OfReal(1, TUdisc.Length());
1607     for(i = 1; i <= UArr->Length(); i++)
1608     {
1609       UArr->ChangeValue(i) = TUdisc(i);
1610     }
1611   }
1612   // proccessing V-surface discontinuities
1613   TColStd_SequenceOfReal TVdisc;
1614
1615   for(k = 2; k <= NbIntSurV; k++)
1616   {
1617     for(i = 1; i <= myNbCurves; i++)
1618     {
1619       //      cout<<"CutPntsV("<<k<<") = "<<CutPntsV(k)<<endl;
1620       for(j = 1; j < mySequence->Value(i)->Length(); j++) {
1621
1622         Vl = mySequence->Value(i)->Value(j).Z();
1623         Vr = mySequence->Value(i)->Value(j + 1).Z();
1624
1625         if(Abs(Vl - CutPntsV(k)) <= myTolV)
1626           TVdisc.Append(mySequence->Value(i)->Value(j).X());
1627         else if (Abs(Vr - CutPntsV(k)) <= myTolV)
1628           TVdisc.Append(mySequence->Value(i)->Value(j + 1).X());
1629         else if((Vl < CutPntsV(k) && CutPntsV(k) < Vr) ||
1630           (Vr < CutPntsV(k) && CutPntsV(k) < Vl))
1631         {
1632           Standard_Real U;
1633           U = (mySequence->Value(i)->Value(j).Y()
1634             + mySequence->Value(i)->Value(j +1).Y())/2;
1635           ProjLib_PrjResolve Solver(myCurve->Curve(), mySurface->Surface(), 3);
1636
1637           gp_Vec2d D;
1638           gp_Pnt Triple;
1639           Triple = mySequence->Value(i)->Value(j);
1640           d1(Triple.X(), Triple.Y(), Triple.Z(), D, myCurve, mySurface);
1641           if (Abs(D.Y()) < Precision::Confusion())
1642             Tol = myTolV;
1643           else
1644             Tol = Min(myTolV, myTolV / Abs(D.Y()));
1645
1646           Tl = mySequence->Value(i)->Value(j).X();
1647           Tr = mySequence->Value(i)->Value(j + 1).X();
1648
1649           Solver.Perform((Tl + Tr)/2, U, CutPntsV(k),
1650             gp_Pnt2d(Tol, myTolV),
1651             gp_Pnt2d(Tl, mySurface->FirstUParameter()),
1652             gp_Pnt2d(Tr, mySurface->LastUParameter()));
1653           //
1654           if(Solver.IsDone())
1655           {
1656             TVdisc.Append(Solver.Solution().X());
1657           }
1658         }
1659       }
1660     }
1661   }
1662
1663   for(i = 2; i <= TVdisc.Length(); i++)
1664   {
1665     if(TVdisc(i) - TVdisc(i-1) < Precision::PConfusion())
1666     {
1667       TVdisc.Remove(i--);
1668     }
1669   }
1670
1671   if(TVdisc.Length())
1672   {
1673     VArr = new TColStd_HArray1OfReal(1, TVdisc.Length());
1674     for(i = 1; i <= VArr->Length(); i++)
1675     {
1676       VArr->ChangeValue(i) = TVdisc(i);
1677     }
1678   }
1679
1680   // fusion
1681   TColStd_SequenceOfReal Fusion;
1682   if(!CArr.IsNull())
1683   {
1684     GeomLib::FuseIntervals(BArr->ChangeArray1(),
1685       CArr->ChangeArray1(),
1686       Fusion, Precision::PConfusion());
1687     BArr = new TColStd_HArray1OfReal(1, Fusion.Length());
1688     for(i = 1; i <= BArr->Length(); i++)
1689     {
1690       BArr->ChangeValue(i) = Fusion(i);
1691     }
1692     Fusion.Clear();
1693   }
1694
1695   if(!UArr.IsNull())
1696   {
1697     GeomLib::FuseIntervals(BArr->ChangeArray1(),
1698       UArr->ChangeArray1(),
1699       Fusion, Precision::PConfusion());
1700     BArr = new TColStd_HArray1OfReal(1, Fusion.Length());
1701     for(i = 1; i <= BArr->Length(); i++)
1702     {
1703       BArr->ChangeValue(i) = Fusion(i);
1704     }
1705     Fusion.Clear();
1706   }
1707
1708   if(!VArr.IsNull())
1709   {
1710     GeomLib::FuseIntervals(BArr->ChangeArray1(),
1711       VArr->ChangeArray1(),
1712       Fusion, Precision::PConfusion());
1713     BArr = new TColStd_HArray1OfReal(1, Fusion.Length());
1714     for(i = 1; i <= BArr->Length(); i++)
1715     {
1716       BArr->ChangeValue(i) = Fusion(i);
1717     }
1718   }
1719
1720   const_cast<ProjLib_CompProjectedCurve*>(this)->myTabInt = new TColStd_HArray1OfReal(1, BArr->Length());
1721   for(i = 1; i <= BArr->Length(); i++)
1722   {
1723     myTabInt->ChangeValue(i) = BArr->Value(i);
1724   }
1725 }
1726
1727 //=======================================================================
1728 //function : Trim
1729 //purpose  :
1730 //=======================================================================
1731
1732 Handle(Adaptor2d_HCurve2d) ProjLib_CompProjectedCurve::Trim
1733   (const Standard_Real First,
1734   const Standard_Real Last,
1735   const Standard_Real Tol) const
1736 {
1737   Handle(ProjLib_HCompProjectedCurve) HCS =
1738     new ProjLib_HCompProjectedCurve(*this);
1739   HCS->ChangeCurve2d().Load(mySurface);
1740   HCS->ChangeCurve2d().Load(myCurve->Trim(First,Last,Tol));
1741   return HCS;
1742 }
1743
1744 //=======================================================================
1745 //function : GetType
1746 //purpose  :
1747 //=======================================================================
1748
1749 GeomAbs_CurveType ProjLib_CompProjectedCurve::GetType() const
1750 {
1751   return GeomAbs_OtherCurve;
1752 }
1753
1754 //=======================================================================
1755 //function : UpdateTripleByTrapCriteria
1756 //purpose  :
1757 //=======================================================================
1758 void ProjLib_CompProjectedCurve::UpdateTripleByTrapCriteria(gp_Pnt &thePoint) const
1759 {
1760   Standard_Boolean isProblemsPossible = Standard_False;
1761   // Check possible traps cases:
1762
1763   // 25892 bug.
1764   if (mySurface->GetType() == GeomAbs_SurfaceOfRevolution)
1765   {
1766     // Compute maximal deviation from 3D and choose the biggest one.
1767     Standard_Real aVRes = mySurface->VResolution(Precision::Confusion());
1768     Standard_Real aMaxTol = Max(Precision::PConfusion(), aVRes);
1769
1770     if (Abs (thePoint.Z() - mySurface->FirstVParameter()) < aMaxTol ||
1771         Abs (thePoint.Z() - mySurface->LastVParameter() ) < aMaxTol )
1772     {
1773       isProblemsPossible = Standard_True;
1774     }
1775   }
1776
1777   // 27135 bug. Trap on degenerated edge.
1778   if (mySurface->GetType() == GeomAbs_Sphere &&
1779      (Abs (thePoint.Z() - mySurface->FirstVParameter()) < Precision::PConfusion() ||
1780       Abs (thePoint.Z() - mySurface->LastVParameter() ) < Precision::PConfusion() ||
1781       Abs (thePoint.Y() - mySurface->FirstUParameter()) < Precision::PConfusion() ||
1782       Abs (thePoint.Y() - mySurface->LastUParameter() ) < Precision::PConfusion() ))
1783   {
1784     isProblemsPossible = Standard_True;
1785   }
1786
1787   if (!isProblemsPossible)
1788     return;
1789
1790   Standard_Real U,V;
1791   Standard_Boolean isDone =
1792     InitialPoint(myCurve->Value(thePoint.X()), thePoint.X(), myCurve, mySurface,
1793                  Precision::PConfusion(), Precision::PConfusion(), U, V);
1794
1795   if (!isDone)
1796     return;
1797
1798   // Restore original position in case of period jump.
1799   if (mySurface->IsUPeriodic() &&
1800       Abs (Abs(U - thePoint.Y()) - mySurface->UPeriod()) < Precision::PConfusion())
1801   {
1802     U = thePoint.Y();
1803   }
1804   if (mySurface->IsVPeriodic() &&
1805       Abs (Abs(V - thePoint.Z()) - mySurface->VPeriod()) < Precision::PConfusion())
1806   {
1807     V = thePoint.Z();
1808   }
1809   thePoint.SetY(U);
1810   thePoint.SetZ(V);
1811 }
1812
1813 //=======================================================================
1814 //function : BuildCurveSplits
1815 //purpose  :
1816 //=======================================================================
1817 void BuildCurveSplits(const Handle(Adaptor3d_HCurve)   &theCurve,
1818                       const Handle(Adaptor3d_HSurface) &theSurface,
1819                       const Standard_Real theTolU,
1820                       const Standard_Real theTolV,
1821                       NCollection_Vector<Standard_Real> &theSplits)
1822 {
1823   SplitDS aDS(theCurve, theSurface, theSplits);
1824
1825   Extrema_ExtPS anExtPS;
1826   anExtPS.Initialize(theSurface->Surface(),
1827                      theSurface->FirstUParameter(), theSurface->LastUParameter(),
1828                      theSurface->FirstVParameter(), theSurface->LastVParameter(),
1829                      theTolU, theTolV);
1830   aDS.myExtPS = &anExtPS;
1831
1832   if (theSurface->IsUPeriodic())
1833   {
1834     aDS.myPeriodicDir = 0;
1835     SplitOnDirection(aDS);
1836   }
1837   if (theSurface->IsVPeriodic())
1838   {
1839     aDS.myPeriodicDir = 1;
1840     SplitOnDirection(aDS);
1841   }
1842
1843   std::sort(aDS.mySplits.begin(), aDS.mySplits.end(), Comparator);
1844 }
1845
1846 //=======================================================================
1847 //function : SplitOnDirection
1848 //purpose  : This method compute points in the parameter space of the curve
1849 //           on which curve should be split since period jump is happen.
1850 //=======================================================================
1851 void SplitOnDirection(SplitDS & theSplitDS)
1852 {
1853   // Algorithm:
1854   // Create 3D curve which is correspond to the periodic bound in 2d space.
1855   // Run curve / curve extrema and run extrema point / surface to check that
1856   // the point will be projected to the periodic bound.
1857   // In this method assumed that the points cannot be closer to each other that 1% of the parameter space.
1858
1859   gp_Pnt2d aStartPnt(theSplitDS.mySurface->FirstUParameter(), theSplitDS.mySurface->FirstVParameter());
1860   gp_Dir2d aDir(theSplitDS.myPeriodicDir, (Standard_Integer)!theSplitDS.myPeriodicDir);
1861
1862   theSplitDS.myPerMinParam = !theSplitDS.myPeriodicDir ? theSplitDS.mySurface->FirstUParameter():
1863                                                          theSplitDS.mySurface->FirstVParameter();
1864   theSplitDS.myPerMaxParam = !theSplitDS.myPeriodicDir ? theSplitDS.mySurface->LastUParameter():
1865                                                          theSplitDS.mySurface->LastVParameter();
1866   Standard_Real aLast2DParam = theSplitDS.myPeriodicDir ?
1867                                theSplitDS.mySurface->LastUParameter() - theSplitDS.mySurface->FirstUParameter():
1868                                theSplitDS.mySurface->LastVParameter() - theSplitDS.mySurface->FirstVParameter();
1869
1870   // Create line which is represent periodic border.
1871   Handle(Geom2d_Curve) aC2GC = new Geom2d_Line(aStartPnt, aDir);
1872   Handle(Geom2dAdaptor_HCurve) aC = new Geom2dAdaptor_HCurve(aC2GC, 0, aLast2DParam);
1873   Adaptor3d_CurveOnSurface  aCOnS(aC, theSplitDS.mySurface);
1874
1875   Extrema_ExtCC anExtCC;
1876   anExtCC.SetCurve(1, aCOnS);
1877   anExtCC.SetCurve(2, theSplitDS.myCurve->Curve());
1878   anExtCC.SetSingleSolutionFlag(Standard_True); // Search only one solution since multiple invocations are needed.
1879   anExtCC.SetRange(1, 0, aLast2DParam);
1880   theSplitDS.myExtCC = &anExtCC;
1881
1882   FindSplitPoint(theSplitDS,
1883                  theSplitDS.myCurve->FirstParameter(), // Initial curve range.
1884                  theSplitDS.myCurve->LastParameter());
1885 }
1886
1887
1888 //=======================================================================
1889 //function : FindSplitPoint
1890 //purpose  :
1891 //=======================================================================
1892 void FindSplitPoint(SplitDS &theSplitDS,
1893                     const Standard_Real theMinParam,
1894                     const Standard_Real theMaxParam)
1895 {
1896   // Make extrema copy to avoid dependencies between different levels of the recursion.
1897   Extrema_ExtCC anExtCC(*theSplitDS.myExtCC);
1898   anExtCC.SetRange(2, theMinParam, theMaxParam);
1899   anExtCC.Perform();
1900
1901   if (anExtCC.IsDone())
1902   {
1903     const Standard_Integer aNbExt = anExtCC.NbExt();
1904     for (Standard_Integer anIdx = 1; anIdx <= aNbExt; ++anIdx)
1905     {
1906       Extrema_POnCurv aPOnC1, aPOnC2;
1907       anExtCC.Points(anIdx, aPOnC1, aPOnC2);
1908
1909       theSplitDS.myExtPS->Perform(aPOnC2.Value());
1910       if (!theSplitDS.myExtPS->IsDone())
1911         return;
1912
1913       // Find point with the minimal Euclidean distance to avoid
1914       // false positive points detection.
1915       Standard_Integer aMinIdx = -1;
1916       Standard_Real aMinSqDist = RealLast();
1917       const Standard_Integer aNbPext = theSplitDS.myExtPS->NbExt();
1918       for(Standard_Integer aPIdx = 1; aPIdx <= aNbPext; ++aPIdx)
1919       {
1920         const Standard_Real aCurrSqDist = theSplitDS.myExtPS->SquareDistance(aPIdx);
1921
1922         if (aCurrSqDist < aMinSqDist)
1923         {
1924           aMinSqDist = aCurrSqDist;
1925           aMinIdx = aPIdx;
1926         }
1927       }
1928
1929       // Check that is point will be projected to the periodic border.
1930       const Extrema_POnSurf &aPOnS = theSplitDS.myExtPS->Point(aMinIdx);
1931       Standard_Real U, V, aProjParam;
1932       aPOnS.Parameter(U, V);
1933       aProjParam = theSplitDS.myPeriodicDir ? V : U;
1934
1935
1936       if (Abs(aProjParam - theSplitDS.myPerMinParam) < Precision::PConfusion() ||
1937           Abs(aProjParam - theSplitDS.myPerMaxParam) < Precision::PConfusion() )
1938       {
1939         const Standard_Real aParam = aPOnC2.Parameter();
1940         const Standard_Real aCFParam = theSplitDS.myCurve->FirstParameter();
1941         const Standard_Real aCLParam = theSplitDS.myCurve->LastParameter();
1942
1943         if (aParam > aCFParam + Precision::PConfusion() &&
1944             aParam < aCLParam  - Precision::PConfusion() )
1945         {
1946           // Add only inner points.
1947           theSplitDS.mySplits.Append(aParam);
1948         }
1949
1950         const Standard_Real aDeltaCoeff = 0.01;
1951         const Standard_Real aDelta = (theMaxParam - theMinParam +
1952                                       aCLParam - aCFParam) * aDeltaCoeff;
1953
1954         if (aParam - aDelta > theMinParam + Precision::PConfusion())
1955         {
1956           FindSplitPoint(theSplitDS,
1957                          theMinParam, aParam - aDelta); // Curve parameters.
1958         }
1959
1960         if (aParam + aDelta < theMaxParam - Precision::PConfusion())
1961         {
1962           FindSplitPoint(theSplitDS,
1963                          aParam + aDelta, theMaxParam); // Curve parameters.
1964         }
1965       }
1966     } // for (Standard_Integer anIdx = 1; anIdx <= aNbExt; ++anIdx)
1967   }
1968 }