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 OCCT_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 OCCT_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 OCCT_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 : myNbCurves(0),
 500   myTolU    (0.0),
 501   myTolV    (0.0),
 502   myMaxDist (0.0)
 503 {
 504 }
 505
 506 //=======================================================================
 507 //function : ProjLib_CompProjectedCurve
 508 //purpose  :
 509 //=======================================================================
 510
 511 ProjLib_CompProjectedCurve::ProjLib_CompProjectedCurve
 512                            (const Handle(Adaptor3d_HSurface)& theSurface,
 513                             const Handle(Adaptor3d_HCurve)&   theCurve,
 514                             const Standard_Real               theTolU,
 515                             const Standard_Real               theTolV)
 516 : mySurface (theSurface),
 517   myCurve   (theCurve),
 518   myNbCurves(0),
 519   mySequence(new ProjLib_HSequenceOfHSequenceOfPnt()),
 520   myTolU    (theTolU),
 521   myTolV    (theTolV),
 522   myMaxDist (-1.0)
 523 {
 524   Init();
 525 }
 526
 527 //=======================================================================
 528 //function : ProjLib_CompProjectedCurve
 529 //purpose  :
 530 //=======================================================================
 531
 532 ProjLib_CompProjectedCurve::ProjLib_CompProjectedCurve
 533                            (const Handle(Adaptor3d_HSurface)& theSurface,
 534                             const Handle(Adaptor3d_HCurve)&   theCurve,
 535                             const Standard_Real               theTolU,
 536                             const Standard_Real               theTolV,
 537                             const Standard_Real               theMaxDist)
 538 : mySurface (theSurface),
 539   myCurve   (theCurve),
 540   myNbCurves(0),
 541   mySequence(new ProjLib_HSequenceOfHSequenceOfPnt()),
 542   myTolU    (theTolU),
 543   myTolV    (theTolV),
 544   myMaxDist (theMaxDist)
 545 {
 546   Init();
 547 }
 548
 549 //=======================================================================
 550 //function : Init
 551 //purpose  :
 552 //=======================================================================
 553
 554 void ProjLib_CompProjectedCurve::Init()
 555 {
 556   myTabInt.Nullify();
 557
 558   Standard_Real Tol;// Tolerance for ExactBound
 559   Standard_Integer i, Nend = 0;
 560   Standard_Boolean FromLastU=Standard_False;
 561
 562   //new part (to discard far solutions)
 563   Standard_Real TolC = Precision::Confusion(), TolS = Precision::Confusion();
 564   Extrema_ExtCS CExt(myCurve->Curve(),
 565     mySurface->Surface(),
 566     TolC,
 567     TolS);
 568   if (CExt.IsDone() && CExt.NbExt())
 569   {
 570     // Search for the minimum solution
 571     Nend = CExt.NbExt();
 572     if(myMaxDist > 0 &&
 573        // Avoid usage of extrema result that can be wrong for extrusion
 574        mySurface->GetType() != GeomAbs_SurfaceOfExtrusion)
 575     {
 576       Standard_Real min_val2;
 577       min_val2 = CExt.SquareDistance(1);
 578       for(i = 2; i <= Nend; i++)
 579         if (CExt.SquareDistance(i) < min_val2) min_val2 = CExt.SquareDistance(i);
 580       if (min_val2 > myMaxDist * myMaxDist)
 581         return;
 582     }
 583   }
 584   // end of new part
 585
 586   Standard_Real FirstU, LastU, Step, SearchStep, WalkStep, t;
 587
 588   FirstU = myCurve->FirstParameter();
 589   LastU  = myCurve->LastParameter();
 590   const Standard_Real GlobalMinStep = 1.e-4;
 591   //<GlobalMinStep> is sufficiently small to provide solving from initial point
 592   //and, on the other hand, it is sufficiently large to avoid too close solutions.
 593   const Standard_Real MinStep = 0.01*(LastU - FirstU),
 594     MaxStep = 0.1*(LastU - FirstU);
 595   SearchStep = 10*MinStep;
 596   Step = SearchStep;
 597
 598   //Initialization of aPrjPS
 599   Standard_Real Uinf = mySurface->FirstUParameter();
 600   Standard_Real Usup = mySurface->LastUParameter();
 601   Standard_Real Vinf = mySurface->FirstVParameter();
 602   Standard_Real Vsup = mySurface->LastVParameter();
 603
 604   ProjLib_PrjResolve aPrjPS(myCurve->Curve(), mySurface->Surface(), 1);
 605
 606   t = FirstU;
 607   Standard_Boolean new_part;
 608   Standard_Real prevDeb=0.;
 609   Standard_Boolean SameDeb=Standard_False;
 610
 611
 612   gp_Pnt Triple, prevTriple;
 613
 614   //Basic loop
 615   while(t <= LastU)
 616   {
 617     //Search for the begining a new continuous part
 618     //To avoid infinite computation in some difficult cases
 619     new_part = Standard_False;
 620     if(t > FirstU && Abs(t-prevDeb) <= Precision::PConfusion()) SameDeb=Standard_True;
 621     while(t <= LastU && !new_part && !FromLastU && !SameDeb)
 622     {
 623       prevDeb=t;
 624       if (t == LastU) FromLastU=Standard_True;
 625       Standard_Boolean initpoint=Standard_False;
 626       Standard_Real U = 0., V = 0.;
 627       gp_Pnt CPoint;
 628       Standard_Real ParT,ParU,ParV;
 629
 630       // Search an initpoint in the list of Extrema Curve-Surface
 631       if(Nend != 0 && !CExt.IsParallel())
 632       {
 633         for (i=1;i<=Nend;i++)
 634         {
 635           Extrema_POnCurv P1;
 636           Extrema_POnSurf P2;
 637           CExt.Points(i,P1,P2);
 638           ParT=P1.Parameter();
 639           P2.Parameter(ParU, ParV);
 640
 641           aPrjPS.Perform(ParT, ParU, ParV, gp_Pnt2d(myTolU, myTolV),
 642             gp_Pnt2d(mySurface->FirstUParameter(),mySurface->FirstVParameter()),
 643             gp_Pnt2d(mySurface->LastUParameter(), mySurface->LastVParameter()),
 644             FuncTol, Standard_True);
 645           if ( aPrjPS.IsDone() && P1.Parameter() > Max(FirstU,t-Step+Precision::PConfusion())
 646             && P1.Parameter() <= t)
 647           {
 648             t=ParT;
 649             U=ParU;
 650             V=ParV;
 651             CPoint=P1.Value();
 652             initpoint = Standard_True;
 653             break;
 654           }
 655         }
 656       }
 657       if (!initpoint)
 658       {
 659         myCurve->D0(t,CPoint);
 660 #ifdef OCCT_DEBUG_CHRONO
 661         InitChron(chr_init_point);
 662 #endif
 663         initpoint=InitialPoint(CPoint, t,myCurve,mySurface, myTolU, myTolV, U, V);
 664 #ifdef OCCT_DEBUG_CHRONO
 665         ResultChron(chr_init_point,t_init_point);
 666         init_point_count++;
 667 #endif
 668       }
 669       if(initpoint)
 670       {
 671         // When U or V lie on surface joint in some cases we cannot use them
 672         // as initial point for aPrjPS, so we switch them
 673         gp_Vec2d D;
 674
 675         if ((mySurface->IsUPeriodic() &&
 676             Abs(Usup - Uinf - mySurface->UPeriod()) < Precision::Confusion()) ||
 677             (mySurface->IsVPeriodic() &&
 678             Abs(Vsup - Vinf - mySurface->VPeriod()) < Precision::Confusion()))
 679         {
 680           if((Abs(U - Uinf) < mySurface->UResolution(Precision::PConfusion())) &&
 681             mySurface->IsUPeriodic())
 682           {
 683             d1(t, U, V, D, myCurve, mySurface);
 684             if (D.X() < 0 ) U = Usup;
 685           }
 686           else if((Abs(U - Usup) < mySurface->UResolution(Precision::PConfusion())) &&
 687             mySurface->IsUPeriodic())
 688           {
 689             d1(t, U, V, D, myCurve, mySurface);
 690             if (D.X() > 0) U = Uinf;
 691           }
 692
 693           if((Abs(V - Vinf) < mySurface->VResolution(Precision::PConfusion())) &&
 694             mySurface->IsVPeriodic())
 695           {
 696             d1(t, U, V, D, myCurve, mySurface);
 697             if (D.Y() < 0) V = Vsup;
 698           }
 699           else if((Abs(V - Vsup) <= mySurface->VResolution(Precision::PConfusion())) &&
 700             mySurface->IsVPeriodic())
 701           {
 702             d1(t, U, V, D, myCurve, mySurface);
 703             if (D.Y() > 0) V = Vinf;
 704           }
 705         }
 706
 707
 708         if (myMaxDist > 0)
 709         {
 710           // Here we are going to stop if the distance between projection and
 711           // corresponding curve point is greater than myMaxDist
 712           gp_Pnt POnS;
 713           Standard_Real d;
 714           mySurface->D0(U, V, POnS);
 715           d = CPoint.Distance(POnS);
 716           if (d > myMaxDist)
 717           {
 718             mySequence->Clear();
 719             myNbCurves = 0;
 720             return;
 721           }
 722         }
 723         Triple = gp_Pnt(t, U, V);
 724         if (t != FirstU)
 725         {
 726           //Search for exact boundary point
 727           Tol = Min(myTolU, myTolV);
 728           gp_Vec2d D;
 729           d1(Triple.X(), Triple.Y(), Triple.Z(), D, myCurve, mySurface);
 730           Tol /= Max(Abs(D.X()), Abs(D.Y()));
 731
 732           if(!ExactBound(Triple, t - Step, Tol,
 733             myTolU, myTolV, myCurve, mySurface))
 734           {
 735 #ifdef OCCT_DEBUG
 736             cout<<"There is a problem with ExactBound computation"<<endl;
 737 #endif
 738             DichExactBound(Triple, t - Step, Tol, myTolU, myTolV,
 739               myCurve, mySurface);
 740           }
 741         }
 742         new_part = Standard_True;
 743       }
 744       else
 745       {
 746         if(t == LastU) break;
 747         t += Step;
 748         if(t>LastU)
 749         {
 750           Step =Step+LastU-t;
 751           t=LastU;
 752         }
 753       }
 754     }
 755     if (!new_part) break;
 756
 757
 758     //We have found a new continuous part
 759     Handle(TColgp_HSequenceOfPnt) hSeq = new TColgp_HSequenceOfPnt();
 760     mySequence->Append(hSeq);
 761     myNbCurves++;
 762     mySequence->Value(myNbCurves)->Append(Triple);
 763     prevTriple = Triple;
 764
 765     if (Triple.X() == LastU) break;//return;
 766
 767     //Computation of WalkStep
 768     gp_Vec D1, D2;
 769     Standard_Real MagnD1, MagnD2;
 770     d2CurvOnSurf(Triple.X(), Triple.Y(), Triple.Z(), D1, D2, myCurve, mySurface);
 771     MagnD1 = D1.Magnitude();
 772     MagnD2 = D2.Magnitude();
 773     if(MagnD2 < Precision::Confusion()) WalkStep = MaxStep;
 774     else WalkStep = Min(MaxStep, Max(MinStep, 0.1*MagnD1/MagnD2));
 775
 776     Step = WalkStep;
 777
 778     t = Triple.X() + Step;
 779     if (t > LastU) t = LastU;
 780     Standard_Real prevStep = Step;
 781     Standard_Real U0, V0;
 782     gp_Pnt2d aLowBorder(mySurface->FirstUParameter(),mySurface->FirstVParameter());
 783     gp_Pnt2d aUppBorder(mySurface->LastUParameter(), mySurface->LastVParameter());
 784     gp_Pnt2d aTol(myTolU, myTolV);
 785     //Here we are trying to prolong continuous part
 786     while (t <= LastU && new_part)
 787     {
 788
 789       U0 = Triple.Y() + (Step / prevStep) * (Triple.Y() - prevTriple.Y());
 790       V0 = Triple.Z() + (Step / prevStep) * (Triple.Z() - prevTriple.Z());
 791       // adjust U0 to be in [mySurface->FirstUParameter(),mySurface->LastUParameter()]
 792       U0 = Min(Max(U0, aLowBorder.X()), aUppBorder.X());
 793       // adjust V0 to be in [mySurface->FirstVParameter(),mySurface->LastVParameter()]
 794       V0 = Min(Max(V0, aLowBorder.Y()), aUppBorder.Y());
 795
 796
 797       aPrjPS.Perform(t, U0, V0, aTol,
 798                      aLowBorder, aUppBorder, FuncTol, Standard_True);
 799       if(!aPrjPS.IsDone())
 800       {
 801         if (Step <= GlobalMinStep)
 802         {
 803           //Search for exact boundary point
 804           Tol = Min(myTolU, myTolV);
 805           gp_Vec2d D;
 806           d1(Triple.X(), Triple.Y(), Triple.Z(), D, myCurve, mySurface);
 807           Tol /= Max(Abs(D.X()), Abs(D.Y()));
 808
 809           if(!ExactBound(Triple, t, Tol, myTolU, myTolV,
 810             myCurve, mySurface))
 811           {
 812 #ifdef OCCT_DEBUG
 813             cout<<"There is a problem with ExactBound computation"<<endl;
 814 #endif
 815             DichExactBound(Triple, t, Tol, myTolU, myTolV,
 816               myCurve, mySurface);
 817           }
 818
 819           if((Triple.X() - mySequence->Value(myNbCurves)->Value(mySequence->Value(myNbCurves)->Length()).X()) > 1.e-10)
 820             mySequence->Value(myNbCurves)->Append(Triple);
 821           if((LastU - Triple.X()) < Tol) {t = LastU + 1; break;}//return;
 822
 823           Step = SearchStep;
 824           t = Triple.X() + Step;
 825           if (t > (LastU-MinStep/2) )
 826           {
 827             Step =Step+LastU-t;
 828             t = LastU;
 829           }
 830           new_part = Standard_False;
 831         }
 832         else
 833         {
 834           // decrease step
 835           Standard_Real SaveStep = Step;
 836           Step /= 2.;
 837           t = Triple .X() + Step;
 838           if (t > (LastU-MinStep/4) )
 839           {
 840             Step =Step+LastU-t;
 841             if (Abs(Step - SaveStep) <= Precision::PConfusion())
 842               Step = GlobalMinStep; //to avoid looping
 843             t = LastU;
 844           }
 845         }
 846       }
 847       // Go further
 848       else
 849       {
 850         prevTriple = Triple;
 851         prevStep = Step;
 852         Triple = gp_Pnt(t, aPrjPS.Solution().X(), aPrjPS.Solution().Y());
 853
 854         if (mySurface->GetType() == GeomAbs_SurfaceOfRevolution &&
 855            (Abs (Triple.Z() - mySurface->FirstVParameter()) < Precision::Confusion() ||
 856             Abs (Triple.Z() - mySurface->LastVParameter() ) < Precision::Confusion() ))
 857         {
 858           // Go out from possible attraktor.
 859
 860           Standard_Real U,V;
 861           InitialPoint(myCurve->Value(t), t, myCurve, mySurface, myTolU, myTolV, U, V);
 862           if (Abs (Abs(U - Triple.Y()) - mySurface->UPeriod()) < Precision::Confusion())
 863           {
 864             // Handle period jump.
 865             U = Triple.Y();
 866           }
 867           Triple.SetY(U);
 868           Triple.SetZ(V);
 869         }
 870
 871         if((Triple.X() - mySequence->Value(myNbCurves)->Value(mySequence->Value(myNbCurves)->Length()).X()) > 1.e-10)
 872           mySequence->Value(myNbCurves)->Append(Triple);
 873         if (t == LastU) {t = LastU + 1; break;}//return;
 874
 875         //Computation of WalkStep
 876         d2CurvOnSurf(Triple.X(), Triple.Y(), Triple.Z(), D1, D2, myCurve, mySurface);
 877         MagnD1 = D1.Magnitude();
 878         MagnD2 = D2.Magnitude();
 879         if(MagnD2 < Precision::Confusion() ) WalkStep = MaxStep;
 880         else WalkStep = Min(MaxStep, Max(MinStep, 0.1*MagnD1/MagnD2));
 881
 882         Step = WalkStep;
 883         t += Step;
 884         if (t > (LastU-MinStep/2) )
 885         {
 886           Step =Step+LastU-t;
 887           t = LastU;
 888         }
 889       }
 890     }
 891   }
 892   // Sequence postproceeding
 893   Standard_Integer j;
 894
 895   // 1. Removing poor parts
 896   Standard_Integer NbPart=myNbCurves;
 897   Standard_Integer ipart=1;
 898   for(i = 1; i <= NbPart; i++) {
 899     //    Standard_Integer NbPoints = mySequence->Value(i)->Length();
 900     if(mySequence->Value(ipart)->Length() < 2) {
 901       mySequence->Remove(ipart);
 902       myNbCurves--;
 903     }
 904     else ipart++;
 905   }
 906
 907   if(myNbCurves == 0) return;
 908
 909   // 2. Removing common parts of bounds
 910   for(i = 1; i < myNbCurves; i++)
 911   {
 912     if(mySequence->Value(i)->Value(mySequence->Value(i)->Length()).X() >=
 913       mySequence->Value(i+1)->Value(1).X())
 914       mySequence->ChangeValue(i+1)->ChangeValue(1).SetX(mySequence->Value(i)->Value(mySequence->Value(i)->Length()).X() + 1.e-12);
 915   }
 916
 917   // 3. Computation of the maximum distance from each part of curve to surface
 918
 919   myMaxDistance = new TColStd_HArray1OfReal(1, myNbCurves);
 920   myMaxDistance->Init(0);
 921   for(i = 1; i <= myNbCurves; i++)
 922     for(j = 1; j <= mySequence->Value(i)->Length(); j++)
 923     {
 924       gp_Pnt POnC, POnS, Triple;
 925       Standard_Real Distance;
 926       Triple = mySequence->Value(i)->Value(j);
 927       myCurve->D0(Triple.X(), POnC);
 928       mySurface->D0(Triple.Y(), Triple.Z(), POnS);
 929       Distance = POnC.Distance(POnS);
 930       if (myMaxDistance->Value(i) < Distance)
 931         myMaxDistance->ChangeValue(i) = Distance;
 932     }
 933
 934
 935     // 4. Check the projection to be a single point
 936
 937     gp_Pnt2d Pmoy, Pcurr, P;
 938     Standard_Real AveU, AveV;
 939     mySnglPnts = new TColStd_HArray1OfBoolean(1, myNbCurves);
 940     for(i = 1; i <= myNbCurves; i++) mySnglPnts->SetValue(i, Standard_True);
 941
 942     for(i = 1; i <= myNbCurves; i++)
 943     {
 944       //compute an average U and V
 945
 946       for(j = 1, AveU = 0., AveV = 0.; j <= mySequence->Value(i)->Length(); j++)
 947       {
 948         AveU += mySequence->Value(i)->Value(j).Y();
 949         AveV += mySequence->Value(i)->Value(j).Z();
 950       }
 951       AveU /= mySequence->Value(i)->Length();
 952       AveV /= mySequence->Value(i)->Length();
 953
 954       Pmoy.SetCoord(AveU,AveV);
 955       for(j = 1; j <= mySequence->Value(i)->Length(); j++)
 956       {
 957         Pcurr =
 958           gp_Pnt2d(mySequence->Value(i)->Value(j).Y(), mySequence->Value(i)->Value(j).Z());
 959         if (Pcurr.Distance(Pmoy) > ((myTolU < myTolV) ? myTolV : myTolU))
 960         {
 961           mySnglPnts->SetValue(i, Standard_False);
 962           break;
 963         }
 964       }
 965     }
 966
 967     // 5. Check the projection to be an isoparametric curve of the surface
 968
 969     myUIso = new TColStd_HArray1OfBoolean(1, myNbCurves);
 970     for(i = 1; i <= myNbCurves; i++) myUIso->SetValue(i, Standard_True);
 971
 972     myVIso = new TColStd_HArray1OfBoolean(1, myNbCurves);
 973     for(i = 1; i <= myNbCurves; i++) myVIso->SetValue(i, Standard_True);
 974
 975     for(i = 1; i <= myNbCurves; i++) {
 976       if (IsSinglePnt(i, P)|| mySequence->Value(i)->Length() <=2) {
 977         myUIso->SetValue(i, Standard_False);
 978         myVIso->SetValue(i, Standard_False);
 979         continue;
 980       }
 981
 982       // new test for isoparametrics
 983
 984       if ( mySequence->Value(i)->Length() > 2) {
 985         //compute an average U and V
 986
 987         for(j = 1, AveU = 0., AveV = 0.; j <= mySequence->Value(i)->Length(); j++) {
 988           AveU += mySequence->Value(i)->Value(j).Y();
 989           AveV += mySequence->Value(i)->Value(j).Z();
 990         }
 991         AveU /= mySequence->Value(i)->Length();
 992         AveV /= mySequence->Value(i)->Length();
 993
 994         // is i-part U-isoparametric ?
 995         for(j = 1; j <= mySequence->Value(i)->Length(); j++)
 996         {
 997           if(Abs(mySequence->Value(i)->Value(j).Y() - AveU) > myTolU)
 998           {
 999             myUIso->SetValue(i, Standard_False);
1000             break;
1001           }
1002         }
1003
1004         // is i-part V-isoparametric ?
1005         for(j = 1; j <= mySequence->Value(i)->Length(); j++)
1006         {
1007           if(Abs(mySequence->Value(i)->Value(j).Z() - AveV) > myTolV)
1008           {
1009             myVIso->SetValue(i, Standard_False);
1010             break;
1011           }
1012         }
1013         //
1014       }
1015     }
1016 }
1017 //=======================================================================
1018 //function : Load
1019 //purpose  :
1020 //=======================================================================
1021
1022 void ProjLib_CompProjectedCurve::Load(const Handle(Adaptor3d_HSurface)& S)
1023 {
1024   mySurface = S;
1025 }
1026
1027 //=======================================================================
1028 //function : Load
1029 //purpose  :
1030 //=======================================================================
1031
1032 void ProjLib_CompProjectedCurve::Load(const Handle(Adaptor3d_HCurve)& C)
1033 {
1034   myCurve = C;
1035 }
1036
1037 //=======================================================================
1038 //function : GetSurface
1039 //purpose  :
1040 //=======================================================================
1041
1042 const Handle(Adaptor3d_HSurface)& ProjLib_CompProjectedCurve::GetSurface() const
1043 {
1044   return mySurface;
1045 }
1046
1047
1048 //=======================================================================
1049 //function : GetCurve
1050 //purpose  :
1051 //=======================================================================
1052
1053 const Handle(Adaptor3d_HCurve)& ProjLib_CompProjectedCurve::GetCurve() const
1054 {
1055   return myCurve;
1056 }
1057
1058 //=======================================================================
1059 //function : GetTolerance
1060 //purpose  :
1061 //=======================================================================
1062
1063 void ProjLib_CompProjectedCurve::GetTolerance(Standard_Real& TolU,
1064   Standard_Real& TolV) const
1065 {
1066   TolU = myTolU;
1067   TolV = myTolV;
1068 }
1069
1070 //=======================================================================
1071 //function : NbCurves
1072 //purpose  :
1073 //=======================================================================
1074
1075 Standard_Integer ProjLib_CompProjectedCurve::NbCurves() const
1076 {
1077   return myNbCurves;
1078 }
1079 //=======================================================================
1080 //function : Bounds
1081 //purpose  :
1082 //=======================================================================
1083
1084 void ProjLib_CompProjectedCurve::Bounds(const Standard_Integer Index,
1085   Standard_Real& Udeb,
1086   Standard_Real& Ufin) const
1087 {
1088   if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1089   Udeb = mySequence->Value(Index)->Value(1).X();
1090   Ufin = mySequence->Value(Index)->Value(mySequence->Value(Index)->Length()).X();
1091 }
1092 //=======================================================================
1093 //function : IsSinglePnt
1094 //purpose  :
1095 //=======================================================================
1096
1097 Standard_Boolean ProjLib_CompProjectedCurve::IsSinglePnt(const Standard_Integer Index, gp_Pnt2d& P) const
1098 {
1099   if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1100   P = gp_Pnt2d(mySequence->Value(Index)->Value(1).Y(), mySequence->Value(Index)->Value(1).Z());
1101   return mySnglPnts->Value(Index);
1102 }
1103
1104 //=======================================================================
1105 //function : IsUIso
1106 //purpose  :
1107 //=======================================================================
1108
1109 Standard_Boolean ProjLib_CompProjectedCurve::IsUIso(const Standard_Integer Index, Standard_Real& U) const
1110 {
1111   if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1112   U = mySequence->Value(Index)->Value(1).Y();
1113   return myUIso->Value(Index);
1114 }
1115 //=======================================================================
1116 //function : IsVIso
1117 //purpose  :
1118 //=======================================================================
1119
1120 Standard_Boolean ProjLib_CompProjectedCurve::IsVIso(const Standard_Integer Index, Standard_Real& V) const
1121 {
1122   if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1123   V = mySequence->Value(Index)->Value(1).Z();
1124   return myVIso->Value(Index);
1125 }
1126 //=======================================================================
1127 //function : Value
1128 //purpose  :
1129 //=======================================================================
1130
1131 gp_Pnt2d ProjLib_CompProjectedCurve::Value(const Standard_Real t) const
1132 {
1133   gp_Pnt2d P;
1134   D0(t, P);
1135   return P;
1136 }
1137 //=======================================================================
1138 //function : D0
1139 //purpose  :
1140 //=======================================================================
1141
1142 void ProjLib_CompProjectedCurve::D0(const Standard_Real U,gp_Pnt2d& P) const
1143 {
1144   Standard_Integer i, j;
1145   Standard_Real Udeb, Ufin;
1146   Standard_Boolean found = Standard_False;
1147
1148   for(i = 1; i <= myNbCurves; i++)
1149   {
1150     Bounds(i, Udeb, Ufin);
1151     if (U >= Udeb && U <= Ufin)
1152     {
1153       found = Standard_True;
1154       break;
1155     }
1156   }
1157   if (!found) Standard_DomainError::Raise("ProjLib_CompProjectedCurve::D0");
1158
1159   Standard_Real U0, V0;
1160
1161   Standard_Integer End = mySequence->Value(i)->Length();
1162   for(j = 1; j < End; j++)
1163     if ((U >= mySequence->Value(i)->Value(j).X()) && (U <= mySequence->Value(i)->Value(j + 1).X())) break;
1164
1165   //  U0 = mySequence->Value(i)->Value(j).Y();
1166   //  V0 = mySequence->Value(i)->Value(j).Z();
1167
1168   //  Cubic Interpolation
1169   if(mySequence->Value(i)->Length() < 4 ||
1170     (Abs(U-mySequence->Value(i)->Value(j).X()) <= Precision::PConfusion()) )
1171   {
1172     U0 = mySequence->Value(i)->Value(j).Y();
1173     V0 = mySequence->Value(i)->Value(j).Z();
1174   }
1175   else if (Abs(U-mySequence->Value(i)->Value(j+1).X())
1176     <= Precision::PConfusion())
1177   {
1178     U0 = mySequence->Value(i)->Value(j+1).Y();
1179     V0 = mySequence->Value(i)->Value(j+1).Z();
1180   }
1181   else
1182   {
1183     if (j == 1) j = 2;
1184     if (j > mySequence->Value(i)->Length() - 2)
1185       j = mySequence->Value(i)->Length() - 2;
1186
1187     gp_Vec2d I1, I2, I3, I21, I22, I31, Y1, Y2, Y3, Y4, Res;
1188     Standard_Real X1, X2, X3, X4;
1189
1190     X1 = mySequence->Value(i)->Value(j - 1).X();
1191     X2 = mySequence->Value(i)->Value(j).X();
1192     X3 = mySequence->Value(i)->Value(j + 1).X();
1193     X4 = mySequence->Value(i)->Value(j + 2).X();
1194
1195     Y1 = gp_Vec2d(mySequence->Value(i)->Value(j - 1).Y(),
1196       mySequence->Value(i)->Value(j - 1).Z());
1197     Y2 = gp_Vec2d(mySequence->Value(i)->Value(j).Y(),
1198       mySequence->Value(i)->Value(j).Z());
1199     Y3 = gp_Vec2d(mySequence->Value(i)->Value(j + 1).Y(),
1200       mySequence->Value(i)->Value(j + 1).Z());
1201     Y4 = gp_Vec2d(mySequence->Value(i)->Value(j + 2).Y(),
1202       mySequence->Value(i)->Value(j + 2).Z());
1203
1204     I1 = (Y1 - Y2)/(X1 - X2);
1205     I2 = (Y2 - Y3)/(X2 - X3);
1206     I3 = (Y3 - Y4)/(X3 - X4);
1207
1208     I21 = (I1 - I2)/(X1 - X3);
1209     I22 = (I2 - I3)/(X2 - X4);
1210
1211     I31 = (I21 - I22)/(X1 - X4);
1212
1213     Res = Y1 + (U - X1)*(I1 + (U - X2)*(I21 + (U - X3)*I31));
1214
1215     U0 = Res.X();
1216     V0 = Res.Y();
1217
1218     if(U0 < mySurface->FirstUParameter()) U0 = mySurface->FirstUParameter();
1219     else if(U0 > mySurface->LastUParameter()) U0 = mySurface->LastUParameter();
1220
1221     if(V0 < mySurface->FirstVParameter()) V0 = mySurface->FirstVParameter();
1222     else if(V0 > mySurface->LastVParameter()) V0 = mySurface->LastVParameter();
1223   }
1224   //End of cubic interpolation
1225
1226   ProjLib_PrjResolve aPrjPS(myCurve->Curve(), mySurface->Surface(), 1);
1227   aPrjPS.Perform(U, U0, V0, gp_Pnt2d(myTolU, myTolV),
1228     gp_Pnt2d(mySurface->FirstUParameter(), mySurface->FirstVParameter()),
1229     gp_Pnt2d(mySurface->LastUParameter(), mySurface->LastVParameter()));
1230   if (aPrjPS.IsDone())
1231     P = aPrjPS.Solution();
1232   else
1233   {
1234     gp_Pnt thePoint = myCurve->Value(U);
1235     Extrema_ExtPS aExtPS(thePoint, mySurface->Surface(), myTolU, myTolV);
1236     if (aExtPS.IsDone() && aExtPS.NbExt())
1237     {
1238       Standard_Integer i, Nend, imin = 1;
1239       // Search for the nearest solution which is also a normal projection
1240       Nend = aExtPS.NbExt();
1241       for(i = 2; i <= Nend; i++)
1242         if (aExtPS.SquareDistance(i) < aExtPS.SquareDistance(imin))
1243           imin = i;
1244       const Extrema_POnSurf& POnS = aExtPS.Point(imin);
1245       Standard_Real ParU,ParV;
1246       POnS.Parameter(ParU, ParV);
1247       P.SetCoord(ParU, ParV);
1248     }
1249     else
1250       P.SetCoord(U0,V0);
1251   }
1252 }
1253 //=======================================================================
1254 //function : D1
1255 //purpose  :
1256 //=======================================================================
1257
1258 void ProjLib_CompProjectedCurve::D1(const Standard_Real t,
1259   gp_Pnt2d& P,
1260   gp_Vec2d& V) const
1261 {
1262   Standard_Real u, v;
1263   D0(t, P);
1264   u = P.X();
1265   v = P.Y();
1266   d1(t, u, v, V, myCurve, mySurface);
1267 }
1268 //=======================================================================
1269 //function : D2
1270 //purpose  :
1271 //=======================================================================
1272
1273 void ProjLib_CompProjectedCurve::D2(const Standard_Real t,
1274   gp_Pnt2d& P,
1275   gp_Vec2d& V1,
1276   gp_Vec2d& V2) const
1277 {
1278   Standard_Real u, v;
1279   D0(t, P);
1280   u = P.X();
1281   v = P.Y();
1282   d2(t, u, v, V1, V2, myCurve, mySurface);
1283 }
1284 //=======================================================================
1285 //function : DN
1286 //purpose  :
1287 //=======================================================================
1288
1289 gp_Vec2d ProjLib_CompProjectedCurve::DN(const Standard_Real t,
1290   const Standard_Integer N) const
1291 {
1292   if (N < 1 ) Standard_OutOfRange::Raise("ProjLib_CompProjectedCurve : N must be greater than 0");
1293   else if (N ==1)
1294   {
1295     gp_Pnt2d P;
1296     gp_Vec2d V;
1297     D1(t,P,V);
1298     return V;
1299   }
1300   else if ( N==2)
1301   {
1302     gp_Pnt2d P;
1303     gp_Vec2d V1,V2;
1304     D2(t,P,V1,V2);
1305     return V2;
1306   }
1307   else if (N > 2 )
1308     Standard_NotImplemented::Raise("ProjLib_CompProjectedCurve::DN");
1309   return gp_Vec2d();
1310 }
1311
1312 //=======================================================================
1313 //function : GetSequence
1314 //purpose  :
1315 //=======================================================================
1316
1317 const Handle(ProjLib_HSequenceOfHSequenceOfPnt)& ProjLib_CompProjectedCurve::GetSequence() const
1318 {
1319   return mySequence;
1320 }
1321 //=======================================================================
1322 //function : FirstParameter
1323 //purpose  :
1324 //=======================================================================
1325
1326 Standard_Real ProjLib_CompProjectedCurve::FirstParameter() const
1327 {
1328   return myCurve->FirstParameter();
1329 }
1330
1331 //=======================================================================
1332 //function : LastParameter
1333 //purpose  :
1334 //=======================================================================
1335
1336 Standard_Real ProjLib_CompProjectedCurve::LastParameter() const
1337 {
1338   return myCurve->LastParameter();
1339 }
1340
1341 //=======================================================================
1342 //function : MaxDistance
1343 //purpose  :
1344 //=======================================================================
1345
1346 Standard_Real ProjLib_CompProjectedCurve::MaxDistance(const Standard_Integer Index) const
1347 {
1348   if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1349   return myMaxDistance->Value(Index);
1350 }
1351
1352 //=======================================================================
1353 //function : NbIntervals
1354 //purpose  :
1355 //=======================================================================
1356
1357 Standard_Integer ProjLib_CompProjectedCurve::NbIntervals(const GeomAbs_Shape S) const
1358 {
1359   const_cast<ProjLib_CompProjectedCurve*>(this)->myTabInt.Nullify();
1360   BuildIntervals(S);
1361   return myTabInt->Length() - 1;
1362 }
1363
1364 //=======================================================================
1365 //function : Intervals
1366 //purpose  :
1367 //=======================================================================
1368
1369 void ProjLib_CompProjectedCurve::Intervals(TColStd_Array1OfReal& T,const GeomAbs_Shape S) const
1370 {
1371   if (myTabInt.IsNull()) BuildIntervals (S);
1372   T = myTabInt->Array1();
1373 }
1374
1375 //=======================================================================
1376 //function : BuildIntervals
1377 //purpose  :
1378 //=======================================================================
1379
1380 void ProjLib_CompProjectedCurve::BuildIntervals(const GeomAbs_Shape S) const
1381 {
1382   GeomAbs_Shape SforS = GeomAbs_CN;
1383   switch(S) {
1384   case GeomAbs_C0:
1385     SforS = GeomAbs_C1;
1386     break;
1387   case GeomAbs_C1:
1388     SforS = GeomAbs_C2;
1389     break;
1390   case GeomAbs_C2:
1391     SforS = GeomAbs_C3;
1392     break;
1393   case GeomAbs_C3:
1394     SforS = GeomAbs_CN;
1395     break;
1396   case GeomAbs_CN:
1397     SforS = GeomAbs_CN;
1398     break;
1399   default:
1400     Standard_OutOfRange::Raise();
1401   }
1402   Standard_Integer i, j, k;
1403   Standard_Integer NbIntCur = myCurve->NbIntervals(S);
1404   Standard_Integer NbIntSurU = mySurface->NbUIntervals(SforS);
1405   Standard_Integer NbIntSurV = mySurface->NbVIntervals(SforS);
1406
1407   TColStd_Array1OfReal CutPntsT(1, NbIntCur+1);
1408   TColStd_Array1OfReal CutPntsU(1, NbIntSurU+1);
1409   TColStd_Array1OfReal CutPntsV(1, NbIntSurV+1);
1410
1411   myCurve->Intervals(CutPntsT, S);
1412   mySurface->UIntervals(CutPntsU, SforS);
1413   mySurface->VIntervals(CutPntsV, SforS);
1414
1415   Standard_Real Tl, Tr, Ul, Ur, Vl, Vr, Tol;
1416
1417   Handle(TColStd_HArray1OfReal) BArr = NULL,
1418     CArr = NULL,
1419     UArr = NULL,
1420     VArr = NULL;
1421
1422   // proccessing projection bounds
1423   BArr = new TColStd_HArray1OfReal(1, 2*myNbCurves);
1424   for(i = 1; i <= myNbCurves; i++)
1425     Bounds(i, BArr->ChangeValue(2*i - 1), BArr->ChangeValue(2*i));
1426
1427   // proccessing curve discontinuities
1428   if(NbIntCur > 1) {
1429     CArr = new TColStd_HArray1OfReal(1, NbIntCur - 1);
1430     for(i = 1; i <= CArr->Length(); i++)
1431       CArr->ChangeValue(i) = CutPntsT(i + 1);
1432   }
1433
1434   // proccessing U-surface discontinuities
1435   TColStd_SequenceOfReal TUdisc;
1436
1437   for(k = 2; k <= NbIntSurU; k++) {
1438     //    cout<<"CutPntsU("<<k<<") = "<<CutPntsU(k)<<endl;
1439     for(i = 1; i <= myNbCurves; i++)
1440       for(j = 1; j < mySequence->Value(i)->Length(); j++) {
1441         Ul = mySequence->Value(i)->Value(j).Y();
1442         Ur = mySequence->Value(i)->Value(j + 1).Y();
1443
1444         if(Abs(Ul - CutPntsU(k)) <= myTolU)
1445           TUdisc.Append(mySequence->Value(i)->Value(j).X());
1446         else if(Abs(Ur - CutPntsU(k)) <= myTolU)
1447           TUdisc.Append(mySequence->Value(i)->Value(j + 1).X());
1448         else if((Ul < CutPntsU(k) && CutPntsU(k) < Ur) ||
1449           (Ur < CutPntsU(k) && CutPntsU(k) < Ul))
1450         {
1451           Standard_Real V;
1452           V = (mySequence->Value(i)->Value(j).Z()
1453             + mySequence->Value(i)->Value(j +1).Z())/2;
1454           ProjLib_PrjResolve Solver(myCurve->Curve(), mySurface->Surface(), 2);
1455
1456           gp_Vec2d D;
1457           gp_Pnt Triple;
1458           Triple = mySequence->Value(i)->Value(j);
1459           d1(Triple.X(), Triple.Y(), Triple.Z(), D, myCurve, mySurface);
1460           if (Abs(D.X()) < Precision::Confusion())
1461             Tol = myTolU;
1462           else
1463             Tol = Min(myTolU, myTolU / Abs(D.X()));
1464
1465           Tl = mySequence->Value(i)->Value(j).X();
1466           Tr = mySequence->Value(i)->Value(j + 1).X();
1467
1468           Solver.Perform((Tl + Tr)/2, CutPntsU(k), V,
1469             gp_Pnt2d(Tol, myTolV),
1470             gp_Pnt2d(Tl, mySurface->FirstVParameter()),
1471             gp_Pnt2d(Tr, mySurface->LastVParameter()));
1472           //
1473           if(Solver.IsDone())
1474           {
1475             TUdisc.Append(Solver.Solution().X());
1476           }
1477         }
1478       }
1479   }
1480   for(i = 2; i <= TUdisc.Length(); i++)
1481     if(TUdisc(i) - TUdisc(i-1) < Precision::PConfusion())
1482       TUdisc.Remove(i--);
1483
1484   if(TUdisc.Length())
1485   {
1486     UArr = new TColStd_HArray1OfReal(1, TUdisc.Length());
1487     for(i = 1; i <= UArr->Length(); i++)
1488       UArr->ChangeValue(i) = TUdisc(i);
1489   }
1490   // proccessing V-surface discontinuities
1491   TColStd_SequenceOfReal TVdisc;
1492
1493   for(k = 2; k <= NbIntSurV; k++)
1494     for(i = 1; i <= myNbCurves; i++)
1495     {
1496       //      cout<<"CutPntsV("<<k<<") = "<<CutPntsV(k)<<endl;
1497       for(j = 1; j < mySequence->Value(i)->Length(); j++) {
1498
1499         Vl = mySequence->Value(i)->Value(j).Z();
1500         Vr = mySequence->Value(i)->Value(j + 1).Z();
1501
1502         if(Abs(Vl - CutPntsV(k)) <= myTolV)
1503           TVdisc.Append(mySequence->Value(i)->Value(j).X());
1504         else if (Abs(Vr - CutPntsV(k)) <= myTolV)
1505           TVdisc.Append(mySequence->Value(i)->Value(j + 1).X());
1506         else if((Vl < CutPntsV(k) && CutPntsV(k) < Vr) ||
1507           (Vr < CutPntsV(k) && CutPntsV(k) < Vl))
1508         {
1509           Standard_Real U;
1510           U = (mySequence->Value(i)->Value(j).Y()
1511             + mySequence->Value(i)->Value(j +1).Y())/2;
1512           ProjLib_PrjResolve Solver(myCurve->Curve(), mySurface->Surface(), 3);
1513
1514           gp_Vec2d D;
1515           gp_Pnt Triple;
1516           Triple = mySequence->Value(i)->Value(j);
1517           d1(Triple.X(), Triple.Y(), Triple.Z(), D, myCurve, mySurface);
1518           if (Abs(D.Y()) < Precision::Confusion())
1519             Tol = myTolV;
1520           else
1521             Tol = Min(myTolV, myTolV / Abs(D.Y()));
1522
1523           Tl = mySequence->Value(i)->Value(j).X();
1524           Tr = mySequence->Value(i)->Value(j + 1).X();
1525
1526           Solver.Perform((Tl + Tr)/2, U, CutPntsV(k),
1527             gp_Pnt2d(Tol, myTolV),
1528             gp_Pnt2d(Tl, mySurface->FirstUParameter()),
1529             gp_Pnt2d(Tr, mySurface->LastUParameter()));
1530           //
1531           if(Solver.IsDone())
1532           {
1533             TVdisc.Append(Solver.Solution().X());
1534           }
1535         }
1536       }
1537     }
1538     for(i = 2; i <= TVdisc.Length(); i++)
1539       if(TVdisc(i) - TVdisc(i-1) < Precision::PConfusion())
1540         TVdisc.Remove(i--);
1541
1542     if(TVdisc.Length())
1543     {
1544       VArr = new TColStd_HArray1OfReal(1, TVdisc.Length());
1545       for(i = 1; i <= VArr->Length(); i++)
1546         VArr->ChangeValue(i) = TVdisc(i);
1547     }
1548
1549     // fusion
1550     TColStd_SequenceOfReal Fusion;
1551     if(!CArr.IsNull())
1552     {
1553       GeomLib::FuseIntervals(BArr->ChangeArray1(),
1554         CArr->ChangeArray1(),
1555         Fusion, Precision::PConfusion());
1556       BArr = new TColStd_HArray1OfReal(1, Fusion.Length());
1557       for(i = 1; i <= BArr->Length(); i++)
1558         BArr->ChangeValue(i) = Fusion(i);
1559       Fusion.Clear();
1560     }
1561
1562     if(!UArr.IsNull())
1563     {
1564       GeomLib::FuseIntervals(BArr->ChangeArray1(),
1565         UArr->ChangeArray1(),
1566         Fusion, Precision::PConfusion());
1567       BArr = new TColStd_HArray1OfReal(1, Fusion.Length());
1568       for(i = 1; i <= BArr->Length(); i++)
1569         BArr->ChangeValue(i) = Fusion(i);
1570       Fusion.Clear();
1571     }
1572
1573     if(!VArr.IsNull())
1574     {
1575       GeomLib::FuseIntervals(BArr->ChangeArray1(),
1576         VArr->ChangeArray1(),
1577         Fusion, Precision::PConfusion());
1578       BArr = new TColStd_HArray1OfReal(1, Fusion.Length());
1579       for(i = 1; i <= BArr->Length(); i++)
1580         BArr->ChangeValue(i) = Fusion(i);
1581     }
1582
1583     const_cast<ProjLib_CompProjectedCurve*>(this)->myTabInt = new TColStd_HArray1OfReal(1, BArr->Length());
1584     for(i = 1; i <= BArr->Length(); i++)
1585       myTabInt->ChangeValue(i) = BArr->Value(i);
1586
1587 }
1588
1589 //=======================================================================
1590 //function : Trim
1591 //purpose  :
1592 //=======================================================================
1593
1594 Handle(Adaptor2d_HCurve2d) ProjLib_CompProjectedCurve::Trim
1595   (const Standard_Real First,
1596   const Standard_Real Last,
1597   const Standard_Real Tol) const
1598 {
1599   Handle(ProjLib_HCompProjectedCurve) HCS =
1600     new ProjLib_HCompProjectedCurve(*this);
1601   HCS->ChangeCurve2d().Load(mySurface);
1602   HCS->ChangeCurve2d().Load(myCurve->Trim(First,Last,Tol));
1603   return HCS;
1604 }
1605
1606 //=======================================================================
1607 //function : GetType
1608 //purpose  :
1609 //=======================================================================
1610
1611 GeomAbs_CurveType ProjLib_CompProjectedCurve::GetType() const
1612 {
1613   return GeomAbs_OtherCurve;
1614 }