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 <Adaptor2d_HCurve2d.hxx>
  19 #include <Adaptor3d_HCurve.hxx>
  20 #include <Adaptor3d_HSurface.hxx>
  21 #include <Extrema_ExtCS.hxx>
  22 #include <Extrema_ExtPS.hxx>
  23 #include <Extrema_GenLocateExtPS.hxx>
  24 #include <Extrema_POnCurv.hxx>
  25 #include <Extrema_POnSurf.hxx>
  26 #include <GeomAbs_CurveType.hxx>
  27 #include <GeomLib.hxx>
  28 #include <gp_Mat2d.hxx>
  29 #include <gp_Pnt2d.hxx>
  30 #include <gp_Vec2d.hxx>
  31 #include <gp_XY.hxx>
  32 #include <Precision.hxx>
  33 #include <ProjLib_CompProjectedCurve.hxx>
  34 #include <ProjLib_HCompProjectedCurve.hxx>
  35 #include <ProjLib_PrjResolve.hxx>
  36 #include <Standard_DomainError.hxx>
  37 #include <Standard_NoSuchObject.hxx>
  38 #include <Standard_NotImplemented.hxx>
  39 #include <Standard_OutOfRange.hxx>
  40 #include <TColgp_HSequenceOfPnt.hxx>
  41
  42 #define FuncTol 1.e-10
  43
  44 #ifdef OCCT_DEBUG_CHRONO
  45 #include <OSD_Timer.hxx>
  46
  47 static OSD_Chronometer chr_init_point, chr_dicho_bound;
  48
  49 Standard_EXPORT Standard_Real t_init_point, t_dicho_bound;
  50 Standard_EXPORT Standard_Integer init_point_count, dicho_bound_count;
  51
  52 static void InitChron(OSD_Chronometer& ch)
  53 {
  54   ch.Reset();
  55   ch.Start();
  56 }
  57
  58 static void ResultChron( OSD_Chronometer & ch, Standard_Real & time)
  59 {
  60   Standard_Real tch ;
  61   ch.Stop();
  62   ch.Show(tch);
  63   time=time +tch;
  64 }
  65 #endif
  66
  67
  68 //=======================================================================
  69 //function : d1
  70 //purpose  : computes first derivative of the projected curve
  71 //=======================================================================
  72
  73 static void d1(const Standard_Real t,
  74   const Standard_Real u,
  75   const Standard_Real v,
  76   gp_Vec2d& V,
  77   const Handle(Adaptor3d_HCurve)& Curve,
  78   const Handle(Adaptor3d_HSurface)& Surface)
  79 {
  80   gp_Pnt S, C;
  81   gp_Vec DS1_u, DS1_v, DS2_u, DS2_uv, DS2_v, DC1_t;
  82   Surface->D2(u, v, S, DS1_u, DS1_v, DS2_u, DS2_v, DS2_uv);
  83   Curve->D1(t, C, DC1_t);
  84   gp_Vec Ort(C, S);// Ort = S - C
  85
  86   gp_Vec2d dE_dt(-DC1_t*DS1_u, -DC1_t*DS1_v);
  87   gp_XY dE_du(DS1_u*DS1_u + Ort*DS2_u,
  88     DS1_u*DS1_v + Ort*DS2_uv);
  89   gp_XY dE_dv(DS1_v*DS1_u + Ort*DS2_uv,
  90     DS1_v*DS1_v + Ort*DS2_v);
  91
  92   Standard_Real det = dE_du.X()*dE_dv.Y() - dE_du.Y()*dE_dv.X();
  93   if (fabs(det) < gp::Resolution()) Standard_ConstructionError::Raise();
  94
  95   gp_Mat2d M(gp_XY(dE_dv.Y()/det, -dE_du.Y()/det),
  96     gp_XY(-dE_dv.X()/det, dE_du.X()/det));
  97
  98   V = - gp_Vec2d(gp_Vec2d(M.Row(1))*dE_dt, gp_Vec2d(M.Row(2))*dE_dt);
  99 }
 100
 101 //=======================================================================
 102 //function : d2
 103 //purpose  : computes second derivative of the projected curve
 104 //=======================================================================
 105
 106 static void d2(const Standard_Real t,
 107   const Standard_Real u,
 108   const Standard_Real v,
 109   gp_Vec2d& V1, gp_Vec2d& V2,
 110   const Handle(Adaptor3d_HCurve)& Curve,
 111   const Handle(Adaptor3d_HSurface)& Surface)
 112 {
 113   gp_Pnt S, C;
 114   gp_Vec DS1_u, DS1_v, DS2_u, DS2_uv, DS2_v,
 115     DS3_u, DS3_v, DS3_uuv, DS3_uvv,
 116     DC1_t, DC2_t;
 117   Surface->D3(u, v, S, DS1_u, DS1_v, DS2_u, DS2_v, DS2_uv,
 118     DS3_u, DS3_v, DS3_uuv, DS3_uvv);
 119   Curve->D2(t, C, DC1_t, DC2_t);
 120   gp_Vec Ort(C, S);
 121
 122   gp_Vec2d dE_dt(-DC1_t*DS1_u, -DC1_t*DS1_v);
 123   gp_XY dE_du(DS1_u*DS1_u + Ort*DS2_u,
 124     DS1_u*DS1_v + Ort*DS2_uv);
 125   gp_XY dE_dv(DS1_v*DS1_u + Ort*DS2_uv,
 126     DS1_v*DS1_v + Ort*DS2_v);
 127
 128   Standard_Real det = dE_du.X()*dE_dv.Y() - dE_du.Y()*dE_dv.X();
 129   if (fabs(det) < gp::Resolution()) Standard_ConstructionError::Raise();
 130
 131   gp_Mat2d M(gp_XY(dE_dv.Y()/det, -dE_du.Y()/det),
 132     gp_XY(-dE_dv.X()/det, dE_du.X()/det));
 133
 134   // First derivative
 135   V1 = - gp_Vec2d(gp_Vec2d(M.Row(1))*dE_dt, gp_Vec2d(M.Row(2))*dE_dt);
 136
 137   /* Second derivative */
 138
 139   // Computation of d2E_dt2 = S1
 140   gp_Vec2d d2E_dt(-DC2_t*DS1_u, -DC2_t*DS1_v);
 141
 142   // Computation of 2*(d2E/dtdX)(dX/dt) = S2
 143   gp_Vec2d d2E1_dtdX(-DC1_t*DS2_u,
 144     -DC1_t*DS2_uv);
 145   gp_Vec2d d2E2_dtdX(-DC1_t*DS2_uv,
 146     -DC1_t*DS2_v);
 147   gp_Vec2d S2 = 2*gp_Vec2d(d2E1_dtdX*V1, d2E2_dtdX*V1);
 148
 149   // Computation of (d2E/dX2)*(dX/dt)2 = S3
 150
 151   // Row11 = (d2E1/du2, d2E1/dudv)
 152   Standard_Real tmp;
 153   gp_Vec2d Row11(3*DS1_u*DS2_u + Ort*DS3_u,
 154     tmp = 2*DS1_u*DS2_uv +
 155     DS1_v*DS2_u + Ort*DS3_uuv);
 156
 157   // Row12 = (d2E1/dudv, d2E1/dv2)
 158   gp_Vec2d Row12(tmp, DS2_v*DS1_u + 2*DS1_v*DS2_uv +
 159     Ort*DS3_uvv);
 160
 161   // Row21 = (d2E2/du2, d2E2/dudv)
 162   gp_Vec2d Row21(DS2_u*DS1_v + 2*DS1_u*DS2_uv + Ort*DS3_uuv,
 163     tmp = 2*DS2_uv*DS1_v + DS1_u*DS2_v + Ort*DS3_uvv);
 164
 165   // Row22 = (d2E2/duv, d2E2/dvdv)
 166   gp_Vec2d Row22(tmp, 3*DS1_v*DS2_v + Ort*DS3_v);
 167
 168   gp_Vec2d S3(V1*gp_Vec2d(Row11*V1, Row12*V1),
 169     V1*gp_Vec2d(Row21*V1, Row22*V1));
 170
 171   gp_Vec2d Sum = d2E_dt + S2 + S3;
 172
 173   V2 = - gp_Vec2d(gp_Vec2d(M.Row(1))*Sum, gp_Vec2d(M.Row(2))*Sum);
 174 }
 175 //=======================================================================
 176 //function : d1CurveOnSurf
 177 //purpose  : computes first derivative of the 3d projected curve
 178 //=======================================================================
 179
 180 #if 0
 181 static void d1CurvOnSurf(const Standard_Real t,
 182   const Standard_Real u,
 183   const Standard_Real v,
 184   gp_Vec& V,
 185   const Handle(Adaptor3d_HCurve)& Curve,
 186   const Handle(Adaptor3d_HSurface)& Surface)
 187 {
 188   gp_Pnt S, C;
 189   gp_Vec2d V2d;
 190   gp_Vec DS1_u, DS1_v, DS2_u, DS2_uv, DS2_v, DC1_t;
 191   Surface->D2(u, v, S, DS1_u, DS1_v, DS2_u, DS2_v, DS2_uv);
 192   Curve->D1(t, C, DC1_t);
 193   gp_Vec Ort(C, S);// Ort = S - C
 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()) Standard_ConstructionError::Raise();
 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   V2d = - gp_Vec2d(gp_Vec2d(M.Row(1))*dE_dt, gp_Vec2d(M.Row(2))*dE_dt);
 208
 209   V = DS1_u * V2d.X() + DS1_v * V2d.Y();
 210
 211 }
 212 #endif
 213
 214 //=======================================================================
 215 //function : d2CurveOnSurf
 216 //purpose  : computes second derivative of the 3D projected curve
 217 //=======================================================================
 218
 219 static void d2CurvOnSurf(const Standard_Real t,
 220   const Standard_Real u,
 221   const Standard_Real v,
 222   gp_Vec& V1 , gp_Vec& V2 ,
 223   const Handle(Adaptor3d_HCurve)& Curve,
 224   const Handle(Adaptor3d_HSurface)& Surface)
 225 {
 226   gp_Pnt S, C;
 227   gp_Vec2d V12d,V22d;
 228   gp_Vec DS1_u, DS1_v, DS2_u, DS2_uv, DS2_v,
 229     DS3_u, DS3_v, DS3_uuv, DS3_uvv,
 230     DC1_t, DC2_t;
 231   Surface->D3(u, v, S, DS1_u, DS1_v, DS2_u, DS2_v, DS2_uv,
 232     DS3_u, DS3_v, DS3_uuv, DS3_uvv);
 233   Curve->D2(t, C, DC1_t, DC2_t);
 234   gp_Vec Ort(C, S);
 235
 236   gp_Vec2d dE_dt(-DC1_t*DS1_u, -DC1_t*DS1_v);
 237   gp_XY dE_du(DS1_u*DS1_u + Ort*DS2_u,
 238     DS1_u*DS1_v + Ort*DS2_uv);
 239   gp_XY dE_dv(DS1_v*DS1_u + Ort*DS2_uv,
 240     DS1_v*DS1_v + Ort*DS2_v);
 241
 242   Standard_Real det = dE_du.X()*dE_dv.Y() - dE_du.Y()*dE_dv.X();
 243   if (fabs(det) < gp::Resolution()) Standard_ConstructionError::Raise();
 244
 245   gp_Mat2d M(gp_XY(dE_dv.Y()/det, -dE_du.Y()/det),
 246     gp_XY(-dE_dv.X()/det, dE_du.X()/det));
 247
 248   // First derivative
 249   V12d = - gp_Vec2d(gp_Vec2d(M.Row(1))*dE_dt, gp_Vec2d(M.Row(2))*dE_dt);
 250
 251   /* Second derivative */
 252
 253   // Computation of d2E_dt2 = S1
 254   gp_Vec2d d2E_dt(-DC2_t*DS1_u, -DC2_t*DS1_v);
 255
 256   // Computation of 2*(d2E/dtdX)(dX/dt) = S2
 257   gp_Vec2d d2E1_dtdX(-DC1_t*DS2_u,
 258     -DC1_t*DS2_uv);
 259   gp_Vec2d d2E2_dtdX(-DC1_t*DS2_uv,
 260     -DC1_t*DS2_v);
 261   gp_Vec2d S2 = 2*gp_Vec2d(d2E1_dtdX*V12d, d2E2_dtdX*V12d);
 262
 263   // Computation of (d2E/dX2)*(dX/dt)2 = S3
 264
 265   // Row11 = (d2E1/du2, d2E1/dudv)
 266   Standard_Real tmp;
 267   gp_Vec2d Row11(3*DS1_u*DS2_u + Ort*DS3_u,
 268     tmp = 2*DS1_u*DS2_uv +
 269     DS1_v*DS2_u + Ort*DS3_uuv);
 270
 271   // Row12 = (d2E1/dudv, d2E1/dv2)
 272   gp_Vec2d Row12(tmp, DS2_v*DS1_u + 2*DS1_v*DS2_uv +
 273     Ort*DS3_uvv);
 274
 275   // Row21 = (d2E2/du2, d2E2/dudv)
 276   gp_Vec2d Row21(DS2_u*DS1_v + 2*DS1_u*DS2_uv + Ort*DS3_uuv,
 277     tmp = 2*DS2_uv*DS1_v + DS1_u*DS2_v + Ort*DS3_uvv);
 278
 279   // Row22 = (d2E2/duv, d2E2/dvdv)
 280   gp_Vec2d Row22(tmp, 3*DS1_v*DS2_v + Ort*DS3_v);
 281
 282   gp_Vec2d S3(V12d*gp_Vec2d(Row11*V12d, Row12*V12d),
 283     V12d*gp_Vec2d(Row21*V12d, Row22*V12d));
 284
 285   gp_Vec2d Sum = d2E_dt + S2 + S3;
 286
 287   V22d = - gp_Vec2d(gp_Vec2d(M.Row(1))*Sum, gp_Vec2d(M.Row(2))*Sum);
 288
 289   V1 = DS1_u * V12d.X() + DS1_v * V12d.Y();
 290   V2 =     DS2_u * V12d.X() *V12d.X()
 291     +     DS1_u * V22d.X()
 292     + 2 * DS2_uv * V12d.X() *V12d.Y()
 293     +     DS2_v * V12d.Y() * V12d.Y()
 294     +     DS1_v * V22d.Y();
 295 }
 296
 297 //=======================================================================
 298 //function : ExactBound
 299 //purpose  : computes exact boundary point
 300 //=======================================================================
 301
 302 static Standard_Boolean ExactBound(gp_Pnt& Sol,
 303   const Standard_Real NotSol,
 304   const Standard_Real Tol,
 305   const Standard_Real TolU,
 306   const Standard_Real TolV,
 307   const Handle(Adaptor3d_HCurve)& Curve,
 308   const Handle(Adaptor3d_HSurface)& Surface)
 309 {
 310   Standard_Real U0, V0, t, t1, t2, FirstU, LastU, FirstV, LastV;
 311   gp_Pnt2d POnS;
 312   U0 = Sol.Y();
 313   V0 = Sol.Z();
 314   FirstU = Surface->FirstUParameter();
 315   LastU = Surface->LastUParameter();
 316   FirstV = Surface->FirstVParameter();
 317   LastV = Surface->LastVParameter();
 318   // Here we have to compute the boundary that projection is going to intersect
 319   gp_Vec2d D2d;
 320   //these variables are to estimate which boundary has more apportunity
 321   //to be intersected
 322   Standard_Real RU1, RU2, RV1, RV2;
 323   d1(Sol.X(), U0, V0, D2d, Curve, Surface);
 324   // Here we assume that D2d != (0, 0)
 325   if(Abs(D2d.X()) < gp::Resolution())
 326   {
 327     RU1 = Precision::Infinite();
 328     RU2 = Precision::Infinite();
 329     RV1 = V0 - FirstV;
 330     RV2 = LastV - V0;
 331   }
 332   else if(Abs(D2d.Y()) < gp::Resolution())
 333   {
 334     RU1 = U0 - FirstU;
 335     RU2 = LastU - U0;
 336     RV1 = Precision::Infinite();
 337     RV2 = Precision::Infinite();
 338   }
 339   else
 340   {
 341     RU1 = gp_Pnt2d(U0, V0).
 342       Distance(gp_Pnt2d(FirstU, V0 + (FirstU - U0)*D2d.Y()/D2d.X()));
 343     RU2 = gp_Pnt2d(U0, V0).
 344       Distance(gp_Pnt2d(LastU, V0 + (LastU - U0)*D2d.Y()/D2d.X()));
 345     RV1 = gp_Pnt2d(U0, V0).
 346       Distance(gp_Pnt2d(U0 + (FirstV - V0)*D2d.X()/D2d.Y(), FirstV));
 347     RV2 = gp_Pnt2d(U0, V0).
 348       Distance(gp_Pnt2d(U0 + (LastV - V0)*D2d.X()/D2d.Y(), LastV));
 349   }
 350   TColgp_SequenceOfPnt Seq;
 351   Seq.Append(gp_Pnt(FirstU, RU1, 2));
 352   Seq.Append(gp_Pnt(LastU, RU2, 2));
 353   Seq.Append(gp_Pnt(FirstV, RV1, 3));
 354   Seq.Append(gp_Pnt(LastV, RV2, 3));
 355   Standard_Integer i, j;
 356   for(i = 1; i <= 3; i++)
 357     for(j = 1; j <= 4-i; j++)
 358       if(Seq(j).Y() < Seq(j+1).Y())
 359       {
 360         gp_Pnt swp;
 361         swp = Seq.Value(j+1);
 362         Seq.ChangeValue(j+1) = Seq.Value(j);
 363         Seq.ChangeValue(j) = swp;
 364       }
 365
 366       t = Sol.X();
 367       t1 = Min(Sol.X(), NotSol);
 368       t2 = Max(Sol.X(), NotSol);
 369
 370       Standard_Boolean isDone = Standard_False;
 371       while (!Seq.IsEmpty())
 372       {
 373         gp_Pnt P;
 374         P = Seq.Last();
 375         Seq.Remove(Seq.Length());
 376         ProjLib_PrjResolve aPrjPS(Curve->Curve(),
 377           Surface->Surface(),
 378           Standard_Integer(P.Z()));
 379         if(Standard_Integer(P.Z()) == 2)
 380         {
 381           aPrjPS.Perform(t, P.X(), V0, gp_Pnt2d(Tol, TolV),
 382             gp_Pnt2d(t1, Surface->FirstVParameter()),
 383             gp_Pnt2d(t2, Surface->LastVParameter()), FuncTol);
 384           if(!aPrjPS.IsDone()) continue;
 385           POnS = aPrjPS.Solution();
 386           Sol = gp_Pnt(POnS.X(), P.X(), POnS.Y());
 387           isDone = Standard_True;
 388           break;
 389         }
 390         else
 391         {
 392           aPrjPS.Perform(t, U0, P.X(), gp_Pnt2d(Tol, TolU),
 393             gp_Pnt2d(t1, Surface->FirstUParameter()),
 394             gp_Pnt2d(t2, Surface->LastUParameter()), FuncTol);
 395           if(!aPrjPS.IsDone()) continue;
 396           POnS = aPrjPS.Solution();
 397           Sol = gp_Pnt(POnS.X(), POnS.Y(), P.X());
 398           isDone = Standard_True;
 399           break;
 400         }
 401       }
 402
 403       return isDone;
 404 }
 405
 406 //=======================================================================
 407 //function : DichExactBound
 408 //purpose  : computes exact boundary point
 409 //=======================================================================
 410
 411 static void DichExactBound(gp_Pnt& Sol,
 412   const Standard_Real NotSol,
 413   const Standard_Real Tol,
 414   const Standard_Real TolU,
 415   const Standard_Real TolV,
 416   const Handle(Adaptor3d_HCurve)& Curve,
 417   const Handle(Adaptor3d_HSurface)& Surface)
 418 {
 419 #ifdef OCCT_DEBUG_CHRONO
 420   InitChron(chr_dicho_bound);
 421 #endif
 422
 423   Standard_Real U0, V0, t;
 424   gp_Pnt2d POnS;
 425   U0 = Sol.Y();
 426   V0 = Sol.Z();
 427   ProjLib_PrjResolve aPrjPS(Curve->Curve(), Surface->Surface(), 1);
 428
 429   Standard_Real aNotSol = NotSol;
 430   while (fabs(Sol.X() - aNotSol) > Tol)
 431   {
 432     t = (Sol.X() + aNotSol)/2;
 433     aPrjPS.Perform(t, U0, V0, gp_Pnt2d(TolU, TolV),
 434       gp_Pnt2d(Surface->FirstUParameter(),Surface->FirstVParameter()),
 435       gp_Pnt2d(Surface->LastUParameter(),Surface->LastVParameter()),
 436       FuncTol, Standard_True);
 437
 438     if (aPrjPS.IsDone())
 439     {
 440       POnS = aPrjPS.Solution();
 441       Sol = gp_Pnt(t, POnS.X(), POnS.Y());
 442       U0=Sol.Y();
 443       V0=Sol.Z();
 444     }
 445     else aNotSol = t;
 446   }
 447 #ifdef OCCT_DEBUG_CHRONO
 448   ResultChron(chr_dicho_bound,t_dicho_bound);
 449   dicho_bound_count++;
 450 #endif
 451 }
 452
 453 //=======================================================================
 454 //function : InitialPoint
 455 //purpose  :
 456 //=======================================================================
 457
 458 static Standard_Boolean InitialPoint(const gp_Pnt& Point,
 459   const Standard_Real t,
 460   const Handle(Adaptor3d_HCurve)& C,
 461   const Handle(Adaptor3d_HSurface)& S,
 462   const Standard_Real TolU,
 463   const Standard_Real TolV,
 464   Standard_Real& U,
 465   Standard_Real& V)
 466 {
 467
 468   ProjLib_PrjResolve aPrjPS(C->Curve(), S->Surface(), 1);
 469   Standard_Real ParU,ParV;
 470   Extrema_ExtPS aExtPS;
 471   aExtPS.Initialize(S->Surface(), S->FirstUParameter(),
 472     S->LastUParameter(), S->FirstVParameter(),
 473     S->LastVParameter(), TolU, TolV);
 474
 475   aExtPS.Perform(Point);
 476   Standard_Integer argmin = 0;
 477   if (aExtPS.IsDone() && aExtPS.NbExt())
 478   {
 479     Standard_Integer i, Nend;
 480     // Search for the nearest solution which is also a normal projection
 481     Nend = aExtPS.NbExt();
 482     for(i = 1; i <= Nend; i++)
 483     {
 484       Extrema_POnSurf POnS = aExtPS.Point(i);
 485       POnS.Parameter(ParU, ParV);
 486       aPrjPS.Perform(t, ParU, ParV, gp_Pnt2d(TolU, TolV),
 487         gp_Pnt2d(S->FirstUParameter(), S->FirstVParameter()),
 488         gp_Pnt2d(S->LastUParameter(), S->LastVParameter()),
 489         FuncTol, Standard_True);
 490       if(aPrjPS.IsDone() )
 491         if  (argmin == 0 || aExtPS.SquareDistance(i) < aExtPS.SquareDistance(argmin)) argmin = i;
 492     }
 493   }
 494   if( argmin == 0 ) return Standard_False;
 495   else
 496   {
 497     Extrema_POnSurf POnS = aExtPS.Point(argmin);
 498     POnS.Parameter(U, V);
 499     return Standard_True;
 500   }
 501 }
 502
 503 //=======================================================================
 504 //function : ProjLib_CompProjectedCurve
 505 //purpose  :
 506 //=======================================================================
 507
 508 ProjLib_CompProjectedCurve::ProjLib_CompProjectedCurve()
 509 : myNbCurves(0),
 510   myTolU    (0.0),
 511   myTolV    (0.0),
 512   myMaxDist (0.0)
 513 {
 514 }
 515
 516 //=======================================================================
 517 //function : ProjLib_CompProjectedCurve
 518 //purpose  :
 519 //=======================================================================
 520
 521 ProjLib_CompProjectedCurve::ProjLib_CompProjectedCurve
 522                            (const Handle(Adaptor3d_HSurface)& theSurface,
 523                             const Handle(Adaptor3d_HCurve)&   theCurve,
 524                             const Standard_Real               theTolU,
 525                             const Standard_Real               theTolV)
 526 : mySurface (theSurface),
 527   myCurve   (theCurve),
 528   myNbCurves(0),
 529   mySequence(new ProjLib_HSequenceOfHSequenceOfPnt()),
 530   myTolU    (theTolU),
 531   myTolV    (theTolV),
 532   myMaxDist (-1.0)
 533 {
 534   Init();
 535 }
 536
 537 //=======================================================================
 538 //function : ProjLib_CompProjectedCurve
 539 //purpose  :
 540 //=======================================================================
 541
 542 ProjLib_CompProjectedCurve::ProjLib_CompProjectedCurve
 543                            (const Handle(Adaptor3d_HSurface)& theSurface,
 544                             const Handle(Adaptor3d_HCurve)&   theCurve,
 545                             const Standard_Real               theTolU,
 546                             const Standard_Real               theTolV,
 547                             const Standard_Real               theMaxDist)
 548 : mySurface (theSurface),
 549   myCurve   (theCurve),
 550   myNbCurves(0),
 551   mySequence(new ProjLib_HSequenceOfHSequenceOfPnt()),
 552   myTolU    (theTolU),
 553   myTolV    (theTolV),
 554   myMaxDist (theMaxDist)
 555 {
 556   Init();
 557 }
 558
 559 //=======================================================================
 560 //function : Init
 561 //purpose  :
 562 //=======================================================================
 563
 564 void ProjLib_CompProjectedCurve::Init()
 565 {
 566   myTabInt.Nullify();
 567
 568   Standard_Real Tol;// Tolerance for ExactBound
 569   Standard_Integer i, Nend = 0;
 570   Standard_Boolean FromLastU=Standard_False;
 571
 572   //new part (to discard far solutions)
 573   Standard_Real TolC = Precision::Confusion(), TolS = Precision::Confusion();
 574   Extrema_ExtCS CExt(myCurve->Curve(),
 575     mySurface->Surface(),
 576     TolC,
 577     TolS);
 578   if (CExt.IsDone() && CExt.NbExt())
 579   {
 580     // Search for the minimum solution
 581     Nend = CExt.NbExt();
 582     if(myMaxDist > 0 &&
 583        // Avoid usage of extrema result that can be wrong for extrusion
 584        mySurface->GetType() != GeomAbs_SurfaceOfExtrusion)
 585     {
 586       Standard_Real min_val2;
 587       min_val2 = CExt.SquareDistance(1);
 588       for(i = 2; i <= Nend; i++)
 589         if (CExt.SquareDistance(i) < min_val2) min_val2 = CExt.SquareDistance(i);
 590       if (min_val2 > myMaxDist * myMaxDist)
 591         return;
 592     }
 593   }
 594   // end of new part
 595
 596   Standard_Real FirstU, LastU, Step, SearchStep, WalkStep, t;
 597
 598   FirstU = myCurve->FirstParameter();
 599   LastU  = myCurve->LastParameter();
 600   const Standard_Real GlobalMinStep = 1.e-4;
 601   //<GlobalMinStep> is sufficiently small to provide solving from initial point
 602   //and, on the other hand, it is sufficiently large to avoid too close solutions.
 603   const Standard_Real MinStep = 0.01*(LastU - FirstU),
 604     MaxStep = 0.1*(LastU - FirstU);
 605   SearchStep = 10*MinStep;
 606   Step = SearchStep;
 607
 608   //Initialization of aPrjPS
 609   Standard_Real Uinf = mySurface->FirstUParameter();
 610   Standard_Real Usup = mySurface->LastUParameter();
 611   Standard_Real Vinf = mySurface->FirstVParameter();
 612   Standard_Real Vsup = mySurface->LastVParameter();
 613
 614   ProjLib_PrjResolve aPrjPS(myCurve->Curve(), mySurface->Surface(), 1);
 615
 616   t = FirstU;
 617   Standard_Boolean new_part;
 618   Standard_Real prevDeb=0.;
 619   Standard_Boolean SameDeb=Standard_False;
 620
 621
 622   gp_Pnt Triple, prevTriple;
 623
 624   //Basic loop
 625   while(t <= LastU)
 626   {
 627     // Search for the beginning of a new continuous part
 628     // to avoid infinite computation in some difficult cases.
 629     new_part = Standard_False;
 630     if(t > FirstU && Abs(t-prevDeb) <= Precision::PConfusion()) SameDeb=Standard_True;
 631     while(t <= LastU && !new_part && !FromLastU && !SameDeb)
 632     {
 633       prevDeb=t;
 634       if (t == LastU) FromLastU=Standard_True;
 635       Standard_Boolean initpoint=Standard_False;
 636       Standard_Real U = 0., V = 0.;
 637       gp_Pnt CPoint;
 638       Standard_Real ParT,ParU,ParV;
 639
 640       // Search an initial point in the list of Extrema Curve-Surface
 641       if(Nend != 0 && !CExt.IsParallel())
 642       {
 643         for (i=1;i<=Nend;i++)
 644         {
 645           Extrema_POnCurv P1;
 646           Extrema_POnSurf P2;
 647           CExt.Points(i,P1,P2);
 648           ParT=P1.Parameter();
 649           P2.Parameter(ParU, ParV);
 650
 651           aPrjPS.Perform(ParT, ParU, ParV, gp_Pnt2d(myTolU, myTolV),
 652             gp_Pnt2d(mySurface->FirstUParameter(),mySurface->FirstVParameter()),
 653             gp_Pnt2d(mySurface->LastUParameter(), mySurface->LastVParameter()),
 654             FuncTol, Standard_True);
 655           if ( aPrjPS.IsDone() && P1.Parameter() > Max(FirstU,t-Step+Precision::PConfusion())
 656             && P1.Parameter() <= t)
 657           {
 658             t=ParT;
 659             U=ParU;
 660             V=ParV;
 661             CPoint=P1.Value();
 662             initpoint = Standard_True;
 663             break;
 664           }
 665         }
 666       }
 667       if (!initpoint)
 668       {
 669         myCurve->D0(t,CPoint);
 670 #ifdef OCCT_DEBUG_CHRONO
 671         InitChron(chr_init_point);
 672 #endif
 673         initpoint=InitialPoint(CPoint, t,myCurve,mySurface, myTolU, myTolV, U, V);
 674 #ifdef OCCT_DEBUG_CHRONO
 675         ResultChron(chr_init_point,t_init_point);
 676         init_point_count++;
 677 #endif
 678       }
 679       if(initpoint)
 680       {
 681         // When U or V lie on surface joint in some cases we cannot use them
 682         // as initial point for aPrjPS, so we switch them
 683         gp_Vec2d D;
 684
 685         if ((mySurface->IsUPeriodic() &&
 686             Abs(Usup - Uinf - mySurface->UPeriod()) < Precision::Confusion()) ||
 687             (mySurface->IsVPeriodic() &&
 688             Abs(Vsup - Vinf - mySurface->VPeriod()) < Precision::Confusion()))
 689         {
 690           if((Abs(U - Uinf) < mySurface->UResolution(Precision::PConfusion())) &&
 691             mySurface->IsUPeriodic())
 692           {
 693             d1(t, U, V, D, myCurve, mySurface);
 694             if (D.X() < 0 ) U = Usup;
 695           }
 696           else if((Abs(U - Usup) < mySurface->UResolution(Precision::PConfusion())) &&
 697             mySurface->IsUPeriodic())
 698           {
 699             d1(t, U, V, D, myCurve, mySurface);
 700             if (D.X() > 0) U = Uinf;
 701           }
 702
 703           if((Abs(V - Vinf) < mySurface->VResolution(Precision::PConfusion())) &&
 704             mySurface->IsVPeriodic())
 705           {
 706             d1(t, U, V, D, myCurve, mySurface);
 707             if (D.Y() < 0) V = Vsup;
 708           }
 709           else if((Abs(V - Vsup) <= mySurface->VResolution(Precision::PConfusion())) &&
 710             mySurface->IsVPeriodic())
 711           {
 712             d1(t, U, V, D, myCurve, mySurface);
 713             if (D.Y() > 0) V = Vinf;
 714           }
 715         }
 716
 717
 718         if (myMaxDist > 0)
 719         {
 720           // Here we are going to stop if the distance between projection and
 721           // corresponding curve point is greater than myMaxDist
 722           gp_Pnt POnS;
 723           Standard_Real d;
 724           mySurface->D0(U, V, POnS);
 725           d = CPoint.Distance(POnS);
 726           if (d > myMaxDist)
 727           {
 728             mySequence->Clear();
 729             myNbCurves = 0;
 730             return;
 731           }
 732         }
 733         Triple = gp_Pnt(t, U, V);
 734         if (t != FirstU)
 735         {
 736           //Search for exact boundary point
 737           Tol = Min(myTolU, myTolV);
 738           gp_Vec2d aD;
 739           d1(Triple.X(), Triple.Y(), Triple.Z(), aD, myCurve, mySurface);
 740           Tol /= Max(Abs(aD.X()), Abs(aD.Y()));
 741
 742           if(!ExactBound(Triple, t - Step, Tol,
 743             myTolU, myTolV, myCurve, mySurface))
 744           {
 745 #ifdef OCCT_DEBUG
 746             cout<<"There is a problem with ExactBound computation"<<endl;
 747 #endif
 748             DichExactBound(Triple, t - Step, Tol, myTolU, myTolV,
 749               myCurve, mySurface);
 750           }
 751         }
 752         new_part = Standard_True;
 753       }
 754       else
 755       {
 756         if(t == LastU) break;
 757         t += Step;
 758         if(t>LastU)
 759         {
 760           Step =Step+LastU-t;
 761           t=LastU;
 762         }
 763       }
 764     }
 765     if (!new_part) break;
 766
 767
 768     //We have found a new continuous part
 769     Handle(TColgp_HSequenceOfPnt) hSeq = new TColgp_HSequenceOfPnt();
 770     mySequence->Append(hSeq);
 771     myNbCurves++;
 772     mySequence->Value(myNbCurves)->Append(Triple);
 773     prevTriple = Triple;
 774
 775     if (Triple.X() == LastU) break;//return;
 776
 777     //Computation of WalkStep
 778     gp_Vec D1, D2;
 779     Standard_Real MagnD1, MagnD2;
 780     d2CurvOnSurf(Triple.X(), Triple.Y(), Triple.Z(), D1, D2, myCurve, mySurface);
 781     MagnD1 = D1.Magnitude();
 782     MagnD2 = D2.Magnitude();
 783     if(MagnD2 < Precision::Confusion()) WalkStep = MaxStep;
 784     else WalkStep = Min(MaxStep, Max(MinStep, 0.1*MagnD1/MagnD2));
 785
 786     Step = WalkStep;
 787
 788     t = Triple.X() + Step;
 789     if (t > LastU) t = LastU;
 790     Standard_Real prevStep = Step;
 791     Standard_Real U0, V0;
 792     gp_Pnt2d aLowBorder(mySurface->FirstUParameter(),mySurface->FirstVParameter());
 793     gp_Pnt2d aUppBorder(mySurface->LastUParameter(), mySurface->LastVParameter());
 794     gp_Pnt2d aTol(myTolU, myTolV);
 795     //Here we are trying to prolong continuous part
 796     while (t <= LastU && new_part)
 797     {
 798
 799       U0 = Triple.Y() + (Step / prevStep) * (Triple.Y() - prevTriple.Y());
 800       V0 = Triple.Z() + (Step / prevStep) * (Triple.Z() - prevTriple.Z());
 801       // adjust U0 to be in [mySurface->FirstUParameter(),mySurface->LastUParameter()]
 802       U0 = Min(Max(U0, aLowBorder.X()), aUppBorder.X());
 803       // adjust V0 to be in [mySurface->FirstVParameter(),mySurface->LastVParameter()]
 804       V0 = Min(Max(V0, aLowBorder.Y()), aUppBorder.Y());
 805
 806
 807       aPrjPS.Perform(t, U0, V0, aTol,
 808                      aLowBorder, aUppBorder, FuncTol, Standard_True);
 809       if(!aPrjPS.IsDone())
 810       {
 811         if (Step <= GlobalMinStep)
 812         {
 813           //Search for exact boundary point
 814           Tol = Min(myTolU, myTolV);
 815           gp_Vec2d D;
 816           d1(Triple.X(), Triple.Y(), Triple.Z(), D, myCurve, mySurface);
 817           Tol /= Max(Abs(D.X()), Abs(D.Y()));
 818
 819           if(!ExactBound(Triple, t, Tol, myTolU, myTolV,
 820             myCurve, mySurface))
 821           {
 822 #ifdef OCCT_DEBUG
 823             cout<<"There is a problem with ExactBound computation"<<endl;
 824 #endif
 825             DichExactBound(Triple, t, Tol, myTolU, myTolV,
 826               myCurve, mySurface);
 827           }
 828
 829           if((Triple.X() - mySequence->Value(myNbCurves)->Value(mySequence->Value(myNbCurves)->Length()).X()) > 1.e-10)
 830             mySequence->Value(myNbCurves)->Append(Triple);
 831           if((LastU - Triple.X()) < Tol) {t = LastU + 1; break;}//return;
 832
 833           Step = SearchStep;
 834           t = Triple.X() + Step;
 835           if (t > (LastU-MinStep/2) )
 836           {
 837             Step =Step+LastU-t;
 838             t = LastU;
 839           }
 840           new_part = Standard_False;
 841         }
 842         else
 843         {
 844           // decrease step
 845           Standard_Real SaveStep = Step;
 846           Step /= 2.;
 847           t = Triple .X() + Step;
 848           if (t > (LastU-MinStep/4) )
 849           {
 850             Step =Step+LastU-t;
 851             if (Abs(Step - SaveStep) <= Precision::PConfusion())
 852               Step = GlobalMinStep; //to avoid looping
 853             t = LastU;
 854           }
 855         }
 856       }
 857       // Go further
 858       else
 859       {
 860         prevTriple = Triple;
 861         prevStep = Step;
 862         Triple = gp_Pnt(t, aPrjPS.Solution().X(), aPrjPS.Solution().Y());
 863
 864         // Check for possible local traps.
 865         UpdateTripleByTrapCriteria(Triple);
 866
 867         if((Triple.X() - mySequence->Value(myNbCurves)->Value(mySequence->Value(myNbCurves)->Length()).X()) > 1.e-10)
 868           mySequence->Value(myNbCurves)->Append(Triple);
 869         if (t == LastU) {t = LastU + 1; break;}//return;
 870         //Computation of WalkStep
 871         d2CurvOnSurf(Triple.X(), Triple.Y(), Triple.Z(), D1, D2, myCurve, mySurface);
 872         MagnD1 = D1.Magnitude();
 873         MagnD2 = D2.Magnitude();
 874         if(MagnD2 < Precision::Confusion() ) WalkStep = MaxStep;
 875         else WalkStep = Min(MaxStep, Max(MinStep, 0.1*MagnD1/MagnD2));
 876
 877         Step = WalkStep;
 878         t += Step;
 879         if (t > (LastU-MinStep/2) )
 880         {
 881           Step =Step+LastU-t;
 882           t = LastU;
 883         }
 884       }
 885     }
 886   }
 887   // Sequence post-proceeding.
 888   Standard_Integer j;
 889
 890   // 1. Removing poor parts
 891   Standard_Integer NbPart=myNbCurves;
 892   Standard_Integer ipart=1;
 893   for(i = 1; i <= NbPart; i++) {
 894     //    Standard_Integer NbPoints = mySequence->Value(i)->Length();
 895     if(mySequence->Value(ipart)->Length() < 2) {
 896       mySequence->Remove(ipart);
 897       myNbCurves--;
 898     }
 899     else ipart++;
 900   }
 901
 902   if(myNbCurves == 0) return;
 903
 904   // 2. Removing common parts of bounds
 905   for(i = 1; i < myNbCurves; i++)
 906   {
 907     if(mySequence->Value(i)->Value(mySequence->Value(i)->Length()).X() >=
 908       mySequence->Value(i+1)->Value(1).X())
 909       mySequence->ChangeValue(i+1)->ChangeValue(1).SetX(mySequence->Value(i)->Value(mySequence->Value(i)->Length()).X() + 1.e-12);
 910   }
 911
 912   // 3. Computation of the maximum distance from each part of curve to surface
 913
 914   myMaxDistance = new TColStd_HArray1OfReal(1, myNbCurves);
 915   myMaxDistance->Init(0);
 916   for(i = 1; i <= myNbCurves; i++)
 917     for(j = 1; j <= mySequence->Value(i)->Length(); j++)
 918     {
 919       gp_Pnt POnC, POnS, aTriple;
 920       Standard_Real Distance;
 921       aTriple = mySequence->Value(i)->Value(j);
 922       myCurve->D0(aTriple.X(), POnC);
 923       mySurface->D0(aTriple.Y(), aTriple.Z(), POnS);
 924       Distance = POnC.Distance(POnS);
 925       if (myMaxDistance->Value(i) < Distance)
 926         myMaxDistance->ChangeValue(i) = Distance;
 927     }
 928
 929
 930     // 4. Check the projection to be a single point
 931
 932     gp_Pnt2d Pmoy, Pcurr, P;
 933     Standard_Real AveU, AveV;
 934     mySnglPnts = new TColStd_HArray1OfBoolean(1, myNbCurves);
 935     for(i = 1; i <= myNbCurves; i++) mySnglPnts->SetValue(i, Standard_True);
 936
 937     for(i = 1; i <= myNbCurves; i++)
 938     {
 939       //compute an average U and V
 940
 941       for(j = 1, AveU = 0., AveV = 0.; j <= mySequence->Value(i)->Length(); j++)
 942       {
 943         AveU += mySequence->Value(i)->Value(j).Y();
 944         AveV += mySequence->Value(i)->Value(j).Z();
 945       }
 946       AveU /= mySequence->Value(i)->Length();
 947       AveV /= mySequence->Value(i)->Length();
 948
 949       Pmoy.SetCoord(AveU,AveV);
 950       for(j = 1; j <= mySequence->Value(i)->Length(); j++)
 951       {
 952         Pcurr =
 953           gp_Pnt2d(mySequence->Value(i)->Value(j).Y(), mySequence->Value(i)->Value(j).Z());
 954         if (Pcurr.Distance(Pmoy) > ((myTolU < myTolV) ? myTolV : myTolU))
 955         {
 956           mySnglPnts->SetValue(i, Standard_False);
 957           break;
 958         }
 959       }
 960     }
 961
 962     // 5. Check the projection to be an isoparametric curve of the surface
 963
 964     myUIso = new TColStd_HArray1OfBoolean(1, myNbCurves);
 965     for(i = 1; i <= myNbCurves; i++) myUIso->SetValue(i, Standard_True);
 966
 967     myVIso = new TColStd_HArray1OfBoolean(1, myNbCurves);
 968     for(i = 1; i <= myNbCurves; i++) myVIso->SetValue(i, Standard_True);
 969
 970     for(i = 1; i <= myNbCurves; i++) {
 971       if (IsSinglePnt(i, P)|| mySequence->Value(i)->Length() <=2) {
 972         myUIso->SetValue(i, Standard_False);
 973         myVIso->SetValue(i, Standard_False);
 974         continue;
 975       }
 976
 977       // new test for isoparametrics
 978
 979       if ( mySequence->Value(i)->Length() > 2) {
 980         //compute an average U and V
 981
 982         for(j = 1, AveU = 0., AveV = 0.; j <= mySequence->Value(i)->Length(); j++) {
 983           AveU += mySequence->Value(i)->Value(j).Y();
 984           AveV += mySequence->Value(i)->Value(j).Z();
 985         }
 986         AveU /= mySequence->Value(i)->Length();
 987         AveV /= mySequence->Value(i)->Length();
 988
 989         // is i-part U-isoparametric ?
 990         for(j = 1; j <= mySequence->Value(i)->Length(); j++)
 991         {
 992           if(Abs(mySequence->Value(i)->Value(j).Y() - AveU) > myTolU)
 993           {
 994             myUIso->SetValue(i, Standard_False);
 995             break;
 996           }
 997         }
 998
 999         // is i-part V-isoparametric ?
1000         for(j = 1; j <= mySequence->Value(i)->Length(); j++)
1001         {
1002           if(Abs(mySequence->Value(i)->Value(j).Z() - AveV) > myTolV)
1003           {
1004             myVIso->SetValue(i, Standard_False);
1005             break;
1006           }
1007         }
1008         //
1009       }
1010     }
1011 }
1012 //=======================================================================
1013 //function : Load
1014 //purpose  :
1015 //=======================================================================
1016
1017 void ProjLib_CompProjectedCurve::Load(const Handle(Adaptor3d_HSurface)& S)
1018 {
1019   mySurface = S;
1020 }
1021
1022 //=======================================================================
1023 //function : Load
1024 //purpose  :
1025 //=======================================================================
1026
1027 void ProjLib_CompProjectedCurve::Load(const Handle(Adaptor3d_HCurve)& C)
1028 {
1029   myCurve = C;
1030 }
1031
1032 //=======================================================================
1033 //function : GetSurface
1034 //purpose  :
1035 //=======================================================================
1036
1037 const Handle(Adaptor3d_HSurface)& ProjLib_CompProjectedCurve::GetSurface() const
1038 {
1039   return mySurface;
1040 }
1041
1042
1043 //=======================================================================
1044 //function : GetCurve
1045 //purpose  :
1046 //=======================================================================
1047
1048 const Handle(Adaptor3d_HCurve)& ProjLib_CompProjectedCurve::GetCurve() const
1049 {
1050   return myCurve;
1051 }
1052
1053 //=======================================================================
1054 //function : GetTolerance
1055 //purpose  :
1056 //=======================================================================
1057
1058 void ProjLib_CompProjectedCurve::GetTolerance(Standard_Real& TolU,
1059   Standard_Real& TolV) const
1060 {
1061   TolU = myTolU;
1062   TolV = myTolV;
1063 }
1064
1065 //=======================================================================
1066 //function : NbCurves
1067 //purpose  :
1068 //=======================================================================
1069
1070 Standard_Integer ProjLib_CompProjectedCurve::NbCurves() const
1071 {
1072   return myNbCurves;
1073 }
1074 //=======================================================================
1075 //function : Bounds
1076 //purpose  :
1077 //=======================================================================
1078
1079 void ProjLib_CompProjectedCurve::Bounds(const Standard_Integer Index,
1080   Standard_Real& Udeb,
1081   Standard_Real& Ufin) const
1082 {
1083   if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1084   Udeb = mySequence->Value(Index)->Value(1).X();
1085   Ufin = mySequence->Value(Index)->Value(mySequence->Value(Index)->Length()).X();
1086 }
1087 //=======================================================================
1088 //function : IsSinglePnt
1089 //purpose  :
1090 //=======================================================================
1091
1092 Standard_Boolean ProjLib_CompProjectedCurve::IsSinglePnt(const Standard_Integer Index, gp_Pnt2d& P) const
1093 {
1094   if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1095   P = gp_Pnt2d(mySequence->Value(Index)->Value(1).Y(), mySequence->Value(Index)->Value(1).Z());
1096   return mySnglPnts->Value(Index);
1097 }
1098
1099 //=======================================================================
1100 //function : IsUIso
1101 //purpose  :
1102 //=======================================================================
1103
1104 Standard_Boolean ProjLib_CompProjectedCurve::IsUIso(const Standard_Integer Index, Standard_Real& U) const
1105 {
1106   if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1107   U = mySequence->Value(Index)->Value(1).Y();
1108   return myUIso->Value(Index);
1109 }
1110 //=======================================================================
1111 //function : IsVIso
1112 //purpose  :
1113 //=======================================================================
1114
1115 Standard_Boolean ProjLib_CompProjectedCurve::IsVIso(const Standard_Integer Index, Standard_Real& V) const
1116 {
1117   if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1118   V = mySequence->Value(Index)->Value(1).Z();
1119   return myVIso->Value(Index);
1120 }
1121 //=======================================================================
1122 //function : Value
1123 //purpose  :
1124 //=======================================================================
1125
1126 gp_Pnt2d ProjLib_CompProjectedCurve::Value(const Standard_Real t) const
1127 {
1128   gp_Pnt2d P;
1129   D0(t, P);
1130   return P;
1131 }
1132 //=======================================================================
1133 //function : D0
1134 //purpose  :
1135 //=======================================================================
1136
1137 void ProjLib_CompProjectedCurve::D0(const Standard_Real U,gp_Pnt2d& P) const
1138 {
1139   Standard_Integer i, j;
1140   Standard_Real Udeb, Ufin;
1141   Standard_Boolean found = Standard_False;
1142
1143   for(i = 1; i <= myNbCurves; i++)
1144   {
1145     Bounds(i, Udeb, Ufin);
1146     if (U >= Udeb && U <= Ufin)
1147     {
1148       found = Standard_True;
1149       break;
1150     }
1151   }
1152   if (!found) Standard_DomainError::Raise("ProjLib_CompProjectedCurve::D0");
1153
1154   Standard_Real U0, V0;
1155
1156   Standard_Integer End = mySequence->Value(i)->Length();
1157   for(j = 1; j < End; j++)
1158     if ((U >= mySequence->Value(i)->Value(j).X()) && (U <= mySequence->Value(i)->Value(j + 1).X())) break;
1159
1160   //  U0 = mySequence->Value(i)->Value(j).Y();
1161   //  V0 = mySequence->Value(i)->Value(j).Z();
1162
1163   //  Cubic Interpolation
1164   if(mySequence->Value(i)->Length() < 4 ||
1165     (Abs(U-mySequence->Value(i)->Value(j).X()) <= Precision::PConfusion()) )
1166   {
1167     U0 = mySequence->Value(i)->Value(j).Y();
1168     V0 = mySequence->Value(i)->Value(j).Z();
1169   }
1170   else if (Abs(U-mySequence->Value(i)->Value(j+1).X())
1171     <= Precision::PConfusion())
1172   {
1173     U0 = mySequence->Value(i)->Value(j+1).Y();
1174     V0 = mySequence->Value(i)->Value(j+1).Z();
1175   }
1176   else
1177   {
1178     if (j == 1) j = 2;
1179     if (j > mySequence->Value(i)->Length() - 2)
1180       j = mySequence->Value(i)->Length() - 2;
1181
1182     gp_Vec2d I1, I2, I3, I21, I22, I31, Y1, Y2, Y3, Y4, Res;
1183     Standard_Real X1, X2, X3, X4;
1184
1185     X1 = mySequence->Value(i)->Value(j - 1).X();
1186     X2 = mySequence->Value(i)->Value(j).X();
1187     X3 = mySequence->Value(i)->Value(j + 1).X();
1188     X4 = mySequence->Value(i)->Value(j + 2).X();
1189
1190     Y1 = gp_Vec2d(mySequence->Value(i)->Value(j - 1).Y(),
1191       mySequence->Value(i)->Value(j - 1).Z());
1192     Y2 = gp_Vec2d(mySequence->Value(i)->Value(j).Y(),
1193       mySequence->Value(i)->Value(j).Z());
1194     Y3 = gp_Vec2d(mySequence->Value(i)->Value(j + 1).Y(),
1195       mySequence->Value(i)->Value(j + 1).Z());
1196     Y4 = gp_Vec2d(mySequence->Value(i)->Value(j + 2).Y(),
1197       mySequence->Value(i)->Value(j + 2).Z());
1198
1199     I1 = (Y1 - Y2)/(X1 - X2);
1200     I2 = (Y2 - Y3)/(X2 - X3);
1201     I3 = (Y3 - Y4)/(X3 - X4);
1202
1203     I21 = (I1 - I2)/(X1 - X3);
1204     I22 = (I2 - I3)/(X2 - X4);
1205
1206     I31 = (I21 - I22)/(X1 - X4);
1207
1208     Res = Y1 + (U - X1)*(I1 + (U - X2)*(I21 + (U - X3)*I31));
1209
1210     U0 = Res.X();
1211     V0 = Res.Y();
1212
1213     if(U0 < mySurface->FirstUParameter()) U0 = mySurface->FirstUParameter();
1214     else if(U0 > mySurface->LastUParameter()) U0 = mySurface->LastUParameter();
1215
1216     if(V0 < mySurface->FirstVParameter()) V0 = mySurface->FirstVParameter();
1217     else if(V0 > mySurface->LastVParameter()) V0 = mySurface->LastVParameter();
1218   }
1219   //End of cubic interpolation
1220
1221   ProjLib_PrjResolve aPrjPS(myCurve->Curve(), mySurface->Surface(), 1);
1222   aPrjPS.Perform(U, U0, V0, gp_Pnt2d(myTolU, myTolV),
1223     gp_Pnt2d(mySurface->FirstUParameter(), mySurface->FirstVParameter()),
1224     gp_Pnt2d(mySurface->LastUParameter(), mySurface->LastVParameter()));
1225   if (aPrjPS.IsDone())
1226     P = aPrjPS.Solution();
1227   else
1228   {
1229     gp_Pnt thePoint = myCurve->Value(U);
1230     Extrema_ExtPS aExtPS(thePoint, mySurface->Surface(), myTolU, myTolV);
1231     if (aExtPS.IsDone() && aExtPS.NbExt())
1232     {
1233       Standard_Integer k, Nend, imin = 1;
1234       // Search for the nearest solution which is also a normal projection
1235       Nend = aExtPS.NbExt();
1236       for(k = 2; k <= Nend; k++)
1237         if (aExtPS.SquareDistance(k) < aExtPS.SquareDistance(imin))
1238           imin = k;
1239       const Extrema_POnSurf& POnS = aExtPS.Point(imin);
1240       Standard_Real ParU,ParV;
1241       POnS.Parameter(ParU, ParV);
1242       P.SetCoord(ParU, ParV);
1243     }
1244     else
1245       P.SetCoord(U0,V0);
1246   }
1247 }
1248 //=======================================================================
1249 //function : D1
1250 //purpose  :
1251 //=======================================================================
1252
1253 void ProjLib_CompProjectedCurve::D1(const Standard_Real t,
1254   gp_Pnt2d& P,
1255   gp_Vec2d& V) const
1256 {
1257   Standard_Real u, v;
1258   D0(t, P);
1259   u = P.X();
1260   v = P.Y();
1261   d1(t, u, v, V, myCurve, mySurface);
1262 }
1263 //=======================================================================
1264 //function : D2
1265 //purpose  :
1266 //=======================================================================
1267
1268 void ProjLib_CompProjectedCurve::D2(const Standard_Real t,
1269   gp_Pnt2d& P,
1270   gp_Vec2d& V1,
1271   gp_Vec2d& V2) const
1272 {
1273   Standard_Real u, v;
1274   D0(t, P);
1275   u = P.X();
1276   v = P.Y();
1277   d2(t, u, v, V1, V2, myCurve, mySurface);
1278 }
1279 //=======================================================================
1280 //function : DN
1281 //purpose  :
1282 //=======================================================================
1283
1284 gp_Vec2d ProjLib_CompProjectedCurve::DN(const Standard_Real t,
1285   const Standard_Integer N) const
1286 {
1287   if (N < 1 ) Standard_OutOfRange::Raise("ProjLib_CompProjectedCurve : N must be greater than 0");
1288   else if (N ==1)
1289   {
1290     gp_Pnt2d P;
1291     gp_Vec2d V;
1292     D1(t,P,V);
1293     return V;
1294   }
1295   else if ( N==2)
1296   {
1297     gp_Pnt2d P;
1298     gp_Vec2d V1,V2;
1299     D2(t,P,V1,V2);
1300     return V2;
1301   }
1302   else if (N > 2 )
1303     Standard_NotImplemented::Raise("ProjLib_CompProjectedCurve::DN");
1304   return gp_Vec2d();
1305 }
1306
1307 //=======================================================================
1308 //function : GetSequence
1309 //purpose  :
1310 //=======================================================================
1311
1312 const Handle(ProjLib_HSequenceOfHSequenceOfPnt)& ProjLib_CompProjectedCurve::GetSequence() const
1313 {
1314   return mySequence;
1315 }
1316 //=======================================================================
1317 //function : FirstParameter
1318 //purpose  :
1319 //=======================================================================
1320
1321 Standard_Real ProjLib_CompProjectedCurve::FirstParameter() const
1322 {
1323   return myCurve->FirstParameter();
1324 }
1325
1326 //=======================================================================
1327 //function : LastParameter
1328 //purpose  :
1329 //=======================================================================
1330
1331 Standard_Real ProjLib_CompProjectedCurve::LastParameter() const
1332 {
1333   return myCurve->LastParameter();
1334 }
1335
1336 //=======================================================================
1337 //function : MaxDistance
1338 //purpose  :
1339 //=======================================================================
1340
1341 Standard_Real ProjLib_CompProjectedCurve::MaxDistance(const Standard_Integer Index) const
1342 {
1343   if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1344   return myMaxDistance->Value(Index);
1345 }
1346
1347 //=======================================================================
1348 //function : NbIntervals
1349 //purpose  :
1350 //=======================================================================
1351
1352 Standard_Integer ProjLib_CompProjectedCurve::NbIntervals(const GeomAbs_Shape S) const
1353 {
1354   const_cast<ProjLib_CompProjectedCurve*>(this)->myTabInt.Nullify();
1355   BuildIntervals(S);
1356   return myTabInt->Length() - 1;
1357 }
1358
1359 //=======================================================================
1360 //function : Intervals
1361 //purpose  :
1362 //=======================================================================
1363
1364 void ProjLib_CompProjectedCurve::Intervals(TColStd_Array1OfReal& T,const GeomAbs_Shape S) const
1365 {
1366   if (myTabInt.IsNull()) BuildIntervals (S);
1367   T = myTabInt->Array1();
1368 }
1369
1370 //=======================================================================
1371 //function : BuildIntervals
1372 //purpose  :
1373 //=======================================================================
1374
1375 void ProjLib_CompProjectedCurve::BuildIntervals(const GeomAbs_Shape S) const
1376 {
1377   GeomAbs_Shape SforS = GeomAbs_CN;
1378   switch(S) {
1379   case GeomAbs_C0:
1380     SforS = GeomAbs_C1;
1381     break;
1382   case GeomAbs_C1:
1383     SforS = GeomAbs_C2;
1384     break;
1385   case GeomAbs_C2:
1386     SforS = GeomAbs_C3;
1387     break;
1388   case GeomAbs_C3:
1389     SforS = GeomAbs_CN;
1390     break;
1391   case GeomAbs_CN:
1392     SforS = GeomAbs_CN;
1393     break;
1394   default:
1395     Standard_OutOfRange::Raise();
1396   }
1397   Standard_Integer i, j, k;
1398   Standard_Integer NbIntCur = myCurve->NbIntervals(S);
1399   Standard_Integer NbIntSurU = mySurface->NbUIntervals(SforS);
1400   Standard_Integer NbIntSurV = mySurface->NbVIntervals(SforS);
1401
1402   TColStd_Array1OfReal CutPntsT(1, NbIntCur+1);
1403   TColStd_Array1OfReal CutPntsU(1, NbIntSurU+1);
1404   TColStd_Array1OfReal CutPntsV(1, NbIntSurV+1);
1405
1406   myCurve->Intervals(CutPntsT, S);
1407   mySurface->UIntervals(CutPntsU, SforS);
1408   mySurface->VIntervals(CutPntsV, SforS);
1409
1410   Standard_Real Tl, Tr, Ul, Ur, Vl, Vr, Tol;
1411
1412   Handle(TColStd_HArray1OfReal) BArr = NULL,
1413     CArr = NULL,
1414     UArr = NULL,
1415     VArr = NULL;
1416
1417   // proccessing projection bounds
1418   BArr = new TColStd_HArray1OfReal(1, 2*myNbCurves);
1419   for(i = 1; i <= myNbCurves; i++)
1420     Bounds(i, BArr->ChangeValue(2*i - 1), BArr->ChangeValue(2*i));
1421
1422   // proccessing curve discontinuities
1423   if(NbIntCur > 1) {
1424     CArr = new TColStd_HArray1OfReal(1, NbIntCur - 1);
1425     for(i = 1; i <= CArr->Length(); i++)
1426       CArr->ChangeValue(i) = CutPntsT(i + 1);
1427   }
1428
1429   // proccessing U-surface discontinuities
1430   TColStd_SequenceOfReal TUdisc;
1431
1432   for(k = 2; k <= NbIntSurU; k++) {
1433     //    cout<<"CutPntsU("<<k<<") = "<<CutPntsU(k)<<endl;
1434     for(i = 1; i <= myNbCurves; i++)
1435       for(j = 1; j < mySequence->Value(i)->Length(); j++) {
1436         Ul = mySequence->Value(i)->Value(j).Y();
1437         Ur = mySequence->Value(i)->Value(j + 1).Y();
1438
1439         if(Abs(Ul - CutPntsU(k)) <= myTolU)
1440           TUdisc.Append(mySequence->Value(i)->Value(j).X());
1441         else if(Abs(Ur - CutPntsU(k)) <= myTolU)
1442           TUdisc.Append(mySequence->Value(i)->Value(j + 1).X());
1443         else if((Ul < CutPntsU(k) && CutPntsU(k) < Ur) ||
1444           (Ur < CutPntsU(k) && CutPntsU(k) < Ul))
1445         {
1446           Standard_Real V;
1447           V = (mySequence->Value(i)->Value(j).Z()
1448             + mySequence->Value(i)->Value(j +1).Z())/2;
1449           ProjLib_PrjResolve Solver(myCurve->Curve(), mySurface->Surface(), 2);
1450
1451           gp_Vec2d D;
1452           gp_Pnt Triple;
1453           Triple = mySequence->Value(i)->Value(j);
1454           d1(Triple.X(), Triple.Y(), Triple.Z(), D, myCurve, mySurface);
1455           if (Abs(D.X()) < Precision::Confusion())
1456             Tol = myTolU;
1457           else
1458             Tol = Min(myTolU, myTolU / Abs(D.X()));
1459
1460           Tl = mySequence->Value(i)->Value(j).X();
1461           Tr = mySequence->Value(i)->Value(j + 1).X();
1462
1463           Solver.Perform((Tl + Tr)/2, CutPntsU(k), V,
1464             gp_Pnt2d(Tol, myTolV),
1465             gp_Pnt2d(Tl, mySurface->FirstVParameter()),
1466             gp_Pnt2d(Tr, mySurface->LastVParameter()));
1467           //
1468           if(Solver.IsDone())
1469           {
1470             TUdisc.Append(Solver.Solution().X());
1471           }
1472         }
1473       }
1474   }
1475   for(i = 2; i <= TUdisc.Length(); i++)
1476     if(TUdisc(i) - TUdisc(i-1) < Precision::PConfusion())
1477       TUdisc.Remove(i--);
1478
1479   if(TUdisc.Length())
1480   {
1481     UArr = new TColStd_HArray1OfReal(1, TUdisc.Length());
1482     for(i = 1; i <= UArr->Length(); i++)
1483       UArr->ChangeValue(i) = TUdisc(i);
1484   }
1485   // proccessing V-surface discontinuities
1486   TColStd_SequenceOfReal TVdisc;
1487
1488   for(k = 2; k <= NbIntSurV; k++)
1489     for(i = 1; i <= myNbCurves; i++)
1490     {
1491       //      cout<<"CutPntsV("<<k<<") = "<<CutPntsV(k)<<endl;
1492       for(j = 1; j < mySequence->Value(i)->Length(); j++) {
1493
1494         Vl = mySequence->Value(i)->Value(j).Z();
1495         Vr = mySequence->Value(i)->Value(j + 1).Z();
1496
1497         if(Abs(Vl - CutPntsV(k)) <= myTolV)
1498           TVdisc.Append(mySequence->Value(i)->Value(j).X());
1499         else if (Abs(Vr - CutPntsV(k)) <= myTolV)
1500           TVdisc.Append(mySequence->Value(i)->Value(j + 1).X());
1501         else if((Vl < CutPntsV(k) && CutPntsV(k) < Vr) ||
1502           (Vr < CutPntsV(k) && CutPntsV(k) < Vl))
1503         {
1504           Standard_Real U;
1505           U = (mySequence->Value(i)->Value(j).Y()
1506             + mySequence->Value(i)->Value(j +1).Y())/2;
1507           ProjLib_PrjResolve Solver(myCurve->Curve(), mySurface->Surface(), 3);
1508
1509           gp_Vec2d D;
1510           gp_Pnt Triple;
1511           Triple = mySequence->Value(i)->Value(j);
1512           d1(Triple.X(), Triple.Y(), Triple.Z(), D, myCurve, mySurface);
1513           if (Abs(D.Y()) < Precision::Confusion())
1514             Tol = myTolV;
1515           else
1516             Tol = Min(myTolV, myTolV / Abs(D.Y()));
1517
1518           Tl = mySequence->Value(i)->Value(j).X();
1519           Tr = mySequence->Value(i)->Value(j + 1).X();
1520
1521           Solver.Perform((Tl + Tr)/2, U, CutPntsV(k),
1522             gp_Pnt2d(Tol, myTolV),
1523             gp_Pnt2d(Tl, mySurface->FirstUParameter()),
1524             gp_Pnt2d(Tr, mySurface->LastUParameter()));
1525           //
1526           if(Solver.IsDone())
1527           {
1528             TVdisc.Append(Solver.Solution().X());
1529           }
1530         }
1531       }
1532     }
1533     for(i = 2; i <= TVdisc.Length(); i++)
1534       if(TVdisc(i) - TVdisc(i-1) < Precision::PConfusion())
1535         TVdisc.Remove(i--);
1536
1537     if(TVdisc.Length())
1538     {
1539       VArr = new TColStd_HArray1OfReal(1, TVdisc.Length());
1540       for(i = 1; i <= VArr->Length(); i++)
1541         VArr->ChangeValue(i) = TVdisc(i);
1542     }
1543
1544     // fusion
1545     TColStd_SequenceOfReal Fusion;
1546     if(!CArr.IsNull())
1547     {
1548       GeomLib::FuseIntervals(BArr->ChangeArray1(),
1549         CArr->ChangeArray1(),
1550         Fusion, Precision::PConfusion());
1551       BArr = new TColStd_HArray1OfReal(1, Fusion.Length());
1552       for(i = 1; i <= BArr->Length(); i++)
1553         BArr->ChangeValue(i) = Fusion(i);
1554       Fusion.Clear();
1555     }
1556
1557     if(!UArr.IsNull())
1558     {
1559       GeomLib::FuseIntervals(BArr->ChangeArray1(),
1560         UArr->ChangeArray1(),
1561         Fusion, Precision::PConfusion());
1562       BArr = new TColStd_HArray1OfReal(1, Fusion.Length());
1563       for(i = 1; i <= BArr->Length(); i++)
1564         BArr->ChangeValue(i) = Fusion(i);
1565       Fusion.Clear();
1566     }
1567
1568     if(!VArr.IsNull())
1569     {
1570       GeomLib::FuseIntervals(BArr->ChangeArray1(),
1571         VArr->ChangeArray1(),
1572         Fusion, Precision::PConfusion());
1573       BArr = new TColStd_HArray1OfReal(1, Fusion.Length());
1574       for(i = 1; i <= BArr->Length(); i++)
1575         BArr->ChangeValue(i) = Fusion(i);
1576     }
1577
1578     const_cast<ProjLib_CompProjectedCurve*>(this)->myTabInt = new TColStd_HArray1OfReal(1, BArr->Length());
1579     for(i = 1; i <= BArr->Length(); i++)
1580       myTabInt->ChangeValue(i) = BArr->Value(i);
1581
1582 }
1583
1584 //=======================================================================
1585 //function : Trim
1586 //purpose  :
1587 //=======================================================================
1588
1589 Handle(Adaptor2d_HCurve2d) ProjLib_CompProjectedCurve::Trim
1590   (const Standard_Real First,
1591   const Standard_Real Last,
1592   const Standard_Real Tol) const
1593 {
1594   Handle(ProjLib_HCompProjectedCurve) HCS =
1595     new ProjLib_HCompProjectedCurve(*this);
1596   HCS->ChangeCurve2d().Load(mySurface);
1597   HCS->ChangeCurve2d().Load(myCurve->Trim(First,Last,Tol));
1598   return HCS;
1599 }
1600
1601 //=======================================================================
1602 //function : GetType
1603 //purpose  :
1604 //=======================================================================
1605
1606 GeomAbs_CurveType ProjLib_CompProjectedCurve::GetType() const
1607 {
1608   return GeomAbs_OtherCurve;
1609 }
1610
1611 //=======================================================================
1612 //function : UpdateTripleByTrapCriteria
1613 //purpose  :
1614 //=======================================================================
1615 void ProjLib_CompProjectedCurve::UpdateTripleByTrapCriteria(gp_Pnt &thePoint) const
1616 {
1617   Standard_Boolean isProblemsPossible = Standard_False;
1618   // Check possible traps cases:
1619
1620   // 25892 bug.
1621   if (mySurface->GetType() == GeomAbs_SurfaceOfRevolution)
1622   {
1623     // Compute maximal deviation from 3D and choose the biggest one.
1624     Standard_Real aVRes = mySurface->VResolution(Precision::Confusion());
1625     Standard_Real aMaxTol = Max(Precision::PConfusion(), aVRes);
1626
1627     if (Abs (thePoint.Z() - mySurface->FirstVParameter()) < aMaxTol ||
1628         Abs (thePoint.Z() - mySurface->LastVParameter() ) < aMaxTol )
1629     {
1630       isProblemsPossible = Standard_True;
1631     }
1632   }
1633
1634   // 27135 bug. Trap on degenerated edge.
1635   if (mySurface->GetType() == GeomAbs_Sphere &&
1636      (Abs (thePoint.Z() - mySurface->FirstVParameter()) < Precision::PConfusion() ||
1637       Abs (thePoint.Z() - mySurface->LastVParameter() ) < Precision::PConfusion() ||
1638       Abs (thePoint.Y() - mySurface->FirstUParameter()) < Precision::PConfusion() ||
1639       Abs (thePoint.Y() - mySurface->LastUParameter() ) < Precision::PConfusion() ))
1640   {
1641     isProblemsPossible = Standard_True;
1642   }
1643
1644   if (!isProblemsPossible)
1645     return;
1646
1647   Standard_Real U,V;
1648   InitialPoint(myCurve->Value(thePoint.X()), thePoint.X(), myCurve, mySurface, myTolU, myTolV, U, V);
1649
1650   // Restore original position in case of period jump.
1651   if (mySurface->IsUPeriodic() &&
1652       Abs (Abs(U - thePoint.Y()) - mySurface->UPeriod()) < Precision::PConfusion())
1653   {
1654     U = thePoint.Y();
1655   }
1656   if (mySurface->IsVPeriodic() &&
1657       Abs (Abs(V - thePoint.Z()) - mySurface->VPeriod()) < Precision::PConfusion())
1658   {
1659     V = thePoint.Z();
1660   }
1661   thePoint.SetY(U);
1662   thePoint.SetZ(V);
1663 }