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()) Standard_ConstructionError::Raise();
 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()) Standard_ConstructionError::Raise();
 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()) Standard_ConstructionError::Raise();
 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()) Standard_ConstructionError::Raise();
 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     for(j = 1; j <= 4-i; j++)
 418       if(Seq(j).Y() < Seq(j+1).Y())
 419       {
 420         gp_Pnt swp;
 421         swp = Seq.Value(j+1);
 422         Seq.ChangeValue(j+1) = Seq.Value(j);
 423         Seq.ChangeValue(j) = swp;
 424       }
 425
 426       t = Sol.X();
 427       t1 = Min(Sol.X(), NotSol);
 428       t2 = Max(Sol.X(), NotSol);
 429
 430       Standard_Boolean isDone = Standard_False;
 431       while (!Seq.IsEmpty())
 432       {
 433         gp_Pnt P;
 434         P = Seq.Last();
 435         Seq.Remove(Seq.Length());
 436         ProjLib_PrjResolve aPrjPS(Curve->Curve(),
 437           Surface->Surface(),
 438           Standard_Integer(P.Z()));
 439         if(Standard_Integer(P.Z()) == 2)
 440         {
 441           aPrjPS.Perform(t, P.X(), V0, gp_Pnt2d(Tol, TolV),
 442             gp_Pnt2d(t1, Surface->FirstVParameter()),
 443             gp_Pnt2d(t2, Surface->LastVParameter()), FuncTol);
 444           if(!aPrjPS.IsDone()) continue;
 445           POnS = aPrjPS.Solution();
 446           Sol = gp_Pnt(POnS.X(), P.X(), POnS.Y());
 447           isDone = Standard_True;
 448           break;
 449         }
 450         else
 451         {
 452           aPrjPS.Perform(t, U0, P.X(), gp_Pnt2d(Tol, TolU),
 453             gp_Pnt2d(t1, Surface->FirstUParameter()),
 454             gp_Pnt2d(t2, Surface->LastUParameter()), FuncTol);
 455           if(!aPrjPS.IsDone()) continue;
 456           POnS = aPrjPS.Solution();
 457           Sol = gp_Pnt(POnS.X(), POnS.Y(), P.X());
 458           isDone = Standard_True;
 459           break;
 460         }
 461       }
 462
 463       return isDone;
 464 }
 465
 466 //=======================================================================
 467 //function : DichExactBound
 468 //purpose  : computes exact boundary point
 469 //=======================================================================
 470
 471 static void DichExactBound(gp_Pnt& Sol,
 472   const Standard_Real NotSol,
 473   const Standard_Real Tol,
 474   const Standard_Real TolU,
 475   const Standard_Real TolV,
 476   const Handle(Adaptor3d_HCurve)& Curve,
 477   const Handle(Adaptor3d_HSurface)& Surface)
 478 {
 479 #ifdef OCCT_DEBUG_CHRONO
 480   InitChron(chr_dicho_bound);
 481 #endif
 482
 483   Standard_Real U0, V0, t;
 484   gp_Pnt2d POnS;
 485   U0 = Sol.Y();
 486   V0 = Sol.Z();
 487   ProjLib_PrjResolve aPrjPS(Curve->Curve(), Surface->Surface(), 1);
 488
 489   Standard_Real aNotSol = NotSol;
 490   while (fabs(Sol.X() - aNotSol) > Tol)
 491   {
 492     t = (Sol.X() + aNotSol)/2;
 493     aPrjPS.Perform(t, U0, V0, gp_Pnt2d(TolU, TolV),
 494       gp_Pnt2d(Surface->FirstUParameter(),Surface->FirstVParameter()),
 495       gp_Pnt2d(Surface->LastUParameter(),Surface->LastVParameter()),
 496       FuncTol, Standard_True);
 497
 498     if (aPrjPS.IsDone())
 499     {
 500       POnS = aPrjPS.Solution();
 501       Sol = gp_Pnt(t, POnS.X(), POnS.Y());
 502       U0=Sol.Y();
 503       V0=Sol.Z();
 504     }
 505     else aNotSol = t;
 506   }
 507 #ifdef OCCT_DEBUG_CHRONO
 508   ResultChron(chr_dicho_bound,t_dicho_bound);
 509   dicho_bound_count++;
 510 #endif
 511 }
 512
 513 //=======================================================================
 514 //function : InitialPoint
 515 //purpose  :
 516 //=======================================================================
 517
 518 static Standard_Boolean InitialPoint(const gp_Pnt& Point,
 519   const Standard_Real t,
 520   const Handle(Adaptor3d_HCurve)& C,
 521   const Handle(Adaptor3d_HSurface)& S,
 522   const Standard_Real TolU,
 523   const Standard_Real TolV,
 524   Standard_Real& U,
 525   Standard_Real& V)
 526 {
 527
 528   ProjLib_PrjResolve aPrjPS(C->Curve(), S->Surface(), 1);
 529   Standard_Real ParU,ParV;
 530   Extrema_ExtPS aExtPS;
 531   aExtPS.Initialize(S->Surface(), S->FirstUParameter(),
 532     S->LastUParameter(), S->FirstVParameter(),
 533     S->LastVParameter(), TolU, TolV);
 534
 535   aExtPS.Perform(Point);
 536   Standard_Integer argmin = 0;
 537   if (aExtPS.IsDone() && aExtPS.NbExt())
 538   {
 539     Standard_Integer i, Nend;
 540     // Search for the nearest solution which is also a normal projection
 541     Nend = aExtPS.NbExt();
 542     for(i = 1; i <= Nend; i++)
 543     {
 544       Extrema_POnSurf POnS = aExtPS.Point(i);
 545       POnS.Parameter(ParU, ParV);
 546       aPrjPS.Perform(t, ParU, ParV, gp_Pnt2d(TolU, TolV),
 547         gp_Pnt2d(S->FirstUParameter(), S->FirstVParameter()),
 548         gp_Pnt2d(S->LastUParameter(), S->LastVParameter()),
 549         FuncTol, Standard_True);
 550       if(aPrjPS.IsDone() )
 551         if  (argmin == 0 || aExtPS.SquareDistance(i) < aExtPS.SquareDistance(argmin)) argmin = i;
 552     }
 553   }
 554   if( argmin == 0 ) return Standard_False;
 555   else
 556   {
 557     Extrema_POnSurf POnS = aExtPS.Point(argmin);
 558     POnS.Parameter(U, V);
 559     return Standard_True;
 560   }
 561 }
 562
 563 //=======================================================================
 564 //function : ProjLib_CompProjectedCurve
 565 //purpose  :
 566 //=======================================================================
 567
 568 ProjLib_CompProjectedCurve::ProjLib_CompProjectedCurve()
 569 : myNbCurves(0),
 570   myTolU    (0.0),
 571   myTolV    (0.0),
 572   myMaxDist (0.0)
 573 {
 574 }
 575
 576 //=======================================================================
 577 //function : ProjLib_CompProjectedCurve
 578 //purpose  :
 579 //=======================================================================
 580
 581 ProjLib_CompProjectedCurve::ProjLib_CompProjectedCurve
 582                            (const Handle(Adaptor3d_HSurface)& theSurface,
 583                             const Handle(Adaptor3d_HCurve)&   theCurve,
 584                             const Standard_Real               theTolU,
 585                             const Standard_Real               theTolV)
 586 : mySurface (theSurface),
 587   myCurve   (theCurve),
 588   myNbCurves(0),
 589   mySequence(new ProjLib_HSequenceOfHSequenceOfPnt()),
 590   myTolU    (theTolU),
 591   myTolV    (theTolV),
 592   myMaxDist (-1.0)
 593 {
 594   Init();
 595 }
 596
 597 //=======================================================================
 598 //function : ProjLib_CompProjectedCurve
 599 //purpose  :
 600 //=======================================================================
 601
 602 ProjLib_CompProjectedCurve::ProjLib_CompProjectedCurve
 603                            (const Handle(Adaptor3d_HSurface)& theSurface,
 604                             const Handle(Adaptor3d_HCurve)&   theCurve,
 605                             const Standard_Real               theTolU,
 606                             const Standard_Real               theTolV,
 607                             const Standard_Real               theMaxDist)
 608 : mySurface (theSurface),
 609   myCurve   (theCurve),
 610   myNbCurves(0),
 611   mySequence(new ProjLib_HSequenceOfHSequenceOfPnt()),
 612   myTolU    (theTolU),
 613   myTolV    (theTolV),
 614   myMaxDist (theMaxDist)
 615 {
 616   Init();
 617 }
 618
 619 //=======================================================================
 620 //function : Init
 621 //purpose  :
 622 //=======================================================================
 623
 624 void ProjLib_CompProjectedCurve::Init()
 625 {
 626   myTabInt.Nullify();
 627   NCollection_Vector<Standard_Real> aSplits;
 628   aSplits.Clear();
 629
 630   Standard_Real Tol;// Tolerance for ExactBound
 631   Standard_Integer i, Nend = 0, aSplitIdx = 0;
 632   Standard_Boolean FromLastU = Standard_False,
 633                    isSplitsComputed = Standard_False;
 634
 635   const Standard_Real aTol3D = Precision::Confusion();
 636   Extrema_ExtCS CExt(myCurve->Curve(), mySurface->Surface(), aTol3D, aTol3D);
 637   if (CExt.IsDone() && CExt.NbExt())
 638   {
 639     // Search for the minimum solution.
 640     // Avoid usage of extrema result that can be wrong for extrusion.
 641     if(myMaxDist > 0 &&
 642
 643        mySurface->GetType() != GeomAbs_SurfaceOfExtrusion)
 644     {
 645       Standard_Real min_val2;
 646       min_val2 = CExt.SquareDistance(1);
 647
 648       Nend = CExt.NbExt();
 649       for(i = 2; i <= Nend; i++)
 650       {
 651         if (CExt.SquareDistance(i) < min_val2)
 652           min_val2 = CExt.SquareDistance(i);
 653       }
 654       if (min_val2 > myMaxDist * myMaxDist)
 655         return; // No near solution -> exit.
 656     }
 657   }
 658
 659   Standard_Real FirstU, LastU, Step, SearchStep, WalkStep, t;
 660
 661   FirstU = myCurve->FirstParameter();
 662   LastU  = myCurve->LastParameter();
 663   const Standard_Real GlobalMinStep = 1.e-4;
 664   //<GlobalMinStep> is sufficiently small to provide solving from initial point
 665   //and, on the other hand, it is sufficiently large to avoid too close solutions.
 666   const Standard_Real MinStep = 0.01*(LastU - FirstU),
 667     MaxStep = 0.1*(LastU - FirstU);
 668   SearchStep = 10*MinStep;
 669   Step = SearchStep;
 670
 671   gp_Pnt2d aLowBorder(mySurface->FirstUParameter(),mySurface->FirstVParameter());
 672   gp_Pnt2d aUppBorder(mySurface->LastUParameter(), mySurface->LastVParameter());
 673   gp_Pnt2d aTol(myTolU, myTolV);
 674   ProjLib_PrjResolve aPrjPS(myCurve->Curve(), mySurface->Surface(), 1);
 675
 676   t = FirstU;
 677   Standard_Boolean new_part;
 678   Standard_Real prevDeb=0.;
 679   Standard_Boolean SameDeb=Standard_False;
 680
 681
 682   gp_Pnt Triple, prevTriple;
 683
 684   //Basic loop
 685   while(t <= LastU)
 686   {
 687     // Search for the beginning of a new continuous part
 688     // to avoid infinite computation in some difficult cases.
 689     new_part = Standard_False;
 690     if(t > FirstU && Abs(t-prevDeb) <= Precision::PConfusion()) SameDeb=Standard_True;
 691     while(t <= LastU && !new_part && !FromLastU && !SameDeb)
 692     {
 693       prevDeb=t;
 694       if (t == LastU) FromLastU=Standard_True;
 695       Standard_Boolean initpoint=Standard_False;
 696       Standard_Real U = 0., V = 0.;
 697       gp_Pnt CPoint;
 698       Standard_Real ParT,ParU,ParV;
 699
 700       // Search an initial point in the list of Extrema Curve-Surface
 701       if(Nend != 0 && !CExt.IsParallel())
 702       {
 703         for (i=1;i<=Nend;i++)
 704         {
 705           Extrema_POnCurv P1;
 706           Extrema_POnSurf P2;
 707           CExt.Points(i,P1,P2);
 708           ParT=P1.Parameter();
 709           P2.Parameter(ParU, ParV);
 710
 711           aPrjPS.Perform(ParT, ParU, ParV, aTol, aLowBorder, aUppBorder, FuncTol, Standard_True);
 712
 713           if ( aPrjPS.IsDone() && P1.Parameter() > Max(FirstU,t-Step+Precision::PConfusion())
 714             && P1.Parameter() <= t)
 715           {
 716             t=ParT;
 717             U=ParU;
 718             V=ParV;
 719             CPoint=P1.Value();
 720             initpoint = Standard_True;
 721             break;
 722           }
 723         }
 724       }
 725       if (!initpoint)
 726       {
 727         myCurve->D0(t,CPoint);
 728 #ifdef OCCT_DEBUG_CHRONO
 729         InitChron(chr_init_point);
 730 #endif
 731         // PConfusion - use geometric tolerances in extrema / optimization.
 732         initpoint=InitialPoint(CPoint, t,myCurve,mySurface, Precision::PConfusion(), Precision::PConfusion(), U, V);
 733 #ifdef OCCT_DEBUG_CHRONO
 734         ResultChron(chr_init_point,t_init_point);
 735         init_point_count++;
 736 #endif
 737       }
 738       if(initpoint)
 739       {
 740         // When U or V lie on surface joint in some cases we cannot use them
 741         // as initial point for aPrjPS, so we switch them
 742         gp_Vec2d D;
 743
 744         if ((mySurface->IsUPeriodic() &&
 745             Abs(aUppBorder.X() - aLowBorder.X() - mySurface->UPeriod()) < Precision::Confusion()) ||
 746             (mySurface->IsVPeriodic() &&
 747             Abs(aUppBorder.Y() - aLowBorder.Y() - mySurface->VPeriod()) < Precision::Confusion()))
 748         {
 749           if((Abs(U - aLowBorder.X()) < mySurface->UResolution(Precision::PConfusion())) &&
 750             mySurface->IsUPeriodic())
 751           {
 752             d1(t, U, V, D, myCurve, mySurface);
 753             if (D.X() < 0 ) U = aUppBorder.X();
 754           }
 755           else if((Abs(U - aUppBorder.X()) < mySurface->UResolution(Precision::PConfusion())) &&
 756             mySurface->IsUPeriodic())
 757           {
 758             d1(t, U, V, D, myCurve, mySurface);
 759             if (D.X() > 0) U = aLowBorder.X();
 760           }
 761
 762           if((Abs(V - aLowBorder.Y()) < mySurface->VResolution(Precision::PConfusion())) &&
 763             mySurface->IsVPeriodic())
 764           {
 765             d1(t, U, V, D, myCurve, mySurface);
 766             if (D.Y() < 0) V = aUppBorder.Y();
 767           }
 768           else if((Abs(V - aUppBorder.Y()) <= mySurface->VResolution(Precision::PConfusion())) &&
 769             mySurface->IsVPeriodic())
 770           {
 771             d1(t, U, V, D, myCurve, mySurface);
 772             if (D.Y() > 0) V = aLowBorder.Y();
 773           }
 774         }
 775
 776         if (myMaxDist > 0)
 777         {
 778           // Here we are going to stop if the distance between projection and
 779           // corresponding curve point is greater than myMaxDist
 780           gp_Pnt POnS;
 781           Standard_Real d;
 782           mySurface->D0(U, V, POnS);
 783           d = CPoint.Distance(POnS);
 784           if (d > myMaxDist)
 785           {
 786             mySequence->Clear();
 787             myNbCurves = 0;
 788             return;
 789           }
 790         }
 791         Triple = gp_Pnt(t, U, V);
 792         if (t != FirstU)
 793         {
 794           //Search for exact boundary point
 795           Tol = Min(myTolU, myTolV);
 796           gp_Vec2d aD;
 797           d1(Triple.X(), Triple.Y(), Triple.Z(), aD, myCurve, mySurface);
 798           Tol /= Max(Abs(aD.X()), Abs(aD.Y()));
 799
 800           if(!ExactBound(Triple, t - Step, Tol,
 801             myTolU, myTolV, myCurve, mySurface))
 802           {
 803 #ifdef OCCT_DEBUG
 804             cout<<"There is a problem with ExactBound computation"<<endl;
 805 #endif
 806             DichExactBound(Triple, t - Step, Tol, myTolU, myTolV,
 807               myCurve, mySurface);
 808           }
 809         }
 810         new_part = Standard_True;
 811       }
 812       else
 813       {
 814         if(t == LastU) break;
 815         t += Step;
 816         if(t>LastU)
 817         {
 818           Step =Step+LastU-t;
 819           t=LastU;
 820         }
 821       }
 822     }
 823     if (!new_part) break;
 824
 825     //We have found a new continuous part
 826     Handle(TColgp_HSequenceOfPnt) hSeq = new TColgp_HSequenceOfPnt();
 827     mySequence->Append(hSeq);
 828     myNbCurves++;
 829     mySequence->Value(myNbCurves)->Append(Triple);
 830     prevTriple = Triple;
 831
 832     if (Triple.X() == LastU) break;//return;
 833
 834     //Computation of WalkStep
 835     gp_Vec D1, D2;
 836     Standard_Real MagnD1, MagnD2;
 837     d2CurvOnSurf(Triple.X(), Triple.Y(), Triple.Z(), D1, D2, myCurve, mySurface);
 838     MagnD1 = D1.Magnitude();
 839     MagnD2 = D2.Magnitude();
 840     if(MagnD2 < Precision::Confusion()) WalkStep = MaxStep;
 841     else WalkStep = Min(MaxStep, Max(MinStep, 0.1*MagnD1/MagnD2));
 842
 843     Step = WalkStep;
 844
 845     t = Triple.X() + Step;
 846     if (t > LastU) t = LastU;
 847     Standard_Real prevStep = Step;
 848     Standard_Real U0, V0;
 849
 850     //Here we are trying to prolong continuous part
 851     while (t <= LastU && new_part)
 852     {
 853
 854       U0 = Triple.Y() + (Step / prevStep) * (Triple.Y() - prevTriple.Y());
 855       V0 = Triple.Z() + (Step / prevStep) * (Triple.Z() - prevTriple.Z());
 856       // adjust U0 to be in [mySurface->FirstUParameter(),mySurface->LastUParameter()]
 857       U0 = Min(Max(U0, aLowBorder.X()), aUppBorder.X());
 858       // adjust V0 to be in [mySurface->FirstVParameter(),mySurface->LastVParameter()]
 859       V0 = Min(Max(V0, aLowBorder.Y()), aUppBorder.Y());
 860
 861
 862       aPrjPS.Perform(t, U0, V0, aTol,
 863                      aLowBorder, aUppBorder, FuncTol, Standard_True);
 864       if(!aPrjPS.IsDone())
 865       {
 866         if (Step <= GlobalMinStep)
 867         {
 868           //Search for exact boundary point
 869           Tol = Min(myTolU, myTolV);
 870           gp_Vec2d D;
 871           d1(Triple.X(), Triple.Y(), Triple.Z(), D, myCurve, mySurface);
 872           Tol /= Max(Abs(D.X()), Abs(D.Y()));
 873
 874           if(!ExactBound(Triple, t, Tol, myTolU, myTolV,
 875             myCurve, mySurface))
 876           {
 877 #ifdef OCCT_DEBUG
 878             cout<<"There is a problem with ExactBound computation"<<endl;
 879 #endif
 880             DichExactBound(Triple, t, Tol, myTolU, myTolV,
 881               myCurve, mySurface);
 882           }
 883
 884           if((Triple.X() - mySequence->Value(myNbCurves)->Value(mySequence->Value(myNbCurves)->Length()).X()) > 1.e-10)
 885             mySequence->Value(myNbCurves)->Append(Triple);
 886           if((LastU - Triple.X()) < Tol) {t = LastU + 1; break;}//return;
 887
 888           Step = SearchStep;
 889           t = Triple.X() + Step;
 890           if (t > (LastU-MinStep/2) )
 891           {
 892             Step =Step+LastU-t;
 893             t = LastU;
 894           }
 895           new_part = Standard_False;
 896         }
 897         else
 898         {
 899           // decrease step
 900           Standard_Real SaveStep = Step;
 901           Step /= 2.;
 902           t = Triple .X() + Step;
 903           if (t > (LastU-MinStep/4) )
 904           {
 905             Step =Step+LastU-t;
 906             if (Abs(Step - SaveStep) <= Precision::PConfusion())
 907               Step = GlobalMinStep; //to avoid looping
 908             t = LastU;
 909           }
 910         }
 911       }
 912       // Go further
 913       else
 914       {
 915         prevTriple = Triple;
 916         prevStep = Step;
 917         Triple = gp_Pnt(t, aPrjPS.Solution().X(), aPrjPS.Solution().Y());
 918
 919         // Check for possible local traps.
 920         UpdateTripleByTrapCriteria(Triple);
 921
 922         // Protection from case when the whole curve lies on a seam.
 923         if (!isSplitsComputed)
 924         {
 925           Standard_Boolean isUPossible = Standard_False;
 926           if (mySurface->IsUPeriodic() &&
 927              (Abs(Triple.Y() - mySurface->FirstUParameter() ) > Precision::PConfusion() &&
 928               Abs(Triple.Y() - mySurface->LastUParameter()  ) > Precision::PConfusion()))
 929           {
 930             isUPossible = Standard_True;
 931           }
 932
 933           Standard_Boolean isVPossible = Standard_False;
 934           if (mySurface->IsVPeriodic() &&
 935              (Abs(Triple.Z() - mySurface->FirstVParameter() ) > Precision::PConfusion() &&
 936               Abs(Triple.Z() - mySurface->LastVParameter()  ) > Precision::PConfusion()))
 937           {
 938             isVPossible = Standard_True;
 939           }
 940
 941           if (isUPossible || isVPossible)
 942           {
 943             // When point is good conditioned.
 944             BuildCurveSplits(myCurve, mySurface, myTolU, myTolV, aSplits);
 945             isSplitsComputed = Standard_True;
 946           }
 947         }
 948
 949         if((Triple.X() - mySequence->Value(myNbCurves)->Value(mySequence->Value(myNbCurves)->Length()).X()) > 1.e-10)
 950           mySequence->Value(myNbCurves)->Append(Triple);
 951         if (t == LastU) {t = LastU + 1; break;}//return;
 952         //Computation of WalkStep
 953         d2CurvOnSurf(Triple.X(), Triple.Y(), Triple.Z(), D1, D2, myCurve, mySurface);
 954         MagnD1 = D1.Magnitude();
 955         MagnD2 = D2.Magnitude();
 956         if(MagnD2 < Precision::Confusion() ) WalkStep = MaxStep;
 957         else WalkStep = Min(MaxStep, Max(MinStep, 0.1*MagnD1/MagnD2));
 958
 959         Step = WalkStep;
 960         t += Step;
 961         if (t > (LastU-MinStep/2))
 962         {
 963           Step = Step + LastU - t;
 964           t = LastU;
 965         }
 966
 967         // We assume at least one point of cache inside of a split.
 968         const Standard_Integer aSize = aSplits.Size();
 969         for(Standard_Integer anIdx = aSplitIdx; anIdx < aSize; ++anIdx)
 970         {
 971           const Standard_Real aParam = aSplits(anIdx);
 972           if (Abs(aParam - Triple.X() ) < Precision::PConfusion())
 973           {
 974             // The current point is equal to a split point.
 975             new_part = Standard_False;
 976
 977             // Move split index to avoid check of the whole list.
 978             ++aSplitIdx;
 979             break;
 980           }
 981           else if (aParam < t + Precision::PConfusion() )
 982           {
 983             // The next point crosses the split point.
 984             t = aParam;
 985             Step = t - prevTriple.X();
 986           }
 987         } // for(Standard_Integer anIdx = aSplitIdx; anIdx < aSize; ++anIdx)
 988       }
 989     }
 990   }
 991
 992   // Sequence post-proceeding.
 993   Standard_Integer j;
 994
 995   // 1. Removing poor parts
 996   Standard_Integer NbPart=myNbCurves;
 997   Standard_Integer ipart=1;
 998   for(i = 1; i <= NbPart; i++) {
 999     //    Standard_Integer NbPoints = mySequence->Value(i)->Length();
1000     if(mySequence->Value(ipart)->Length() < 2) {
1001       mySequence->Remove(ipart);
1002       myNbCurves--;
1003     }
1004     else ipart++;
1005   }
1006
1007   if(myNbCurves == 0) return;
1008
1009   // 2. Removing common parts of bounds
1010   for(i = 1; i < myNbCurves; i++)
1011   {
1012     if(mySequence->Value(i)->Value(mySequence->Value(i)->Length()).X() >=
1013       mySequence->Value(i+1)->Value(1).X())
1014       mySequence->ChangeValue(i+1)->ChangeValue(1).SetX(mySequence->Value(i)->Value(mySequence->Value(i)->Length()).X() + 1.e-12);
1015   }
1016
1017   // 3. Computation of the maximum distance from each part of curve to surface
1018
1019   myMaxDistance = new TColStd_HArray1OfReal(1, myNbCurves);
1020   myMaxDistance->Init(0);
1021   for(i = 1; i <= myNbCurves; i++)
1022     for(j = 1; j <= mySequence->Value(i)->Length(); j++)
1023     {
1024       gp_Pnt POnC, POnS, aTriple;
1025       Standard_Real Distance;
1026       aTriple = mySequence->Value(i)->Value(j);
1027       myCurve->D0(aTriple.X(), POnC);
1028       mySurface->D0(aTriple.Y(), aTriple.Z(), POnS);
1029       Distance = POnC.Distance(POnS);
1030       if (myMaxDistance->Value(i) < Distance)
1031         myMaxDistance->ChangeValue(i) = Distance;
1032     }
1033
1034
1035     // 4. Check the projection to be a single point
1036
1037     gp_Pnt2d Pmoy, Pcurr, P;
1038     Standard_Real AveU, AveV;
1039     mySnglPnts = new TColStd_HArray1OfBoolean(1, myNbCurves);
1040     for(i = 1; i <= myNbCurves; i++) mySnglPnts->SetValue(i, Standard_True);
1041
1042     for(i = 1; i <= myNbCurves; i++)
1043     {
1044       //compute an average U and V
1045
1046       for(j = 1, AveU = 0., AveV = 0.; j <= mySequence->Value(i)->Length(); j++)
1047       {
1048         AveU += mySequence->Value(i)->Value(j).Y();
1049         AveV += mySequence->Value(i)->Value(j).Z();
1050       }
1051       AveU /= mySequence->Value(i)->Length();
1052       AveV /= mySequence->Value(i)->Length();
1053
1054       Pmoy.SetCoord(AveU,AveV);
1055       for(j = 1; j <= mySequence->Value(i)->Length(); j++)
1056       {
1057         Pcurr =
1058           gp_Pnt2d(mySequence->Value(i)->Value(j).Y(), mySequence->Value(i)->Value(j).Z());
1059         if (Pcurr.Distance(Pmoy) > ((myTolU < myTolV) ? myTolV : myTolU))
1060         {
1061           mySnglPnts->SetValue(i, Standard_False);
1062           break;
1063         }
1064       }
1065     }
1066
1067     // 5. Check the projection to be an isoparametric curve of the surface
1068
1069     myUIso = new TColStd_HArray1OfBoolean(1, myNbCurves);
1070     for(i = 1; i <= myNbCurves; i++) myUIso->SetValue(i, Standard_True);
1071
1072     myVIso = new TColStd_HArray1OfBoolean(1, myNbCurves);
1073     for(i = 1; i <= myNbCurves; i++) myVIso->SetValue(i, Standard_True);
1074
1075     for(i = 1; i <= myNbCurves; i++) {
1076       if (IsSinglePnt(i, P)|| mySequence->Value(i)->Length() <=2) {
1077         myUIso->SetValue(i, Standard_False);
1078         myVIso->SetValue(i, Standard_False);
1079         continue;
1080       }
1081
1082       // new test for isoparametrics
1083
1084       if ( mySequence->Value(i)->Length() > 2) {
1085         //compute an average U and V
1086
1087         for(j = 1, AveU = 0., AveV = 0.; j <= mySequence->Value(i)->Length(); j++) {
1088           AveU += mySequence->Value(i)->Value(j).Y();
1089           AveV += mySequence->Value(i)->Value(j).Z();
1090         }
1091         AveU /= mySequence->Value(i)->Length();
1092         AveV /= mySequence->Value(i)->Length();
1093
1094         // is i-part U-isoparametric ?
1095         for(j = 1; j <= mySequence->Value(i)->Length(); j++)
1096         {
1097           if(Abs(mySequence->Value(i)->Value(j).Y() - AveU) > myTolU)
1098           {
1099             myUIso->SetValue(i, Standard_False);
1100             break;
1101           }
1102         }
1103
1104         // is i-part V-isoparametric ?
1105         for(j = 1; j <= mySequence->Value(i)->Length(); j++)
1106         {
1107           if(Abs(mySequence->Value(i)->Value(j).Z() - AveV) > myTolV)
1108           {
1109             myVIso->SetValue(i, Standard_False);
1110             break;
1111           }
1112         }
1113         //
1114       }
1115     }
1116 }
1117 //=======================================================================
1118 //function : Load
1119 //purpose  :
1120 //=======================================================================
1121
1122 void ProjLib_CompProjectedCurve::Load(const Handle(Adaptor3d_HSurface)& S)
1123 {
1124   mySurface = S;
1125 }
1126
1127 //=======================================================================
1128 //function : Load
1129 //purpose  :
1130 //=======================================================================
1131
1132 void ProjLib_CompProjectedCurve::Load(const Handle(Adaptor3d_HCurve)& C)
1133 {
1134   myCurve = C;
1135 }
1136
1137 //=======================================================================
1138 //function : GetSurface
1139 //purpose  :
1140 //=======================================================================
1141
1142 const Handle(Adaptor3d_HSurface)& ProjLib_CompProjectedCurve::GetSurface() const
1143 {
1144   return mySurface;
1145 }
1146
1147
1148 //=======================================================================
1149 //function : GetCurve
1150 //purpose  :
1151 //=======================================================================
1152
1153 const Handle(Adaptor3d_HCurve)& ProjLib_CompProjectedCurve::GetCurve() const
1154 {
1155   return myCurve;
1156 }
1157
1158 //=======================================================================
1159 //function : GetTolerance
1160 //purpose  :
1161 //=======================================================================
1162
1163 void ProjLib_CompProjectedCurve::GetTolerance(Standard_Real& TolU,
1164   Standard_Real& TolV) const
1165 {
1166   TolU = myTolU;
1167   TolV = myTolV;
1168 }
1169
1170 //=======================================================================
1171 //function : NbCurves
1172 //purpose  :
1173 //=======================================================================
1174
1175 Standard_Integer ProjLib_CompProjectedCurve::NbCurves() const
1176 {
1177   return myNbCurves;
1178 }
1179 //=======================================================================
1180 //function : Bounds
1181 //purpose  :
1182 //=======================================================================
1183
1184 void ProjLib_CompProjectedCurve::Bounds(const Standard_Integer Index,
1185   Standard_Real& Udeb,
1186   Standard_Real& Ufin) const
1187 {
1188   if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1189   Udeb = mySequence->Value(Index)->Value(1).X();
1190   Ufin = mySequence->Value(Index)->Value(mySequence->Value(Index)->Length()).X();
1191 }
1192 //=======================================================================
1193 //function : IsSinglePnt
1194 //purpose  :
1195 //=======================================================================
1196
1197 Standard_Boolean ProjLib_CompProjectedCurve::IsSinglePnt(const Standard_Integer Index, gp_Pnt2d& P) const
1198 {
1199   if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1200   P = gp_Pnt2d(mySequence->Value(Index)->Value(1).Y(), mySequence->Value(Index)->Value(1).Z());
1201   return mySnglPnts->Value(Index);
1202 }
1203
1204 //=======================================================================
1205 //function : IsUIso
1206 //purpose  :
1207 //=======================================================================
1208
1209 Standard_Boolean ProjLib_CompProjectedCurve::IsUIso(const Standard_Integer Index, Standard_Real& U) const
1210 {
1211   if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1212   U = mySequence->Value(Index)->Value(1).Y();
1213   return myUIso->Value(Index);
1214 }
1215 //=======================================================================
1216 //function : IsVIso
1217 //purpose  :
1218 //=======================================================================
1219
1220 Standard_Boolean ProjLib_CompProjectedCurve::IsVIso(const Standard_Integer Index, Standard_Real& V) const
1221 {
1222   if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1223   V = mySequence->Value(Index)->Value(1).Z();
1224   return myVIso->Value(Index);
1225 }
1226 //=======================================================================
1227 //function : Value
1228 //purpose  :
1229 //=======================================================================
1230
1231 gp_Pnt2d ProjLib_CompProjectedCurve::Value(const Standard_Real t) const
1232 {
1233   gp_Pnt2d P;
1234   D0(t, P);
1235   return P;
1236 }
1237 //=======================================================================
1238 //function : D0
1239 //purpose  :
1240 //=======================================================================
1241
1242 void ProjLib_CompProjectedCurve::D0(const Standard_Real U,gp_Pnt2d& P) const
1243 {
1244   Standard_Integer i, j;
1245   Standard_Real Udeb, Ufin;
1246   Standard_Boolean found = Standard_False;
1247
1248   for(i = 1; i <= myNbCurves; i++)
1249   {
1250     Bounds(i, Udeb, Ufin);
1251     if (U >= Udeb && U <= Ufin)
1252     {
1253       found = Standard_True;
1254       break;
1255     }
1256   }
1257   if (!found) Standard_DomainError::Raise("ProjLib_CompProjectedCurve::D0");
1258
1259   Standard_Real U0, V0;
1260
1261   Standard_Integer End = mySequence->Value(i)->Length();
1262   for(j = 1; j < End; j++)
1263     if ((U >= mySequence->Value(i)->Value(j).X()) && (U <= mySequence->Value(i)->Value(j + 1).X())) break;
1264
1265   //  U0 = mySequence->Value(i)->Value(j).Y();
1266   //  V0 = mySequence->Value(i)->Value(j).Z();
1267
1268   //  Cubic Interpolation
1269   if(mySequence->Value(i)->Length() < 4 ||
1270     (Abs(U-mySequence->Value(i)->Value(j).X()) <= Precision::PConfusion()) )
1271   {
1272     U0 = mySequence->Value(i)->Value(j).Y();
1273     V0 = mySequence->Value(i)->Value(j).Z();
1274   }
1275   else if (Abs(U-mySequence->Value(i)->Value(j+1).X())
1276     <= Precision::PConfusion())
1277   {
1278     U0 = mySequence->Value(i)->Value(j+1).Y();
1279     V0 = mySequence->Value(i)->Value(j+1).Z();
1280   }
1281   else
1282   {
1283     if (j == 1) j = 2;
1284     if (j > mySequence->Value(i)->Length() - 2)
1285       j = mySequence->Value(i)->Length() - 2;
1286
1287     gp_Vec2d I1, I2, I3, I21, I22, I31, Y1, Y2, Y3, Y4, Res;
1288     Standard_Real X1, X2, X3, X4;
1289
1290     X1 = mySequence->Value(i)->Value(j - 1).X();
1291     X2 = mySequence->Value(i)->Value(j).X();
1292     X3 = mySequence->Value(i)->Value(j + 1).X();
1293     X4 = mySequence->Value(i)->Value(j + 2).X();
1294
1295     Y1 = gp_Vec2d(mySequence->Value(i)->Value(j - 1).Y(),
1296       mySequence->Value(i)->Value(j - 1).Z());
1297     Y2 = gp_Vec2d(mySequence->Value(i)->Value(j).Y(),
1298       mySequence->Value(i)->Value(j).Z());
1299     Y3 = gp_Vec2d(mySequence->Value(i)->Value(j + 1).Y(),
1300       mySequence->Value(i)->Value(j + 1).Z());
1301     Y4 = gp_Vec2d(mySequence->Value(i)->Value(j + 2).Y(),
1302       mySequence->Value(i)->Value(j + 2).Z());
1303
1304     I1 = (Y1 - Y2)/(X1 - X2);
1305     I2 = (Y2 - Y3)/(X2 - X3);
1306     I3 = (Y3 - Y4)/(X3 - X4);
1307
1308     I21 = (I1 - I2)/(X1 - X3);
1309     I22 = (I2 - I3)/(X2 - X4);
1310
1311     I31 = (I21 - I22)/(X1 - X4);
1312
1313     Res = Y1 + (U - X1)*(I1 + (U - X2)*(I21 + (U - X3)*I31));
1314
1315     U0 = Res.X();
1316     V0 = Res.Y();
1317
1318     if(U0 < mySurface->FirstUParameter()) U0 = mySurface->FirstUParameter();
1319     else if(U0 > mySurface->LastUParameter()) U0 = mySurface->LastUParameter();
1320
1321     if(V0 < mySurface->FirstVParameter()) V0 = mySurface->FirstVParameter();
1322     else if(V0 > mySurface->LastVParameter()) V0 = mySurface->LastVParameter();
1323   }
1324   //End of cubic interpolation
1325
1326   ProjLib_PrjResolve aPrjPS(myCurve->Curve(), mySurface->Surface(), 1);
1327   aPrjPS.Perform(U, U0, V0, gp_Pnt2d(myTolU, myTolV),
1328     gp_Pnt2d(mySurface->FirstUParameter(), mySurface->FirstVParameter()),
1329     gp_Pnt2d(mySurface->LastUParameter(), mySurface->LastVParameter()));
1330   if (aPrjPS.IsDone())
1331     P = aPrjPS.Solution();
1332   else
1333   {
1334     gp_Pnt thePoint = myCurve->Value(U);
1335     Extrema_ExtPS aExtPS(thePoint, mySurface->Surface(), myTolU, myTolV);
1336     if (aExtPS.IsDone() && aExtPS.NbExt())
1337     {
1338       Standard_Integer k, Nend, imin = 1;
1339       // Search for the nearest solution which is also a normal projection
1340       Nend = aExtPS.NbExt();
1341       for(k = 2; k <= Nend; k++)
1342         if (aExtPS.SquareDistance(k) < aExtPS.SquareDistance(imin))
1343           imin = k;
1344       const Extrema_POnSurf& POnS = aExtPS.Point(imin);
1345       Standard_Real ParU,ParV;
1346       POnS.Parameter(ParU, ParV);
1347       P.SetCoord(ParU, ParV);
1348     }
1349     else
1350       P.SetCoord(U0,V0);
1351   }
1352 }
1353 //=======================================================================
1354 //function : D1
1355 //purpose  :
1356 //=======================================================================
1357
1358 void ProjLib_CompProjectedCurve::D1(const Standard_Real t,
1359   gp_Pnt2d& P,
1360   gp_Vec2d& V) const
1361 {
1362   Standard_Real u, v;
1363   D0(t, P);
1364   u = P.X();
1365   v = P.Y();
1366   d1(t, u, v, V, myCurve, mySurface);
1367 }
1368 //=======================================================================
1369 //function : D2
1370 //purpose  :
1371 //=======================================================================
1372
1373 void ProjLib_CompProjectedCurve::D2(const Standard_Real t,
1374   gp_Pnt2d& P,
1375   gp_Vec2d& V1,
1376   gp_Vec2d& V2) const
1377 {
1378   Standard_Real u, v;
1379   D0(t, P);
1380   u = P.X();
1381   v = P.Y();
1382   d2(t, u, v, V1, V2, myCurve, mySurface);
1383 }
1384 //=======================================================================
1385 //function : DN
1386 //purpose  :
1387 //=======================================================================
1388
1389 gp_Vec2d ProjLib_CompProjectedCurve::DN(const Standard_Real t,
1390   const Standard_Integer N) const
1391 {
1392   if (N < 1 ) Standard_OutOfRange::Raise("ProjLib_CompProjectedCurve : N must be greater than 0");
1393   else if (N ==1)
1394   {
1395     gp_Pnt2d P;
1396     gp_Vec2d V;
1397     D1(t,P,V);
1398     return V;
1399   }
1400   else if ( N==2)
1401   {
1402     gp_Pnt2d P;
1403     gp_Vec2d V1,V2;
1404     D2(t,P,V1,V2);
1405     return V2;
1406   }
1407   else if (N > 2 )
1408     Standard_NotImplemented::Raise("ProjLib_CompProjectedCurve::DN");
1409   return gp_Vec2d();
1410 }
1411
1412 //=======================================================================
1413 //function : GetSequence
1414 //purpose  :
1415 //=======================================================================
1416
1417 const Handle(ProjLib_HSequenceOfHSequenceOfPnt)& ProjLib_CompProjectedCurve::GetSequence() const
1418 {
1419   return mySequence;
1420 }
1421 //=======================================================================
1422 //function : FirstParameter
1423 //purpose  :
1424 //=======================================================================
1425
1426 Standard_Real ProjLib_CompProjectedCurve::FirstParameter() const
1427 {
1428   return myCurve->FirstParameter();
1429 }
1430
1431 //=======================================================================
1432 //function : LastParameter
1433 //purpose  :
1434 //=======================================================================
1435
1436 Standard_Real ProjLib_CompProjectedCurve::LastParameter() const
1437 {
1438   return myCurve->LastParameter();
1439 }
1440
1441 //=======================================================================
1442 //function : MaxDistance
1443 //purpose  :
1444 //=======================================================================
1445
1446 Standard_Real ProjLib_CompProjectedCurve::MaxDistance(const Standard_Integer Index) const
1447 {
1448   if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1449   return myMaxDistance->Value(Index);
1450 }
1451
1452 //=======================================================================
1453 //function : NbIntervals
1454 //purpose  :
1455 //=======================================================================
1456
1457 Standard_Integer ProjLib_CompProjectedCurve::NbIntervals(const GeomAbs_Shape S) const
1458 {
1459   const_cast<ProjLib_CompProjectedCurve*>(this)->myTabInt.Nullify();
1460   BuildIntervals(S);
1461   return myTabInt->Length() - 1;
1462 }
1463
1464 //=======================================================================
1465 //function : Intervals
1466 //purpose  :
1467 //=======================================================================
1468
1469 void ProjLib_CompProjectedCurve::Intervals(TColStd_Array1OfReal& T,const GeomAbs_Shape S) const
1470 {
1471   if (myTabInt.IsNull()) BuildIntervals (S);
1472   T = myTabInt->Array1();
1473 }
1474
1475 //=======================================================================
1476 //function : BuildIntervals
1477 //purpose  :
1478 //=======================================================================
1479
1480 void ProjLib_CompProjectedCurve::BuildIntervals(const GeomAbs_Shape S) const
1481 {
1482   GeomAbs_Shape SforS = GeomAbs_CN;
1483   switch(S) {
1484   case GeomAbs_C0:
1485     SforS = GeomAbs_C1;
1486     break;
1487   case GeomAbs_C1:
1488     SforS = GeomAbs_C2;
1489     break;
1490   case GeomAbs_C2:
1491     SforS = GeomAbs_C3;
1492     break;
1493   case GeomAbs_C3:
1494     SforS = GeomAbs_CN;
1495     break;
1496   case GeomAbs_CN:
1497     SforS = GeomAbs_CN;
1498     break;
1499   default:
1500     Standard_OutOfRange::Raise();
1501   }
1502   Standard_Integer i, j, k;
1503   Standard_Integer NbIntCur = myCurve->NbIntervals(S);
1504   Standard_Integer NbIntSurU = mySurface->NbUIntervals(SforS);
1505   Standard_Integer NbIntSurV = mySurface->NbVIntervals(SforS);
1506
1507   TColStd_Array1OfReal CutPntsT(1, NbIntCur+1);
1508   TColStd_Array1OfReal CutPntsU(1, NbIntSurU+1);
1509   TColStd_Array1OfReal CutPntsV(1, NbIntSurV+1);
1510
1511   myCurve->Intervals(CutPntsT, S);
1512   mySurface->UIntervals(CutPntsU, SforS);
1513   mySurface->VIntervals(CutPntsV, SforS);
1514
1515   Standard_Real Tl, Tr, Ul, Ur, Vl, Vr, Tol;
1516
1517   Handle(TColStd_HArray1OfReal) BArr = NULL,
1518     CArr = NULL,
1519     UArr = NULL,
1520     VArr = NULL;
1521
1522   // proccessing projection bounds
1523   BArr = new TColStd_HArray1OfReal(1, 2*myNbCurves);
1524   for(i = 1; i <= myNbCurves; i++)
1525     Bounds(i, BArr->ChangeValue(2*i - 1), BArr->ChangeValue(2*i));
1526
1527   // proccessing curve discontinuities
1528   if(NbIntCur > 1) {
1529     CArr = new TColStd_HArray1OfReal(1, NbIntCur - 1);
1530     for(i = 1; i <= CArr->Length(); i++)
1531       CArr->ChangeValue(i) = CutPntsT(i + 1);
1532   }
1533
1534   // proccessing U-surface discontinuities
1535   TColStd_SequenceOfReal TUdisc;
1536
1537   for(k = 2; k <= NbIntSurU; k++) {
1538     //    cout<<"CutPntsU("<<k<<") = "<<CutPntsU(k)<<endl;
1539     for(i = 1; i <= myNbCurves; i++)
1540       for(j = 1; j < mySequence->Value(i)->Length(); j++) {
1541         Ul = mySequence->Value(i)->Value(j).Y();
1542         Ur = mySequence->Value(i)->Value(j + 1).Y();
1543
1544         if(Abs(Ul - CutPntsU(k)) <= myTolU)
1545           TUdisc.Append(mySequence->Value(i)->Value(j).X());
1546         else if(Abs(Ur - CutPntsU(k)) <= myTolU)
1547           TUdisc.Append(mySequence->Value(i)->Value(j + 1).X());
1548         else if((Ul < CutPntsU(k) && CutPntsU(k) < Ur) ||
1549           (Ur < CutPntsU(k) && CutPntsU(k) < Ul))
1550         {
1551           Standard_Real V;
1552           V = (mySequence->Value(i)->Value(j).Z()
1553             + mySequence->Value(i)->Value(j +1).Z())/2;
1554           ProjLib_PrjResolve Solver(myCurve->Curve(), mySurface->Surface(), 2);
1555
1556           gp_Vec2d D;
1557           gp_Pnt Triple;
1558           Triple = mySequence->Value(i)->Value(j);
1559           d1(Triple.X(), Triple.Y(), Triple.Z(), D, myCurve, mySurface);
1560           if (Abs(D.X()) < Precision::Confusion())
1561             Tol = myTolU;
1562           else
1563             Tol = Min(myTolU, myTolU / Abs(D.X()));
1564
1565           Tl = mySequence->Value(i)->Value(j).X();
1566           Tr = mySequence->Value(i)->Value(j + 1).X();
1567
1568           Solver.Perform((Tl + Tr)/2, CutPntsU(k), V,
1569             gp_Pnt2d(Tol, myTolV),
1570             gp_Pnt2d(Tl, mySurface->FirstVParameter()),
1571             gp_Pnt2d(Tr, mySurface->LastVParameter()));
1572           //
1573           if(Solver.IsDone())
1574           {
1575             TUdisc.Append(Solver.Solution().X());
1576           }
1577         }
1578       }
1579   }
1580   for(i = 2; i <= TUdisc.Length(); i++)
1581     if(TUdisc(i) - TUdisc(i-1) < Precision::PConfusion())
1582       TUdisc.Remove(i--);
1583
1584   if(TUdisc.Length())
1585   {
1586     UArr = new TColStd_HArray1OfReal(1, TUdisc.Length());
1587     for(i = 1; i <= UArr->Length(); i++)
1588       UArr->ChangeValue(i) = TUdisc(i);
1589   }
1590   // proccessing V-surface discontinuities
1591   TColStd_SequenceOfReal TVdisc;
1592
1593   for(k = 2; k <= NbIntSurV; k++)
1594     for(i = 1; i <= myNbCurves; i++)
1595     {
1596       //      cout<<"CutPntsV("<<k<<") = "<<CutPntsV(k)<<endl;
1597       for(j = 1; j < mySequence->Value(i)->Length(); j++) {
1598
1599         Vl = mySequence->Value(i)->Value(j).Z();
1600         Vr = mySequence->Value(i)->Value(j + 1).Z();
1601
1602         if(Abs(Vl - CutPntsV(k)) <= myTolV)
1603           TVdisc.Append(mySequence->Value(i)->Value(j).X());
1604         else if (Abs(Vr - CutPntsV(k)) <= myTolV)
1605           TVdisc.Append(mySequence->Value(i)->Value(j + 1).X());
1606         else if((Vl < CutPntsV(k) && CutPntsV(k) < Vr) ||
1607           (Vr < CutPntsV(k) && CutPntsV(k) < Vl))
1608         {
1609           Standard_Real U;
1610           U = (mySequence->Value(i)->Value(j).Y()
1611             + mySequence->Value(i)->Value(j +1).Y())/2;
1612           ProjLib_PrjResolve Solver(myCurve->Curve(), mySurface->Surface(), 3);
1613
1614           gp_Vec2d D;
1615           gp_Pnt Triple;
1616           Triple = mySequence->Value(i)->Value(j);
1617           d1(Triple.X(), Triple.Y(), Triple.Z(), D, myCurve, mySurface);
1618           if (Abs(D.Y()) < Precision::Confusion())
1619             Tol = myTolV;
1620           else
1621             Tol = Min(myTolV, myTolV / Abs(D.Y()));
1622
1623           Tl = mySequence->Value(i)->Value(j).X();
1624           Tr = mySequence->Value(i)->Value(j + 1).X();
1625
1626           Solver.Perform((Tl + Tr)/2, U, CutPntsV(k),
1627             gp_Pnt2d(Tol, myTolV),
1628             gp_Pnt2d(Tl, mySurface->FirstUParameter()),
1629             gp_Pnt2d(Tr, mySurface->LastUParameter()));
1630           //
1631           if(Solver.IsDone())
1632           {
1633             TVdisc.Append(Solver.Solution().X());
1634           }
1635         }
1636       }
1637     }
1638     for(i = 2; i <= TVdisc.Length(); i++)
1639       if(TVdisc(i) - TVdisc(i-1) < Precision::PConfusion())
1640         TVdisc.Remove(i--);
1641
1642     if(TVdisc.Length())
1643     {
1644       VArr = new TColStd_HArray1OfReal(1, TVdisc.Length());
1645       for(i = 1; i <= VArr->Length(); i++)
1646         VArr->ChangeValue(i) = TVdisc(i);
1647     }
1648
1649     // fusion
1650     TColStd_SequenceOfReal Fusion;
1651     if(!CArr.IsNull())
1652     {
1653       GeomLib::FuseIntervals(BArr->ChangeArray1(),
1654         CArr->ChangeArray1(),
1655         Fusion, Precision::PConfusion());
1656       BArr = new TColStd_HArray1OfReal(1, Fusion.Length());
1657       for(i = 1; i <= BArr->Length(); i++)
1658         BArr->ChangeValue(i) = Fusion(i);
1659       Fusion.Clear();
1660     }
1661
1662     if(!UArr.IsNull())
1663     {
1664       GeomLib::FuseIntervals(BArr->ChangeArray1(),
1665         UArr->ChangeArray1(),
1666         Fusion, Precision::PConfusion());
1667       BArr = new TColStd_HArray1OfReal(1, Fusion.Length());
1668       for(i = 1; i <= BArr->Length(); i++)
1669         BArr->ChangeValue(i) = Fusion(i);
1670       Fusion.Clear();
1671     }
1672
1673     if(!VArr.IsNull())
1674     {
1675       GeomLib::FuseIntervals(BArr->ChangeArray1(),
1676         VArr->ChangeArray1(),
1677         Fusion, Precision::PConfusion());
1678       BArr = new TColStd_HArray1OfReal(1, Fusion.Length());
1679       for(i = 1; i <= BArr->Length(); i++)
1680         BArr->ChangeValue(i) = Fusion(i);
1681     }
1682
1683     const_cast<ProjLib_CompProjectedCurve*>(this)->myTabInt = new TColStd_HArray1OfReal(1, BArr->Length());
1684     for(i = 1; i <= BArr->Length(); i++)
1685       myTabInt->ChangeValue(i) = BArr->Value(i);
1686
1687 }
1688
1689 //=======================================================================
1690 //function : Trim
1691 //purpose  :
1692 //=======================================================================
1693
1694 Handle(Adaptor2d_HCurve2d) ProjLib_CompProjectedCurve::Trim
1695   (const Standard_Real First,
1696   const Standard_Real Last,
1697   const Standard_Real Tol) const
1698 {
1699   Handle(ProjLib_HCompProjectedCurve) HCS =
1700     new ProjLib_HCompProjectedCurve(*this);
1701   HCS->ChangeCurve2d().Load(mySurface);
1702   HCS->ChangeCurve2d().Load(myCurve->Trim(First,Last,Tol));
1703   return HCS;
1704 }
1705
1706 //=======================================================================
1707 //function : GetType
1708 //purpose  :
1709 //=======================================================================
1710
1711 GeomAbs_CurveType ProjLib_CompProjectedCurve::GetType() const
1712 {
1713   return GeomAbs_OtherCurve;
1714 }
1715
1716 //=======================================================================
1717 //function : UpdateTripleByTrapCriteria
1718 //purpose  :
1719 //=======================================================================
1720 void ProjLib_CompProjectedCurve::UpdateTripleByTrapCriteria(gp_Pnt &thePoint) const
1721 {
1722   Standard_Boolean isProblemsPossible = Standard_False;
1723   // Check possible traps cases:
1724
1725   // 25892 bug.
1726   if (mySurface->GetType() == GeomAbs_SurfaceOfRevolution)
1727   {
1728     // Compute maximal deviation from 3D and choose the biggest one.
1729     Standard_Real aVRes = mySurface->VResolution(Precision::Confusion());
1730     Standard_Real aMaxTol = Max(Precision::PConfusion(), aVRes);
1731
1732     if (Abs (thePoint.Z() - mySurface->FirstVParameter()) < aMaxTol ||
1733         Abs (thePoint.Z() - mySurface->LastVParameter() ) < aMaxTol )
1734     {
1735       isProblemsPossible = Standard_True;
1736     }
1737   }
1738
1739   // 27135 bug. Trap on degenerated edge.
1740   if (mySurface->GetType() == GeomAbs_Sphere &&
1741      (Abs (thePoint.Z() - mySurface->FirstVParameter()) < Precision::PConfusion() ||
1742       Abs (thePoint.Z() - mySurface->LastVParameter() ) < Precision::PConfusion() ||
1743       Abs (thePoint.Y() - mySurface->FirstUParameter()) < Precision::PConfusion() ||
1744       Abs (thePoint.Y() - mySurface->LastUParameter() ) < Precision::PConfusion() ))
1745   {
1746     isProblemsPossible = Standard_True;
1747   }
1748
1749   if (!isProblemsPossible)
1750     return;
1751
1752   Standard_Real U,V;
1753   Standard_Boolean isDone =
1754     InitialPoint(myCurve->Value(thePoint.X()), thePoint.X(), myCurve, mySurface,
1755                  Precision::PConfusion(), Precision::PConfusion(), U, V);
1756
1757   if (!isDone)
1758     return;
1759
1760   // Restore original position in case of period jump.
1761   if (mySurface->IsUPeriodic() &&
1762       Abs (Abs(U - thePoint.Y()) - mySurface->UPeriod()) < Precision::PConfusion())
1763   {
1764     U = thePoint.Y();
1765   }
1766   if (mySurface->IsVPeriodic() &&
1767       Abs (Abs(V - thePoint.Z()) - mySurface->VPeriod()) < Precision::PConfusion())
1768   {
1769     V = thePoint.Z();
1770   }
1771   thePoint.SetY(U);
1772   thePoint.SetZ(V);
1773 }
1774
1775 //=======================================================================
1776 //function : BuildCurveSplits
1777 //purpose  :
1778 //=======================================================================
1779 void BuildCurveSplits(const Handle(Adaptor3d_HCurve)   &theCurve,
1780                       const Handle(Adaptor3d_HSurface) &theSurface,
1781                       const Standard_Real theTolU,
1782                       const Standard_Real theTolV,
1783                       NCollection_Vector<Standard_Real> &theSplits)
1784 {
1785   SplitDS aDS(theCurve, theSurface, theSplits);
1786
1787   Extrema_ExtPS anExtPS;
1788   anExtPS.Initialize(theSurface->Surface(),
1789                      theSurface->FirstUParameter(), theSurface->LastUParameter(),
1790                      theSurface->FirstVParameter(), theSurface->LastVParameter(),
1791                      theTolU, theTolV);
1792   aDS.myExtPS = &anExtPS;
1793
1794   if (theSurface->IsUPeriodic())
1795   {
1796     aDS.myPeriodicDir = 0;
1797     SplitOnDirection(aDS);
1798   }
1799   if (theSurface->IsVPeriodic())
1800   {
1801     aDS.myPeriodicDir = 1;
1802     SplitOnDirection(aDS);
1803   }
1804
1805   std::sort(aDS.mySplits.begin(), aDS.mySplits.end(), Comparator);
1806 }
1807
1808 //=======================================================================
1809 //function : SplitOnDirection
1810 //purpose  : This method compute points in the parameter space of the curve
1811 //           on which curve should be split since period jump is happen.
1812 //=======================================================================
1813 void SplitOnDirection(SplitDS & theSplitDS)
1814 {
1815   // Algorithm:
1816   // Create 3D curve which is correspond to the periodic bound in 2d space.
1817   // Run curve / curve extrema and run extrema point / surface to check that
1818   // the point will be projected to the periodic bound.
1819   // In this method assumed that the points cannot be closer to each other that 1% of the parameter space.
1820
1821   gp_Pnt2d aStartPnt(theSplitDS.mySurface->FirstUParameter(), theSplitDS.mySurface->FirstVParameter());
1822   gp_Dir2d aDir(theSplitDS.myPeriodicDir, (Standard_Integer)!theSplitDS.myPeriodicDir);
1823
1824   theSplitDS.myPerMinParam = !theSplitDS.myPeriodicDir ? theSplitDS.mySurface->FirstUParameter():
1825                                                          theSplitDS.mySurface->FirstVParameter();
1826   theSplitDS.myPerMaxParam = !theSplitDS.myPeriodicDir ? theSplitDS.mySurface->LastUParameter():
1827                                                          theSplitDS.mySurface->LastVParameter();
1828   Standard_Real aLast2DParam = theSplitDS.myPeriodicDir ?
1829                                theSplitDS.mySurface->LastUParameter() - theSplitDS.mySurface->FirstUParameter():
1830                                theSplitDS.mySurface->LastVParameter() - theSplitDS.mySurface->FirstVParameter();
1831
1832   // Create line which is represent periodic border.
1833   Handle(Geom2d_Curve) aC2GC = new Geom2d_Line(aStartPnt, aDir);
1834   Handle(Geom2dAdaptor_HCurve) aC = new Geom2dAdaptor_HCurve(aC2GC, 0, aLast2DParam);
1835   Adaptor3d_CurveOnSurface  aCOnS(aC, theSplitDS.mySurface);
1836
1837   Extrema_ExtCC anExtCC;
1838   anExtCC.SetCurve(1, aCOnS);
1839   anExtCC.SetCurve(2, theSplitDS.myCurve->Curve());
1840   anExtCC.SetSingleSolutionFlag(Standard_True); // Search only one solution since multiple invocations are needed.
1841   anExtCC.SetRange(1, 0, aLast2DParam);
1842   theSplitDS.myExtCC = &anExtCC;
1843
1844   FindSplitPoint(theSplitDS,
1845                  theSplitDS.myCurve->FirstParameter(), // Initial curve range.
1846                  theSplitDS.myCurve->LastParameter());
1847 }
1848
1849
1850 //=======================================================================
1851 //function : FindSplitPoint
1852 //purpose  :
1853 //=======================================================================
1854 void FindSplitPoint(SplitDS &theSplitDS,
1855                     const Standard_Real theMinParam,
1856                     const Standard_Real theMaxParam)
1857 {
1858   // Make extrema copy to avoid dependencies between different levels of the recursion.
1859   Extrema_ExtCC anExtCC(*theSplitDS.myExtCC);
1860   anExtCC.SetRange(2, theMinParam, theMaxParam);
1861   anExtCC.Perform();
1862
1863   if (anExtCC.IsDone())
1864   {
1865     const Standard_Integer aNbExt = anExtCC.NbExt();
1866     for (Standard_Integer anIdx = 1; anIdx <= aNbExt; ++anIdx)
1867     {
1868       Extrema_POnCurv aPOnC1, aPOnC2;
1869       anExtCC.Points(anIdx, aPOnC1, aPOnC2);
1870
1871       theSplitDS.myExtPS->Perform(aPOnC2.Value());
1872       if (!theSplitDS.myExtPS->IsDone())
1873         return;
1874
1875       // Find point with the minimal Euclidean distance to avoid
1876       // false positive points detection.
1877       Standard_Integer aMinIdx = -1;
1878       Standard_Real aMinSqDist = RealLast();
1879       const Standard_Integer aNbPext = theSplitDS.myExtPS->NbExt();
1880       for(Standard_Integer aPIdx = 1; aPIdx <= aNbPext; ++aPIdx)
1881       {
1882         const Standard_Real aCurrSqDist = theSplitDS.myExtPS->SquareDistance(aPIdx);
1883
1884         if (aCurrSqDist < aMinSqDist)
1885         {
1886           aMinSqDist = aCurrSqDist;
1887           aMinIdx = aPIdx;
1888         }
1889       }
1890
1891       // Check that is point will be projected to the periodic border.
1892       const Extrema_POnSurf &aPOnS = theSplitDS.myExtPS->Point(aMinIdx);
1893       Standard_Real U, V, aProjParam;
1894       aPOnS.Parameter(U, V);
1895       aProjParam = theSplitDS.myPeriodicDir ? V : U;
1896
1897
1898       if (Abs(aProjParam - theSplitDS.myPerMinParam) < Precision::PConfusion() ||
1899           Abs(aProjParam - theSplitDS.myPerMaxParam) < Precision::PConfusion() )
1900       {
1901         const Standard_Real aParam = aPOnC2.Parameter();
1902         const Standard_Real aCFParam = theSplitDS.myCurve->FirstParameter();
1903         const Standard_Real aCLParam = theSplitDS.myCurve->LastParameter();
1904
1905         if (aParam > aCFParam + Precision::PConfusion() &&
1906             aParam < aCLParam  - Precision::PConfusion() )
1907         {
1908           // Add only inner points.
1909           theSplitDS.mySplits.Append(aParam);
1910         }
1911
1912         const Standard_Real aDeltaCoeff = 0.01;
1913         const Standard_Real aDelta = (theMaxParam - theMinParam +
1914                                       aCLParam - aCFParam) * aDeltaCoeff;
1915
1916         if (aParam - aDelta > theMinParam + Precision::PConfusion())
1917         {
1918           FindSplitPoint(theSplitDS,
1919                          theMinParam, aParam - aDelta); // Curve parameters.
1920         }
1921
1922         if (aParam + aDelta < theMaxParam - Precision::PConfusion())
1923         {
1924           FindSplitPoint(theSplitDS,
1925                          aParam + aDelta, theMaxParam); // Curve parameters.
1926         }
1927       }
1928     } // for (Standard_Integer anIdx = 1; anIdx <= aNbExt; ++anIdx)
1929   }
1930 }