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