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