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