src/ShapeConstruct/ShapeConstruct_ProjectCurveOnSurface.cxx

   1 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   2 //
   3 // This file is part of Open CASCADE Technology software library.
   4 //
   5 // This library is free software; you can redistribute it and/or modify it under
   6 // the terms of the GNU Lesser General Public License version 2.1 as published
   7 // by the Free Software Foundation, with special exception defined in the file
   8 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
   9 // distribution for complete text of the license and disclaimer of any warranty.
  10 //
  11 // Alternatively, this file may be used under the terms of Open CASCADE
  12 // commercial license or contractual agreement.
  13
  14 //:k1 abv 16.12.98 K4L PRO10107, PRO10108, PRO10109
  15 //:j8 abv 10.12.98 TR10 r0501_db.stp #9423
  16 //:S4030 abv, pdn: new methods - interface to standard ProjLib_CompProjectedCurve
  17 //%12 pdn 15.02.99 PRO9234 optimizing
  18 //%12 pdn 15.02.99 PRO9234 using improved ProjectDegenerated method
  19 //    rln 03.03.99 S4135: bm2_sd_t4-A.stp treatment of Geom_SphericalSurface together with V-closed surfaces
  20 //:p9 abv 11.03.99 PRO7226 #489490: make IsAnIsoparametric to find nearest case
  21 //:q1 abv 15.03.99 (pdn) PRO7226 #525030: limit NextValueOfUV() by tolerance
  22 //:q5 abv 19.03.99 code improvement
  23 //:q9 abv 23.03.99 PRO7226.stp #489490: cashe for projecting end points
  24 //#78 rln 12.03.99 S4135: checking spatial closure with myPreci
  25 //    pdn 12.03.99 S4135: creating pcurve with minimal length in the case of densed points
  26 //    abv 29.03.99 IsAnIsoparametric with Precision::Confusion
  27 //    pdn 09.04.99 IsAnisoparametric uses already computed parameters (S4030, fix PRO14323)
  28 //szv#4 S4163
  29 //:s5 abv 22.04.99  Adding debug printouts in catch {} blocks
  30 //#1  svv 11.01.00  Porting on DEC
  31
  32 #include <Approx_CurveOnSurface.hxx>
  33 #include <Geom2d_BSplineCurve.hxx>
  34 #include <Geom2d_Circle.hxx>
  35 #include <Geom2d_Curve.hxx>
  36 #include <Geom2d_Ellipse.hxx>
  37 #include <Geom2d_Hyperbola.hxx>
  38 #include <Geom2d_Line.hxx>
  39 #include <Geom2d_Parabola.hxx>
  40 #include <Geom2d_TrimmedCurve.hxx>
  41 #include <Geom2dAdaptor.hxx>
  42 #include <Geom2dAPI_Interpolate.hxx>
  43 #include <Geom_BezierSurface.hxx>
  44 #include <Geom_BoundedCurve.hxx>
  45 #include <Geom_BSplineCurve.hxx>
  46 #include <Geom_BSplineSurface.hxx>
  47 #include <Geom_Curve.hxx>
  48 #include <Geom_OffsetSurface.hxx>
  49 #include <Geom_Plane.hxx>
  50 #include <Geom_RectangularTrimmedSurface.hxx>
  51 #include <Geom_SphericalSurface.hxx>
  52 #include <Geom_Surface.hxx>
  53 #include <Geom_SurfaceOfLinearExtrusion.hxx>
  54 #include <Geom_TrimmedCurve.hxx>
  55 #include <GeomAdaptor_HCurve.hxx>
  56 #include <GeomAdaptor_HSurface.hxx>
  57 #include <GeomAPI_Interpolate.hxx>
  58 #include <GeomAPI_PointsToBSpline.hxx>
  59 #include <GeomProjLib.hxx>
  60 #include <gp_Pnt2d.hxx>
  61 #include <NCollection_Sequence.hxx>
  62 #include <Precision.hxx>
  63 #include <ProjLib_CompProjectedCurve.hxx>
  64 #include <ProjLib_HCompProjectedCurve.hxx>
  65 #include <ProjLib_ProjectedCurve.hxx>
  66 #include <ShapeAnalysis_Curve.hxx>
  67 #include <ShapeAnalysis_Surface.hxx>
  68 #include <ShapeConstruct_ProjectCurveOnSurface.hxx>
  69 #include <ShapeExtend.hxx>
  70 #include <Standard_ErrorHandler.hxx>
  71 #include <Standard_Failure.hxx>
  72 #include <Standard_Type.hxx>
  73 #include <TColgp_Array1OfPnt.hxx>
  74 #include <TColStd_Array1OfInteger.hxx>
  75
  76 #include <algorithm>
  77 #define NCONTROL 23
  78
  79 //=======================================================================
  80 //function : ShapeConstruct_ProjectCurveOnSurface
  81 //purpose  :
  82 //=======================================================================
  83
  84 ShapeConstruct_ProjectCurveOnSurface::ShapeConstruct_ProjectCurveOnSurface()
  85 {
  86   myPreci = Precision::Confusion();
  87   myBuild = Standard_False;
  88   myAdjustOverDegen = 1; //:c0 //szv#4:S4163:12Mar99 was boolean
  89   myNbCashe = 0; //:q9
  90 }
  91
  92 //=======================================================================
  93 //function : Init
  94 //purpose  :
  95 //=======================================================================
  96
  97  void ShapeConstruct_ProjectCurveOnSurface::Init(const Handle(Geom_Surface)& surf,const Standard_Real preci)
  98 {
  99   Init (new ShapeAnalysis_Surface (surf), preci);
 100 }
 101
 102 //=======================================================================
 103 //function : Init
 104 //purpose  :
 105 //=======================================================================
 106
 107  void ShapeConstruct_ProjectCurveOnSurface::Init(const Handle(ShapeAnalysis_Surface)& surf,const Standard_Real preci)
 108 {
 109   SetSurface (surf);
 110   SetPrecision (preci);
 111 }
 112
 113 //=======================================================================
 114 //function : SetSurface
 115 //purpose  :
 116 //=======================================================================
 117
 118  void ShapeConstruct_ProjectCurveOnSurface::SetSurface(const Handle(Geom_Surface)& surf)
 119 {
 120   SetSurface (new ShapeAnalysis_Surface (surf));
 121 }
 122
 123 //=======================================================================
 124 //function : SetSurface
 125 //purpose  :
 126 //=======================================================================
 127
 128  void ShapeConstruct_ProjectCurveOnSurface::SetSurface(const Handle(ShapeAnalysis_Surface)& surf)
 129 {
 130   if ( mySurf == surf ) return;
 131   mySurf = surf;
 132   myNbCashe = 0; //:q9
 133 }
 134
 135 //=======================================================================
 136 //function : SetPrecision
 137 //purpose  :
 138 //=======================================================================
 139
 140  void ShapeConstruct_ProjectCurveOnSurface::SetPrecision(const Standard_Real preci)
 141 {
 142   myPreci = preci;
 143 }
 144
 145 //=======================================================================
 146 //function : BuildCurveMode
 147 //purpose  :
 148 //=======================================================================
 149
 150  Standard_Boolean& ShapeConstruct_ProjectCurveOnSurface::BuildCurveMode()
 151 {
 152   return myBuild;
 153 }
 154
 155 //=======================================================================
 156 //function : AdjustOverDegenMode
 157 //purpose  :
 158 //=======================================================================
 159 //:c0
 160
 161 //szv#4:S4163:12Mar99 was Boolean
 162  Standard_Integer& ShapeConstruct_ProjectCurveOnSurface::AdjustOverDegenMode()
 163  {
 164    return myAdjustOverDegen;
 165  }
 166
 167
 168 //=======================================================================
 169 //function : NbSurfIntervals
 170 //purpose  : work-around of bug in standard method
 171 //           GeomAdaptor_Surface->NbIntervals() (PRO16346)
 172 //=======================================================================
 173
 174 static Standard_Integer NbSurfIntervals(const Handle(GeomAdaptor_HSurface)& GAS, const GeomAbs_Shape cont)
 175 {
 176   Standard_Integer NbU = 0;
 177   if (GAS->GetType() == GeomAbs_SurfaceOfExtrusion) {
 178     // extract the surface
 179     Handle(Geom_SurfaceOfLinearExtrusion) surf = Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(GAS->ChangeSurface().Surface());
 180     // build a 3d adaptor curve
 181     GeomAdaptor_Curve Adaptor3dCurve(surf->BasisCurve(), GAS->FirstUParameter(), GAS->LastUParameter());
 182     if (Adaptor3dCurve.GetType() == GeomAbs_BSplineCurve)
 183       NbU = Adaptor3dCurve.NbIntervals(cont);
 184   }
 185   if (NbU == 0)
 186     NbU = GAS->NbUIntervals(cont);
 187   return NbU * (GAS->NbVIntervals(cont));
 188 }
 189
 190 //=======================================================================
 191 //function : Status
 192 //purpose  :
 193 //=======================================================================
 194
 195  Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::Status (const ShapeExtend_Status Status) const
 196 {
 197   return ShapeExtend::DecodeStatus (myStatus, Status);
 198 }
 199
 200 //=======================================================================
 201 //function : Perform
 202 //purpose  :
 203 //=======================================================================
 204 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::Perform (Handle(Geom_Curve)& c3d,
 205                                                                 const Standard_Real First,
 206                                                                 const Standard_Real Last,
 207                                                                 Handle(Geom2d_Curve)& c2d,
 208                                                                 const GeomAbs_Shape,
 209                                                                 const Standard_Integer,
 210                                                                 const Standard_Integer)
 211 {
 212   myStatus = ShapeExtend::EncodeStatus (ShapeExtend_OK);
 213   //Standard_Boolean OK = Standard_True; //szv#4:S4163:12Mar99 not needed
 214
 215   if (mySurf.IsNull()) {
 216     c2d.Nullify();
 217     myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_FAIL1);
 218     return Standard_False;
 219   }
 220 //  Projection Analytique
 221   Handle(Geom_Curve) crv3dtrim = c3d;
 222   if ( ! c3d->IsKind(STANDARD_TYPE(Geom_BoundedCurve)) )
 223     crv3dtrim = new Geom_TrimmedCurve ( c3d, First, Last );
 224   c2d = ProjectAnalytic ( crv3dtrim );
 225   if (!c2d.IsNull()) {
 226     myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE1);
 227     return Standard_True;
 228   }
 229
 230 //  Projection par approximation
 231
 232   // discretize the 3d curve
 233
 234   Standard_Integer nbrPnt;
 235
 236 //   $$$$    :92 abv 28 Jan 98   see PRO10107, big BSplineCurve C0
 237   Standard_Integer nbPini = NCONTROL;  // as in BRepCheck_Edge (RLN/Nijni)
 238  // 20; // number of points for interpolation, should be "parametric dependent"
 239
 240   //:92 abv 28 Jan 98: if curve is BSpline with many intervals,
 241   // increase number of points to provide at least Degree()+1 points per interval
 242   Handle(Geom_BSplineCurve) bspl;
 243   if ( c3d->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)) ) {
 244     Handle(Geom_TrimmedCurve) ctrim = Handle(Geom_TrimmedCurve)::DownCast(c3d);
 245     bspl = Handle(Geom_BSplineCurve)::DownCast ( ctrim->BasisCurve() );
 246   }
 247   else bspl = Handle(Geom_BSplineCurve)::DownCast ( c3d );
 248   if ( ! bspl.IsNull() ) {
 249     Standard_Integer nint = 0;
 250     for ( Standard_Integer i=1; i < bspl->NbKnots(); i++ )
 251       if ( bspl->Knot(i+1) > First && bspl->Knot(i) < Last ) nint++;
 252     Standard_Integer minPnt = nint * ( bspl->Degree() + 1 );
 253     while ( nbPini < minPnt ) nbPini += NCONTROL - 1;
 254 #ifdef OCCT_DEBUG
 255     if ( nbPini > NCONTROL )
 256       cout << "Warning: number of points for projecting is " << nbPini << endl;
 257 #endif
 258   }
 259
 260 //    $$$$    end :92 (big BSplineCurve C0)
 261
 262   // this number should be "parametric dependent"
 263   TColgp_Array1OfPnt points(1, nbPini);
 264   TColStd_Array1OfReal params(1, nbPini);
 265   NCollection_Sequence<Standard_Real> aKnotCoeffs;
 266   gp_Pnt p3d;
 267   Standard_Integer iPnt;
 268
 269   // In case of bspline compute parametrization speed on each
 270   // knot interval inside [aFirstParam, aLastParam].
 271   // If quotient = (MaxSpeed / MinSpeed) >= aMaxQuotientCoeff then
 272   // use PerformByProjLib algorithm.
 273   if(!bspl.IsNull())
 274   {
 275     Standard_Integer aNbIntPnts = NCONTROL;
 276     Standard_Real aFirstParam = First; // First parameter of current interval.
 277     Standard_Real aLastParam = Last; // Last parameter of current interval.
 278
 279     // First index computation.
 280     Standard_Integer anIdx = 1;
 281     for(; anIdx <= bspl->NbKnots() && aFirstParam < Last; anIdx++)
 282     {
 283       if(bspl->Knot(anIdx) > First)
 284       {
 285         break;
 286       }
 287     }
 288
 289     for(; anIdx <= bspl->NbKnots() && aFirstParam < Last; anIdx++)
 290     {
 291       // Fill current knot interval.
 292       aLastParam = Min(Last, bspl->Knot(anIdx));
 293       Standard_Real aStep = (aLastParam - aFirstParam) / (aNbIntPnts - 1);
 294       Standard_Integer anIntIdx;
 295       gp_Pnt p3d1, p3d2;
 296       Standard_Real aLength3d = 0.0;
 297       for(anIntIdx = 0; anIntIdx < aNbIntPnts - 1; anIntIdx++)
 298       {
 299         // Start filling from first point.
 300         Standard_Real aParam1 = aFirstParam + aStep * anIntIdx;
 301         Standard_Real aParam2 = aFirstParam + aStep * (anIntIdx + 1);
 302         c3d->D0 (aParam1, p3d1);
 303         c3d->D0 (aParam2, p3d2);
 304         aLength3d += p3d2.Distance(p3d1);
 305       }
 306       aKnotCoeffs.Append(aLength3d / (aLastParam - aFirstParam));
 307       aFirstParam = aLastParam;
 308     }
 309
 310     Standard_Real anEvenlyCoeff = 0;
 311     if (aKnotCoeffs.Size() > 0)
 312     {
 313       anEvenlyCoeff = *std::max_element(aKnotCoeffs.begin(), aKnotCoeffs.end()) /
 314                       *std::min_element(aKnotCoeffs.begin(), aKnotCoeffs.end());
 315     }
 316
 317     const Standard_Real aMaxQuotientCoeff = 1500.0;
 318     if (anEvenlyCoeff > aMaxQuotientCoeff)
 319     {
 320       PerformByProjLib(c3d, First, Last, c2d);
 321       // PerformByProjLib fail detection:
 322       if (!c2d.IsNull())
 323       {
 324         return Status (ShapeExtend_DONE);
 325       }
 326     }
 327   }
 328
 329   Standard_Real deltaT, t;
 330   deltaT = (Last - First) / (nbPini-1);
 331   nbrPnt = nbPini;
 332   for (iPnt = 1; iPnt <= nbPini; iPnt ++)
 333   {
 334     if      (iPnt == 1)      t = First;
 335     else if (iPnt == nbPini) t = Last;
 336     else                     t = First + (iPnt - 1) * deltaT;
 337
 338     c3d->D0 (t, p3d);
 339     points(iPnt) = p3d;
 340     params(iPnt) = t;
 341   }
 342
 343   //  CALCUL par approximation
 344   TColgp_Array1OfPnt2d pnt2d(1, nbrPnt);
 345   ApproxPCurve (nbrPnt,points,params,pnt2d,c2d); //szv#4:S4163:12Mar99 OK not needed
 346   if (!c2d.IsNull()) {
 347     myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE2);
 348     return Standard_True;
 349   }// cas particulier d iso
 350
 351 //  INTERPOLATION du resultat
 352
 353   if ( myBuild ) {
 354     Handle(TColgp_HArray1OfPnt) thePnts = new TColgp_HArray1OfPnt  (1, nbPini);
 355     Handle(TColStd_HArray1OfReal) theParams = new TColStd_HArray1OfReal(1, nbPini);
 356     for (iPnt = 1; iPnt <= nbPini ; iPnt ++) {
 357       thePnts->SetValue(iPnt, points(iPnt));
 358       theParams->SetValue(iPnt, params(iPnt));
 359     }
 360
 361     Handle(Geom_Curve) newc3d = InterpolateCurve3d (nbPini,thePnts,theParams, c3d);
 362     if ( newc3d.IsNull() ) myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_FAIL2);
 363     else {
 364       myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE3);
 365       c3d = newc3d;
 366     }
 367   }
 368
 369   Handle(TColgp_HArray1OfPnt2d) thePnts2d = new TColgp_HArray1OfPnt2d(1, nbPini);
 370   Handle(TColStd_HArray1OfReal) theParams2d = new TColStd_HArray1OfReal(1, nbPini);
 371   for (iPnt = 1; iPnt <= nbPini ; iPnt ++) {
 372     theParams2d->SetValue(iPnt, params(iPnt));
 373     thePnts2d->SetValue(iPnt, pnt2d(iPnt));
 374   }
 375
 376   c2d = InterpolatePCurve (nbPini, thePnts2d, theParams2d, c3d);
 377 //  c2d = ApproximatePCurve (nbPini, thePnts2d, theParams2d, c3d);
 378 //  Faut-il aussi reprendre la C3D ?
 379   myStatus |= ShapeExtend::EncodeStatus (c2d.IsNull() ? ShapeExtend_FAIL1 : ShapeExtend_DONE2);
 380   return Status (ShapeExtend_DONE);
 381 }
 382
 383 //=======================================================================
 384 //function : PerformByProjLib
 385 //purpose  :
 386 //=======================================================================
 387
 388 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformByProjLib(Handle(Geom_Curve)& c3d,
 389                                                                         const Standard_Real First,
 390                                                                         const Standard_Real Last,
 391                                                                         Handle(Geom2d_Curve)& c2d,
 392                                                                         const GeomAbs_Shape /*continuity*/,
 393                                                                         const Standard_Integer /*maxdeg */,
 394                                                                         const Standard_Integer /*nbinterval */)
 395 {
 396   //Standard_Boolean OK = Standard_True; //szv#4:S4163:12Mar99 unused
 397   c2d.Nullify();
 398   if (mySurf.IsNull()) {
 399     myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL1);
 400     return Standard_False;
 401   }
 402
 403   try
 404   {
 405     OCC_CATCH_SIGNALS
 406     Handle(GeomAdaptor_HSurface) GAS = mySurf->Adaptor3d();
 407     Handle(GeomAdaptor_HCurve) GAC = new GeomAdaptor_HCurve (c3d,First,Last);
 408     ProjLib_ProjectedCurve Projector(GAS, GAC);
 409
 410     switch (Projector.GetType())
 411     {
 412     case GeomAbs_Line :
 413       c2d = new Geom2d_Line(Projector.Line());
 414       break;
 415     case GeomAbs_Circle :
 416       c2d = new Geom2d_Circle(Projector.Circle());
 417       break;
 418     case GeomAbs_Ellipse :
 419       c2d = new Geom2d_Ellipse(Projector.Ellipse());
 420       break;
 421     case GeomAbs_Parabola :
 422       c2d = new Geom2d_Parabola(Projector.Parabola());
 423       break;
 424     case GeomAbs_Hyperbola :
 425       c2d = new Geom2d_Hyperbola(Projector.Hyperbola());
 426       break;
 427     case GeomAbs_BSplineCurve :
 428       c2d = Projector.BSpline();
 429       break;
 430     case GeomAbs_BezierCurve :
 431     case GeomAbs_OtherCurve :
 432       // Not possible, handling added to avoid gcc warning.
 433       break;
 434     }
 435
 436     if(c2d.IsNull())
 437     {
 438       myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL2);
 439       return Standard_False;
 440     }
 441     else
 442     {
 443       myStatus = ShapeExtend::EncodeStatus (ShapeExtend_DONE1);
 444       return Standard_True;
 445     }
 446
 447   }
 448   catch(Standard_Failure)
 449   {
 450 #ifdef OCCT_DEBUG
 451     cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::PerformByProjLib(): Exception: ";
 452     Standard_Failure::Caught()->Print(cout); cout << endl;
 453 #endif
 454     myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL3);
 455     c2d.Nullify();
 456   }
 457   return Standard_False;
 458 }
 459
 460 //=======================================================================
 461 //function : PerformAdvanced
 462 //purpose  :
 463 //=======================================================================
 464
 465 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(Geom_Curve)& c3d,
 466                                                                         const Standard_Real First,
 467                                                                         const Standard_Real Last,
 468                                                                         Handle(Geom2d_Curve)& c2d)
 469 {
 470   Standard_Boolean hasResult = Standard_False;
 471   Standard_Integer nbintervals;
 472
 473   Standard_Boolean isStandard = (mySurf->Adaptor3d()->GetType() != GeomAbs_Cylinder);
 474 // &&   (mySurf->Adaptor3d()->GetType() != GeomAbs_SurfaceOfRevolution);
 475
 476   if (isStandard) isStandard = !mySurf->HasSingularities(myPreci);
 477   if (isStandard) {
 478     Handle(GeomAdaptor_HSurface) GAS = mySurf->Adaptor3d();
 479     Handle(GeomAdaptor_HCurve) GAC = new GeomAdaptor_HCurve (c3d,First,Last);
 480     nbintervals = NbSurfIntervals(GAS, GeomAbs_C1);//+GAC->NbIntervals(GeomAbs_C3);
 481     isStandard = (nbintervals < 2);
 482   }
 483   if (isStandard) {
 484     hasResult = PerformByProjLib(c3d, First, Last, c2d);
 485   }
 486   if (!hasResult) hasResult = Perform (c3d, First, Last, c2d);
 487   return hasResult;
 488 }
 489
 490 //=======================================================================
 491 //function : ProjectAnalytic
 492 //purpose  :
 493 //=======================================================================
 494
 495  Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::ProjectAnalytic(const Handle(Geom_Curve)& c3d) const
 496 {
 497   Handle(Geom2d_Curve) result;
 498
 499   //:k1 abv 16 Dec 98: limit analytic cases by Plane surfaces only
 500   // This is necessary for K4L since it fails on other surfaces
 501   // when general method GeomProjLib::Curve2d() is used
 502   // Projection is done as in BRep_Tool and BRepCheck_Edge
 503   Handle(Geom_Surface) surf = mySurf->Surface();
 504   Handle(Geom_Plane) Plane = Handle(Geom_Plane)::DownCast ( surf );
 505   if ( Plane.IsNull() ) {
 506     Handle(Geom_RectangularTrimmedSurface) RTS =
 507       Handle(Geom_RectangularTrimmedSurface)::DownCast ( surf );
 508     if ( ! RTS.IsNull() ) Plane = Handle(Geom_Plane)::DownCast ( RTS->BasisSurface() );
 509     else {
 510       Handle(Geom_OffsetSurface) OS =
 511         Handle(Geom_OffsetSurface)::DownCast ( surf );
 512       if ( ! OS.IsNull() )
 513         Plane = Handle(Geom_Plane)::DownCast ( OS->BasisSurface() );
 514     }
 515   }
 516   if ( ! Plane.IsNull() ) {
 517     Handle(Geom_Curve) ProjOnPlane =
 518       GeomProjLib::ProjectOnPlane (c3d, Plane,
 519                                    Plane->Position().Direction(), Standard_True);
 520     Handle(GeomAdaptor_HCurve) HC = new GeomAdaptor_HCurve ( ProjOnPlane );
 521     ProjLib_ProjectedCurve Proj ( mySurf->Adaptor3d(), HC );
 522
 523     result = Geom2dAdaptor::MakeCurve(Proj);
 524     if ( result.IsNull() ) return result;
 525     if ( result->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve)) ) {
 526       Handle(Geom2d_TrimmedCurve) TC = Handle(Geom2d_TrimmedCurve)::DownCast ( result );
 527       result = TC->BasisCurve();
 528     }
 529
 530     return result;
 531   }
 532
 533   return result;
 534 }
 535
 536  //! Fix possible period jump and handle walking period parameter.
 537  static Standard_Boolean fixPeriodictyTroubles(gp_Pnt2d *thePnt, // pointer to gp_Pnt2d[4] beginning
 538                                                Standard_Integer theIdx, // Index of objective coord: 1 ~ X, 2 ~ Y
 539                                                Standard_Real thePeriod) // Period on objective coord
 540  {
 541    Standard_Integer i;
 542
 543    Standard_Boolean isNeedToFix = Standard_True;
 544    for (i = 0; i < 3; i++)
 545    {
 546      Standard_Real aDiff = Abs (thePnt[i].Coord(theIdx) - thePnt[i + 1].Coord(theIdx));
 547      if ( aDiff > Precision::PConfusion() &&
 548           aDiff < thePeriod - Precision::PConfusion())
 549      {
 550        // Walk over period coord -> not walking on another isoline in parameter space.
 551        isNeedToFix = Standard_False;
 552      }
 553    }
 554
 555    if (isNeedToFix)
 556    {
 557      // Walking on isoline on another parameter. Fix period paramter to obtained minimum.
 558      Standard_Real aFixParam = Min (thePnt[0].Coord(theIdx), thePnt[3].Coord(theIdx));
 559      for(i = 0; i < 4; i++)
 560        thePnt[i].SetCoord(theIdx, aFixParam);
 561    }
 562
 563    // Fix possible period jump on first point.
 564    if ( Abs(thePnt[0].Coord(theIdx) - thePnt[1].Coord(theIdx) ) >  thePeriod / 2.01)
 565    {
 566      Standard_Real aMult = thePnt[0].Coord(theIdx) < thePnt[1].Coord(theIdx) ? 1.0 : -1.0;
 567      Standard_Real aNewParam = thePnt[0].Coord(theIdx) + aMult * thePeriod;
 568      thePnt[0].SetCoord(theIdx, aNewParam);
 569      return Standard_False;
 570    }
 571
 572    // Fix possible period jump on last point.
 573    if ( Abs(thePnt[2].Coord(theIdx) - thePnt[3].Coord(theIdx) ) >  thePeriod / 2.01)
 574    {
 575      Standard_Real aMult = thePnt[3].Coord(theIdx) < thePnt[2].Coord(theIdx) ? 1.0 : -1.0;
 576      Standard_Real aNewParam = thePnt[3].Coord(theIdx) + aMult * thePeriod;
 577      thePnt[3].SetCoord(theIdx, aNewParam);
 578      return Standard_False;
 579    }
 580
 581    return Standard_True;
 582  }
 583
 584 //=======================================================================
 585 //function : getLine
 586 //purpose  :
 587 //=======================================================================
 588
 589  Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::getLine(
 590    const TColgp_Array1OfPnt& thepoints,
 591    const TColStd_Array1OfReal& theparams,
 592    TColgp_Array1OfPnt2d& thePnt2ds,
 593    Standard_Real theTol,
 594    Standard_Boolean &isRecompute) const
 595  {
 596    Standard_Integer nb =  thepoints.Length();
 597    gp_Pnt aP[4];
 598    aP[0] = thepoints(1);
 599    aP[1] = thepoints(2);
 600    aP[2] = thepoints(nb - 1);
 601    aP[3] = thepoints(nb);
 602    gp_Pnt2d aP2d[4];
 603    Standard_Integer i = 0;
 604
 605    Standard_Real aTol2 = theTol * theTol;
 606    Standard_Boolean isPeriodicU = mySurf->Surface()->IsUPeriodic();
 607    Standard_Boolean isPeriodicV = mySurf->Surface()->IsVPeriodic();
 608
 609    // Workaround:
 610    // Protection against bad "tolerance" shapes.
 611    if (aTol2 > 1.0)
 612    {
 613      theTol = Precision::Confusion();
 614      aTol2 = theTol * theTol;
 615    }
 616    Standard_Real anOldTol2 = aTol2;
 617
 618    // project first and last points
 619    for( ; i < 4; i +=3)
 620    {
 621      Standard_Integer j;
 622      for ( j=0; j < myNbCashe; j++ )
 623        if ( myCashe3d[j].SquareDistance (aP[i] ) < aTol2)
 624        {
 625          aP2d[i] = mySurf->NextValueOfUV (myCashe2d[j], aP[i], theTol,
 626            theTol);
 627          break;
 628        }
 629        if ( j >= myNbCashe )
 630          aP2d[i] = mySurf->ValueOfUV(aP[i], theTol);
 631
 632        Standard_Real aDist = mySurf->Gap();
 633        Standard_Real aCurDist = aDist * aDist;
 634        if( aTol2 < aDist * aDist)
 635          aTol2 = aCurDist;
 636    }
 637
 638    if ( isPeriodicU || isPeriodicV )
 639    {
 640      // Compute second and last but one c2d points.
 641      for(i = 1; i < 3; i++)
 642      {
 643        Standard_Integer j;
 644        for ( j=0; j < myNbCashe; j++ )
 645          if ( myCashe3d[j].SquareDistance (aP[i] ) < aTol2)
 646          {
 647            aP2d[i] = mySurf->NextValueOfUV (myCashe2d[j], aP[i], theTol, theTol);
 648            break;
 649          }
 650          if ( j >= myNbCashe )
 651            aP2d[i] = mySurf->ValueOfUV(aP[i], theTol);
 652
 653          Standard_Real aDist = mySurf->Gap();
 654          Standard_Real aCurDist = aDist * aDist;
 655          if( aTol2 < aDist * aDist)
 656            aTol2 = aCurDist;
 657      }
 658
 659      if (isPeriodicU)
 660      {
 661        isRecompute = fixPeriodictyTroubles(&aP2d[0], 1 /* X Coord */, mySurf->Surface()->UPeriod());
 662      }
 663
 664      if (isPeriodicV)
 665      {
 666        isRecompute = fixPeriodictyTroubles(&aP2d[0], 2 /* Y Coord */, mySurf->Surface()->VPeriod());
 667      }
 668    }
 669
 670    thePnt2ds.SetValue(1, aP2d[0]);
 671    thePnt2ds.SetValue(nb, aP2d[3]);
 672
 673    // Restore old tolerance in 2d space to avoid big gap cases.
 674    aTol2 = anOldTol2;
 675    // Check that straight line in 2d with parameterisation as in 3d will fit
 676    // fit 3d curve at all points.
 677    Standard_Real dPar = theparams(nb) - theparams(1);
 678    if ( Abs(dPar) < Precision::PConfusion() )
 679      return 0;
 680    gp_Vec2d aVec0 (aP2d[0], aP2d[3]);
 681    gp_Vec2d aVec = aVec0 / dPar;
 682    Standard_Real aFirstPointDist = mySurf->Surface()->Value(aP2d[0].X(), aP2d[0].Y()).
 683                                    SquareDistance(thepoints(1));
 684    for(i = 2; i <  nb; i++)
 685    {
 686      gp_XY aCurPoint = aP2d[0].XY() + aVec.XY() * (theparams(i) - theparams(1));
 687      gp_Pnt aCurP;
 688      mySurf->Surface()->D0(aCurPoint.X(), aCurPoint.Y(), aCurP);
 689      Standard_Real aDist1 = aCurP.SquareDistance(thepoints(i));
 690
 691      if(Abs (aFirstPointDist - aDist1) > aTol2)
 692        return 0;
 693    }
 694
 695    // check if pcurve can be represented by Geom2d_Line (parameterised by length)
 696    Standard_Real aLLength = aVec0.Magnitude();
 697    if ( Abs (aLLength - dPar) <= Precision::PConfusion() )
 698    {
 699      gp_XY aDirL = aVec0.XY() / aLLength;
 700      gp_Pnt2d aPL (aP2d[0].XY() - theparams(1) * aDirL);
 701      return new Geom2d_Line (aPL, gp_Dir2d(aDirL));
 702    }
 703
 704    // create straight bspline
 705    TColgp_Array1OfPnt2d aPoles(1, 2);
 706    aPoles(1) = aP2d[0];
 707    aPoles(2) = aP2d[3];
 708
 709    TColStd_Array1OfReal aKnots(1,2);
 710    aKnots(1) = theparams(1);
 711    aKnots(2) = theparams(theparams.Length());
 712
 713    TColStd_Array1OfInteger aMults(1,2);
 714    aMults(1) = 2;
 715    aMults(2) = 2;
 716    Standard_Integer aDegree = 1;
 717    Handle(Geom2d_BSplineCurve) abspl2d =
 718      new Geom2d_BSplineCurve (aPoles, aKnots, aMults, aDegree);
 719    return abspl2d;
 720  }
 721
 722 //=======================================================================
 723 //function : ApproxPCurve
 724 //purpose  :
 725 //=======================================================================
 726
 727   Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::ApproxPCurve(const Standard_Integer nbrPnt,
 728                                                                      const TColgp_Array1OfPnt& points,
 729                                                                      const TColStd_Array1OfReal& params,
 730                                                                      TColgp_Array1OfPnt2d& pnt2d,
 731                                                                      Handle(Geom2d_Curve)& c2d)
 732 {
 733   // for performance, first try to handle typical case when pcurve is straight
 734   Standard_Boolean isRecompute = Standard_False;
 735   c2d = getLine(points, params, pnt2d, myPreci, isRecompute);
 736   if(!c2d.IsNull())
 737   {
 738     return Standard_True;
 739   }
 740     Standard_Boolean isDone = Standard_True;
 741   // test if the curve 3d is a boundary of the surface
 742   // (only for Bezier or BSpline surface)
 743
 744   Standard_Boolean isoParam, isoPar2d3d, isoTypeU, p1OnIso, p2OnIso, isoclosed;
 745   gp_Pnt2d valueP1, valueP2;
 746   Handle(Geom_Curve) cIso;
 747   Standard_Real t1, t2;
 748
 749   Handle(Standard_Type) sType = mySurf->Surface()->DynamicType();
 750   Standard_Boolean isAnalytic = Standard_True;
 751   if (sType == STANDARD_TYPE(Geom_BezierSurface) || sType == STANDARD_TYPE(Geom_BSplineSurface)) isAnalytic = Standard_False;
 752   Standard_Real uf, ul, vf, vl;
 753   mySurf->Surface()->Bounds(uf, ul, vf, vl);
 754   isoclosed = Standard_False;
 755   TColStd_Array1OfReal pout(1, nbrPnt);
 756
 757   isoParam = IsAnIsoparametric(nbrPnt, points, params,
 758                                isoTypeU, p1OnIso, valueP1, p2OnIso, valueP2,
 759                                isoPar2d3d, cIso, t1, t2, pout);
 760
 761   // projection of the points on surfaces
 762
 763   gp_Pnt p3d;
 764   gp_Pnt2d p2d;
 765   Standard_Integer i;
 766   Standard_Real isoValue=0., isoPar1=0., isoPar2=0., tPar=0., tdeb,tfin;
 767   Standard_Real Cf, Cl, parf, parl; //szv#4:S4163:12Mar99 dist not needed
 768
 769   //  Le calcul part-il dans le bon sens, c-a-d deb et fin dans le bon ordre ?
 770   //  Si uclosed et iso en V, attention isoPar1 ET/OU 2 peut toucher la fermeture
 771   if(isoParam){
 772     if(isoTypeU){
 773       isoValue = valueP1.X();
 774       isoPar1 = valueP1.Y();
 775       isoPar2 = valueP2.Y();
 776       isoclosed = mySurf->IsVClosed(myPreci);//#78 rln 12.03.99 S4135
 777       parf = vf;  parl = vl;
 778     }
 779     else {
 780       isoValue = valueP1.Y();
 781       isoPar1 = valueP1.X();
 782       isoPar2 = valueP2.X();
 783       isoclosed = mySurf->IsUClosed(myPreci);//#78 rln 12.03.99 S4135
 784       parf = uf;  parl = ul;
 785     }
 786     if (!isoPar2d3d && !isAnalytic) {
 787       Cf = cIso->FirstParameter();
 788       Cl = cIso->LastParameter();
 789       if (Precision::IsInfinite(Cf))    Cf = -1000;
 790       if (Precision::IsInfinite(Cl))    Cl = +1000;
 791       //pdn S4030 optimizing and fix isopar case on PRO41323
 792       tdeb = pout(2);
 793       //    dist = ShapeAnalysis_Curve().Project (cIso,points(2),myPreci,pt,tdeb,Cf,Cl);
 794       //  Chacun des par1 ou par2 est-il sur un bord. Attention first/last : recaler
 795       if (isoclosed && (isoPar1 == parf || isoPar1 == parl)) {
 796         if (Abs(tdeb-parf) < Abs(tdeb-parl)) isoPar1 = parf;
 797         else isoPar1 = parl;
 798         if (isoTypeU) valueP1.SetY (isoPar1);
 799         else          valueP1.SetX (isoPar1);
 800       }
 801       if (isoclosed && (isoPar2 == parf || isoPar2 == parl)) {
 802         //pdn S4030 optimizing and fix isopar case on PRO41323
 803         tfin = pout(nbrPnt-1);
 804         //dist =  ShapeAnalysis_Curve().Project (cIso,points(nbrPnt-1),myPreci,pt,tfin,Cf,Cl);
 805         if (Abs(tfin-parf) < Abs(tfin-parl)) isoPar2 = parf;
 806         else isoPar2 = parl;
 807         if (isoTypeU) valueP2.SetY (isoPar2);
 808         else          valueP2.SetX (isoPar2);
 809       }
 810
 811       //  Interversion Par1/Par2 (ne veut que si les 2 sont sur les bords ...)
 812       //  Est-ce encore necessaire apres ce qui vient d etre fait ?
 813
 814       // PTV 05.02.02 fix for translation face from 12_hp_mouse (PARASOLID) face 24008
 815       // if curve is periodic do not change the points
 816       // skl change "if" for pout(nbrPnt-1) 19.11.2003
 817       if (!isoclosed) {
 818         if( (Abs(tdeb-isoPar1)>Abs(tdeb-isoPar2)) &&
 819             (Abs(pout(nbrPnt-1)-isoPar2)>Abs(pout(nbrPnt-1)-isoPar1)) ) {
 820           gp_Pnt2d valueTmp = valueP1;
 821           valueP1 = valueP2;  valueP2 = valueTmp;
 822           if (isoTypeU) {
 823             isoValue = valueP1.X();
 824             isoPar1 = valueP1.Y();
 825             isoPar2 = valueP2.Y();
 826           }
 827           else {
 828             isoValue = valueP1.Y();
 829             isoPar1 = valueP1.X();
 830             isoPar2 = valueP2.X();
 831           }
 832           //  Fin calcul sens de courbe iso
 833         }
 834       } // end of fix check 05.02.02
 835     }
 836   }
 837
 838   //  Si pas isoParam, on a quand meme du p1OnIso/p2OnIso possible ... !!!
 839   //  (utile pour detromper bug de projection). Mais detromper aussi circularite
 840   //else {
 841     //if (p1OnIso) valueP1 =
 842     //BestExtremum (valueP1,points(1),points(2));
 843     //if (p2OnIso) valueP2 =
 844     //BestExtremum (valueP2,points(nbrPnt),points(nbrPnt-1));
 845     //}
 846
 847   Standard_Real gap = myPreci; //:q1
 848   Standard_Boolean ChangeCycle = Standard_False; //skl for OCC3430
 849   if( myNbCashe>0 && myCashe3d[0].Distance(points(1))>myCashe3d[0].Distance(points(nbrPnt)) )
 850     //if(myCashe3d[0].Distance(points(nbrPnt))<myPreci)
 851     if(myCashe3d[0].Distance(points(nbrPnt))<Precision::Confusion())
 852       ChangeCycle = Standard_True;
 853   //for( i = 1; i <= nbrPnt; i ++) {
 854   for(Standard_Integer ii=1; ii<=nbrPnt; ii++) {
 855     if(ChangeCycle) //skl for OCC3430
 856       i=nbrPnt-ii+1;
 857     else
 858       i=ii;
 859     p3d = points(i);
 860     if (isoParam) {
 861
 862       if (isoPar2d3d) {
 863         if (isoPar2 > isoPar1) tPar = params(i);
 864         else                   tPar = t1 + t2 - params(i);
 865       } else if (!isAnalytic) {
 866         // projection to iso
 867         if      (i==1)      tPar = isoPar1;
 868         else if (i==nbrPnt) tPar = isoPar2;
 869         else {
 870           tPar = pout(i);
 871           //:S4030  ShapeAnalysis_Curve().Project (cIso,p3d,myPreci,pt,tPar,Cf,Cl); //szv#4:S4163:12Mar99 `dist=` not needed
 872         }
 873       }
 874
 875       if (!isoPar2d3d && isAnalytic) {
 876         if      (i == 1)      p2d = valueP1;
 877         else if (i == nbrPnt) p2d = valueP2;
 878         else {
 879           p2d = mySurf->NextValueOfUV(p2d,p3d, myPreci, //%12 pdn 15.02.99 optimizing
 880                                       Precision::Confusion()+1000*gap); //:q1
 881           gap = mySurf->Gap();
 882         }
 883       } else {
 884         if(isoTypeU)  {  p2d.SetX(isoValue);  p2d.SetY(tPar);     }
 885         else          {  p2d.SetX(tPar);      p2d.SetY(isoValue); }
 886       }
 887     }
 888
 889     else {
 890       if     ( (i == 1)      && p1OnIso)  p2d = valueP1;
 891       else if( (i == nbrPnt) && p2OnIso)  p2d = valueP2;
 892       else  {// general case (not an iso)  mais attention aux singularites !
 893         // first and last points are already computed by getLine()
 894         if ( (i == 1 || i == nbrPnt))
 895         {
 896           if (!isRecompute)
 897           {
 898             p2d = pnt2d(i);
 899             gap = mySurf->Gap();
 900             continue;
 901           }
 902           else
 903           {
 904             //:q9 abv 23 Mar 99: use cashe as 1st approach
 905             Standard_Integer j; // svv #1
 906             for ( j=0; j < myNbCashe; j++ )
 907               if ( myCashe3d[j].SquareDistance ( p3d ) < myPreci*myPreci )
 908               {
 909                 p2d = mySurf->NextValueOfUV (myCashe2d[j], p3d, myPreci,
 910                   Precision::Confusion()+gap);
 911                 break;
 912               }
 913               if ( j >= myNbCashe ) p2d = mySurf->ValueOfUV(p3d, myPreci);
 914           }
 915         }
 916         else {
 917           p2d = mySurf->NextValueOfUV (p2d, p3d, myPreci, //:S4030: optimizing
 918                                        Precision::Confusion()+1000*gap); //:q1
 919         }
 920         gap = mySurf->Gap();
 921       }
 922     }
 923     pnt2d (i) = p2d;
 924     if ( ii > 1 ) {
 925       if(ChangeCycle)
 926         p2d.SetXY ( 2. * p2d.XY() - pnt2d(i+1).XY() );
 927       else
 928         p2d.SetXY ( 2. * p2d.XY() - pnt2d(i-1).XY() );
 929     }
 930   }
 931
 932   //pdn %12 11.02.99 PRO9234 entity 15402
 933   if (!isoPar2d3d) {
 934     mySurf->ProjectDegenerated(nbrPnt,points,pnt2d,myPreci,Standard_True);
 935     mySurf->ProjectDegenerated(nbrPnt,points,pnt2d,myPreci,Standard_False);
 936   }
 937
 938   //  attention aux singularites ... (hors cas iso qui les traite deja)
 939   //  if (!isoParam) {
 940     //    p2d = pnt2d (1);
 941     //    if (mySurf->ProjectDegenerated (points(1),myPreci,pnt2d (2),p2d))
 942     //      pnt2d (1) = p2d;
 943     //    p2d = pnt2d (nbrPnt);
 944     //    if (mySurf->ProjectDegenerated (points(nbrPnt),myPreci,pnt2d (nbrPnt-1),p2d))
 945     //      pnt2d (nbrPnt) = p2d;
 946     //  }
 947
 948   // Si la surface est UCLosed et VClosed, on recadre les points
 949   // algo un peu complique, on retarde l implementation
 950   Standard_Real Up = ul - uf;
 951   Standard_Real Vp = vl - vf;
 952   Standard_Real dist2d;
 953 #ifdef OCCT_DEBUG
 954   if (mySurf->IsUClosed(myPreci) && mySurf->IsVClosed(myPreci)) {//#78 rln 12.03.99 S4135
 955     cout << "WARNING : Recadrage incertain sur U & VClosed" << endl;
 956   }
 957 #endif
 958   // Si la surface est UCLosed, on recadre les points
 959   if (mySurf->IsUClosed(myPreci)) {//#78 rln 12.03.99 S4135
 960     // Premier point dans le domain [uf, ul]
 961     Standard_Real prevX, firstX = pnt2d (1).X();
 962     while (firstX < uf)  {  firstX += Up;   pnt2d (1).SetX(firstX);  }
 963     while (firstX > ul)  {  firstX -= Up;   pnt2d (1).SetX(firstX);  }
 964     prevX = firstX;
 965
 966     //:97 abv 1 Feb 98: treat case when curve is whole out of surface bounds
 967     Standard_Real minX = firstX, maxX = firstX;
 968
 969     // On decalle toujours le suivant
 970     for (i = 2; i <= nbrPnt; i++) {
 971       //      dist2d = pnt2d (i-1).Distance(pnt2d (i));
 972       Standard_Real CurX = pnt2d (i).X();
 973       dist2d = Abs (CurX - prevX);
 974       if (dist2d > ( Up / 2) ) {
 975         if        (CurX > prevX + Up/2) {
 976           while (CurX > prevX + Up/2) {  CurX -= Up;  pnt2d (i).SetX (CurX);  }
 977         } else if (CurX < prevX - Up/2) {
 978           while (CurX < prevX - Up/2) {  CurX += Up;  pnt2d (i).SetX (CurX);  }
 979         }
 980
 981       }
 982       prevX = CurX;
 983       if ( minX > CurX ) minX = CurX;      //:97
 984       else if ( maxX < CurX ) maxX = CurX; //:97
 985     }
 986
 987     //:97
 988     Standard_Real midX = 0.5 * ( minX + maxX );
 989     Standard_Real shiftX=0.;
 990     if ( midX > ul ) shiftX = -Up;
 991     else if ( midX < uf ) shiftX = Up;
 992     if ( shiftX != 0. )
 993       for ( i=1; i <= nbrPnt; i++ ) pnt2d(i).SetX ( pnt2d(i).X() + shiftX );
 994   }
 995   // Si la surface est VCLosed, on recadre les points
 996   // Same code as UClosed : optimisation souhaitable !!
 997   // CKY : d abord un code IDENTIQUE A UClosed; PUIS le special Seam ...
 998   // Si la surface est UCLosed, on recadre les points
 999   //
1000   //#69 rln 01.03.99 S4135 bm2_sd_t4-A.stp entity 30
1001   //#78 rln 12.03.99 S4135
1002   if (mySurf->IsVClosed(myPreci) || mySurf->Surface()->IsKind (STANDARD_TYPE (Geom_SphericalSurface))) {
1003     // Premier point dans le domain [vf, vl]
1004     Standard_Real prevY, firstY = pnt2d (1).Y();
1005     while (firstY < vf)  {  firstY += Vp;  pnt2d (1).SetY(firstY);  }
1006     while (firstY > vl)  {  firstY -= Vp;  pnt2d (1).SetY(firstY);  }
1007     prevY = firstY;
1008
1009     //:97 abv 1 Feb 98: treat case when curve is whole out of surface bounds
1010     Standard_Real minY = firstY, maxY = firstY;
1011
1012     // On decalle toujours le suivant
1013     for (i = 2; i <= nbrPnt; i ++) {
1014       //      dist2d = pnt2d (i-1).Distance(pnt2d (i));
1015       Standard_Real CurY = pnt2d (i).Y();
1016       dist2d = Abs (CurY - prevY);
1017       if (dist2d > ( Vp / 2) ) {
1018         if        (CurY > prevY + Vp/2) {
1019           while (CurY > prevY + Vp/2) {  CurY -= Vp;  pnt2d (i).SetY (CurY);  }
1020         } else if (CurY < prevY - Vp/2) {
1021           while (CurY < prevY - Vp/2) {  CurY += Vp;  pnt2d (i).SetY (CurY);  }
1022         }
1023       }
1024       prevY = CurY;
1025       if ( minY > CurY ) minY = CurY;      //:97
1026       else if ( maxY < CurY ) maxY = CurY; //:97
1027     }
1028
1029     //:97
1030     Standard_Real midY = 0.5 * ( minY + maxY );
1031     Standard_Real shiftY=0.;
1032     if ( midY > vl ) shiftY = -Vp;
1033     else if ( midY < vf ) shiftY = Vp;
1034     if ( shiftY != 0. )
1035       for ( i=1; i <= nbrPnt; i++ ) pnt2d(i).SetY ( pnt2d(i).Y() + shiftY );
1036   }
1037
1038   //#69 rln 01.03.99 S4135 bm2_sd_t4-A.stp entity 30
1039   //#78 rln 12.03.99 S4135
1040   if (mySurf->IsVClosed(myPreci) || mySurf->Surface()->IsKind (STANDARD_TYPE (Geom_SphericalSurface))) {
1041     for (i = 2; i <= nbrPnt; i++) {
1042       //#1 rln 11/02/98 ca_exhaust.stp entity #9869 dist2d = pnt2d (i-1).Distance(pnt2d (i));
1043       dist2d = Abs (pnt2d(i).Y() - pnt2d(i - 1).Y());
1044       if (dist2d > ( Vp / 2) ) {
1045         // ATTENTION : il faut regarder ou le decalage se fait.
1046         // si plusieurs points sont decalles, il faut plusieurs passes
1047         // pour obtenir un resultat correct.
1048         // NOT YET IMPLEMENTED
1049
1050         // one of those point is incorrectly placed
1051         // i.e on the wrong side of the "seam"
1052         // on prend le point le plus pres des bords vf ou vl
1053         Standard_Boolean prevOnFirst = Standard_False;
1054         Standard_Boolean prevOnLast  = Standard_False;
1055         Standard_Boolean currOnFirst = Standard_False;
1056         Standard_Boolean currOnLast  = Standard_False;
1057
1058         //  .X ?  plutot .Y ,  non ?
1059         Standard_Real distPrevVF = Abs(pnt2d (i-1).Y() - vf);
1060         Standard_Real distPrevVL = Abs(pnt2d (i-1).Y() - vl);
1061         Standard_Real distCurrVF = Abs(pnt2d (i).Y() - vf);
1062         Standard_Real distCurrVL = Abs(pnt2d (i).Y() - vl);
1063
1064         Standard_Real theMin = distPrevVF;
1065         prevOnFirst = Standard_True;
1066         if (distPrevVL < theMin) {
1067           theMin = distPrevVL;
1068           prevOnFirst = Standard_False;
1069           prevOnLast  = Standard_True;
1070         }
1071         if (distCurrVF < theMin) {
1072           theMin = distCurrVF;
1073           prevOnFirst = Standard_False;
1074           prevOnLast  = Standard_False;
1075           currOnFirst = Standard_True;
1076         }
1077         if (distCurrVL < theMin) {
1078           theMin = distCurrVL;
1079           prevOnFirst = Standard_False;
1080           prevOnLast  = Standard_False;
1081           currOnFirst = Standard_False;
1082           currOnLast  = Standard_True;
1083         }
1084         //  Modifs RLN/Nijni  3-DEC-1997
1085         if (prevOnFirst) {
1086           // on decalle le point (i-1) en V Last
1087           gp_Pnt2d newPrev(pnt2d (i-1).X(), vf); // instead of  vl RLN/Nijni
1088           pnt2d (i-1) = newPrev;
1089         }
1090         else if (prevOnLast) {
1091           // on decalle le point (i-1) en V first
1092           gp_Pnt2d newPrev(pnt2d (i-1).X(), vl); // instead of  vf RLN/Nijni
1093           pnt2d (i-1) = newPrev;
1094         }
1095         else if (currOnFirst) {
1096           // on decalle le point (i) en V Last
1097           gp_Pnt2d newCurr(pnt2d (i).X(),vf);  // instead of vl  RLN/Nijni
1098           pnt2d (i) = newCurr;
1099         }
1100         else if (currOnLast) {
1101           // on decalle le point (i) en V First
1102           gp_Pnt2d newCurr(pnt2d (i).X(), vl); // instead of vf  RLN/Nijni
1103           pnt2d (i) = newCurr;
1104         }
1105         // on verifie
1106 #ifdef OCCT_DEBUG
1107         dist2d = pnt2d (i-1).Distance(pnt2d (i));
1108         if (dist2d > ( Vp / 2) ) {
1109           cout << "Echec dans le recadrage" << endl;
1110         }
1111 #endif
1112       }
1113     }
1114   }
1115
1116   //:c0 abv 20 Feb 98: treat very special case when 3d curve
1117   // go over the pole of, e.g., sphere, and partly lies along seam.
1118   // 2d representation of such a curve should consist of 3 parts - one on
1119   // regular part of surface (interior), one part along degenerated boundary
1120   // and one along seam.
1121   // Since it cannot be adjusted later by arranging pcurves (curve is single),
1122   // to fix it it is nesessary to have a possibility of adjusting seam
1123   // part of such curve either to left or right boundary of surface.
1124   // Test is performed only if flag AdjustOverDegen is not -1.
1125   // If AdjustOverDegen is True, seam part of curve is adjusted to
1126   // the left, and if False - to the right parametric boundary
1127   // If treated case is detected, flag DONE4 is set to status
1128   // NOTE: currently, precision is Precision::PConfusion() since it
1129   // is enough on encountered example
1130   // (ug_turbine-A.stp from ProSTEP Benchmark #3, entities ##2470 & 5680)
1131   // (r1001_ac.stp from Test Rally #10, face #35027 and others)
1132   if ( myAdjustOverDegen != -1 ) {
1133     if ( mySurf->IsUClosed(myPreci) ) {//#78 rln 12.03.99 S4135
1134       mySurf->IsDegenerated ( gp_Pnt(0,0,0), myPreci );  // pour calculer les dgnr
1135       if ( mySurf->NbSingularities(myPreci) > 0 ) { //rln S4135
1136         // 1st, find gap point (degenerated pole)
1137         Standard_Real PrevX=0.;
1138         Standard_Integer OnBound=0, PrevOnBound=0;
1139         Standard_Integer ind; // svv #1
1140         Standard_Boolean start = Standard_True;
1141         for ( ind=1; ind <= nbrPnt; ind++ ) {
1142           Standard_Real CurX = pnt2d(ind).X();
1143           // abv 16 Mar 00: trj3_s1-ug.stp #697: ignore points in singularity
1144           if ( mySurf->IsDegenerated ( points(ind), Precision::Confusion() ) )
1145             continue;
1146           OnBound = ( Abs ( Abs ( CurX - 0.5 * ( ul + uf ) ) - Up/2 ) <=
1147                      Precision::PConfusion() );
1148           if ( ! start &&  Abs ( Abs ( CurX - PrevX ) - Up/2 ) <= 0.01*Up )
1149             break;
1150           start = Standard_False;
1151           PrevX = CurX;
1152           PrevOnBound = OnBound;
1153         }
1154         // if found, adjust seam part
1155         if ( ind <= nbrPnt ) {
1156           PrevX = ( myAdjustOverDegen ? uf : ul );
1157           Standard_Real dU = Up/2 + Precision::PConfusion();
1158           if ( PrevOnBound ) {
1159             pnt2d(ind-1).SetX ( PrevX );
1160             for ( Standard_Integer j=ind-2; j >0; j-- ) {
1161               Standard_Real CurX = pnt2d(j).X();
1162               while ( CurX < PrevX - dU ) pnt2d(j).SetX ( CurX += Up );
1163               while ( CurX > PrevX + dU ) pnt2d(j).SetX ( CurX -= Up );
1164             }
1165           }
1166           else if ( OnBound ) {
1167             pnt2d(ind).SetX ( PrevX );
1168             for ( Standard_Integer j=ind+1; j <= nbrPnt; j++ ) {
1169               Standard_Real CurX = pnt2d(j).X();
1170               while ( CurX < PrevX - dU ) pnt2d(j).SetX ( CurX += Up );
1171               while ( CurX > PrevX + dU ) pnt2d(j).SetX ( CurX -= Up );
1172             }
1173           }
1174           myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE4);
1175         }
1176       }
1177     }
1178     else if ( mySurf->IsVClosed(myPreci) ) {//#78 rln 12.03.99 S4135
1179       mySurf->IsDegenerated ( gp_Pnt(0,0,0), myPreci );  // pour calculer les dgnr
1180       if ( mySurf->NbSingularities(myPreci) > 0 ) { //rln S4135
1181         // 1st, find gap point (degenerated pole)
1182         Standard_Real PrevY=0.;
1183         Standard_Integer OnBound=0, PrevOnBound=0;
1184         Standard_Integer ind; // svv #1
1185         Standard_Boolean start = Standard_True;
1186         for ( ind=1; ind <= nbrPnt; ind++ ) {
1187           Standard_Real CurY = pnt2d(ind).Y();
1188           // abv 16 Mar 00: trj3_s1-ug.stp #697: ignore points in singularity
1189           if ( mySurf->IsDegenerated ( points(ind), Precision::Confusion() ) )
1190             continue;
1191           OnBound = ( Abs ( Abs ( CurY - 0.5 * ( vl + vf ) ) - Vp/2 ) <=
1192                      Precision::PConfusion() );
1193           if ( ! start &&  Abs ( Abs ( CurY - PrevY ) - Vp/2 ) <= 0.01*Vp )
1194             break;
1195           start = Standard_False;
1196           PrevY = CurY;
1197           PrevOnBound = OnBound;
1198         }
1199         // if found, adjust seam part
1200         if ( ind <= nbrPnt ) {
1201           PrevY = ( myAdjustOverDegen ? vf : vl );
1202           Standard_Real dV = Vp/2 + Precision::PConfusion();
1203           if ( PrevOnBound ) {
1204             pnt2d(ind-1).SetY ( PrevY );
1205             for ( Standard_Integer j=ind-2; j >0; j-- ) {
1206               Standard_Real CurY = pnt2d(j).Y();
1207               while ( CurY < PrevY - dV ) pnt2d(j).SetY ( CurY += Vp );
1208               while ( CurY > PrevY + dV ) pnt2d(j).SetY ( CurY -= Vp );
1209             }
1210           }
1211           else if ( OnBound ) {
1212             pnt2d(ind).SetY ( PrevY );
1213             for ( Standard_Integer j=ind+1; j <= nbrPnt; j++ ) {
1214               Standard_Real CurY = pnt2d(j).Y();
1215               while ( CurY < PrevY - dV ) pnt2d(j).SetY ( CurY += Vp );
1216               while ( CurY > PrevY + dV ) pnt2d(j).SetY ( CurY -= Vp );
1217             }
1218           }
1219           myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE4);
1220         }
1221       }
1222     }
1223   }
1224
1225   //:q9: fill cashe
1226   myNbCashe = 2;
1227   if(ChangeCycle) {  // msv 10.08.04: avoid using of uninitialised field
1228   //if(myCashe3d[0].Distance(points(1))>Precision::Confusion() &&
1229   //   myCashe3d[1].Distance(points(1))>Precision::Confusion()) {
1230     myCashe3d[0] = points(1);
1231     myCashe3d[1] = points(nbrPnt);
1232     myCashe2d[0] = pnt2d(1);
1233     myCashe2d[1] = pnt2d(nbrPnt);
1234   }
1235   else {
1236     myCashe3d[1] = points(1);
1237     myCashe3d[0] = points(nbrPnt);
1238     myCashe2d[1] = pnt2d(1);
1239     myCashe2d[0] = pnt2d(nbrPnt);
1240   }
1241   return isDone;
1242 }
1243
1244 //=======================================================================
1245 //function : ApproximatePCurve
1246 //purpose  :
1247 //=======================================================================
1248
1249 Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(const Standard_Integer /*nbrPnt*/,
1250                                                                              Handle(TColgp_HArray1OfPnt2d)& points2d,
1251                                                                              Handle(TColStd_HArray1OfReal)& params,
1252                                                                              const Handle(Geom_Curve)& /*orig*/) const
1253 {
1254 //  Standard_Real resol = Min(mySurf->Adaptor3d()->VResolution(myPreci), mySurf->Adaptor3d()->UResolution(myPreci));
1255   Standard_Real theTolerance2d = myPreci; // (100*nbrPnt);//resol;
1256   Handle(Geom2d_Curve) C2d;
1257   try {
1258     OCC_CATCH_SIGNALS
1259     CheckPoints2d (points2d, params, theTolerance2d);
1260     Standard_Integer numberPnt = points2d->Length();
1261
1262     TColgp_Array1OfPnt points3d(1,numberPnt);
1263     gp_Pnt2d pnt2d;
1264     gp_Pnt pnt;
1265     Standard_Integer i; // svv #1
1266     for( i = 1; i <= numberPnt; i++) {
1267       pnt2d = points2d->Value(i);
1268       pnt.SetCoord(pnt2d.X(),pnt2d.Y(),0);
1269       points3d(i) = pnt;
1270     }
1271
1272     GeomAPI_PointsToBSpline appr(points3d, params->Array1(), 1, 10, GeomAbs_C1, theTolerance2d);
1273     Handle(Geom_BSplineCurve) crv3d = appr.Curve();
1274     Standard_Integer NbPoles = crv3d->NbPoles();
1275     TColgp_Array1OfPnt poles3d (1, NbPoles);
1276     TColgp_Array1OfPnt2d poles2d (1, NbPoles);
1277     crv3d->Poles(poles3d);
1278     for( i = 1; i <= NbPoles; i++) {
1279       pnt2d.SetCoord(poles3d(i).X(),poles3d(i).Y());
1280       poles2d(i) = pnt2d;
1281     }
1282     TColStd_Array1OfReal weights (1,NbPoles);
1283     TColStd_Array1OfInteger multiplicities (1,crv3d->NbKnots());
1284     TColStd_Array1OfReal knots(1,crv3d->NbKnots());
1285     crv3d->Knots(knots);
1286     crv3d->Weights(weights);
1287     crv3d->Multiplicities(multiplicities);
1288     C2d = new Geom2d_BSplineCurve  ( poles2d, weights, knots, multiplicities, crv3d->Degree(), crv3d->IsPeriodic());
1289     return C2d;
1290   }
1291   catch(Standard_Failure) {
1292 #ifdef OCCT_DEBUG //:s5
1293     //    debug ...
1294     Standard_Integer nbp = params->Length();
1295     Standard_Integer nb2 = points2d->Length();
1296     cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(): Exception: ";
1297     Standard_Failure::Caught()->Print(cout);
1298     cout<<"Pb Geom2dAPI_Approximate, tol2d="<<theTolerance2d<<" NbParams="<<nbp<<" NbPnts="<<nb2<<endl;
1299 //     if (nb2 > nbp) nb2 = nbp;
1300 //     Standard_Real rbp,rb2; rbp = nbp; rb2 = nb2;
1301 //     //    dbl.AddString ("NbP2d/NbParams puis  X Y Param -> mini");
1302 //     dbl.AddReals (rb2,rbp);
1303 //     for (Standard_Integer i = 1; i <= nb2; i ++) {
1304 //       gp_XYZ quoi (points2d->Value(i).X(),points2d->Value(i).Y(),params->Value(i) );
1305 //       dbl.AddXYZ (quoi);
1306 //     }
1307 #endif
1308      C2d.Nullify();
1309   }
1310   return C2d;
1311 }
1312
1313 //=======================================================================
1314 //function : InterpolatePCurve
1315 //purpose  :
1316 //=======================================================================
1317
1318  Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolatePCurve(const Standard_Integer nbrPnt,
1319                                                                               Handle(TColgp_HArray1OfPnt2d)& points2d,
1320                                                                               Handle(TColStd_HArray1OfReal)& params,
1321                                                                               const Handle(Geom_Curve)& /*orig*/) const
1322 {
1323   Handle(Geom2d_Curve) C2d;    // NULL si echec
1324   Standard_Real theTolerance2d = myPreci / (100 * nbrPnt);
1325   try {
1326     OCC_CATCH_SIGNALS
1327     // on verifie d abord s il n y a pas de points confondus
1328     // si besoin on retouche les valeurs ...
1329     CheckPoints2d (points2d, params, theTolerance2d);
1330     Geom2dAPI_Interpolate myInterPol2d (points2d, params,
1331                                         Standard_False, theTolerance2d);
1332     myInterPol2d.Perform();
1333     if (myInterPol2d.IsDone()) C2d = myInterPol2d.Curve();
1334   }
1335   catch(Standard_Failure) {
1336 #ifdef OCCT_DEBUG //:s5
1337 // //    debug ...
1338     Standard_Integer nbp = params->Length();
1339     Standard_Integer nb2 = points2d->Length();
1340     cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::InterpolatePCurve(): Exception: ";
1341     Standard_Failure::Caught()->Print(cout);
1342     cout<<"Pb Geom2dAPI_Interpolate, tol2d="<<theTolerance2d<<" NbParams="<<nbp<<" NbPnts="<<nb2<<endl;
1343 //     if (nb2 > nbp) nb2 = nbp;
1344 //     Standard_Real rbp,rb2; rbp = nbp; rb2 = nb2;
1345 // //    dbl.AddString ("NbP2d/NbParams puis  X Y Param -> mini");
1346 //     dbl.AddReals (rb2,rbp);
1347 //     for (Standard_Integer i = 1; i <= nb2; i ++) {
1348 //       gp_XYZ quoi (points2d->Value(i).X(),points2d->Value(i).Y(),params->Value(i) );
1349 //       dbl.AddXYZ (quoi);
1350 //     }
1351 #endif
1352     C2d.Nullify();
1353   }
1354   return C2d;
1355 }
1356
1357 //=======================================================================
1358 //function : InterpolateCurve3d
1359 //purpose  :
1360 //=======================================================================
1361
1362 Handle(Geom_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(const Standard_Integer,
1363                                                                             Handle(TColgp_HArray1OfPnt)& points,
1364                                                                             Handle(TColStd_HArray1OfReal)& params,
1365                                                                             const Handle(Geom_Curve)& /*orig*/) const
1366 {
1367   Handle(Geom_Curve) C3d;    // NULL si echec
1368   try {
1369     OCC_CATCH_SIGNALS
1370     Standard_Real Tol = myPreci;
1371     CheckPoints(points, params, Tol);
1372     GeomAPI_Interpolate myInterPol(points, params, Standard_False, Tol);
1373     myInterPol.Perform();
1374     if (myInterPol.IsDone()) C3d = myInterPol.Curve();
1375   }
1376   catch(Standard_Failure) {
1377     C3d.Nullify();
1378 #ifdef OCCT_DEBUG //:s5
1379     cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(): Exception: ";
1380     Standard_Failure::Caught()->Print(cout); cout << endl;
1381 #endif
1382   }
1383   return C3d;
1384 }
1385
1386 //=======================================================================
1387 //function : CheckPoints
1388 //purpose  :
1389 //=======================================================================
1390
1391  void ShapeConstruct_ProjectCurveOnSurface::CheckPoints(Handle(TColgp_HArray1OfPnt)& points,Handle(TColStd_HArray1OfReal)& params,Standard_Real& preci) const
1392 {
1393   Standard_Integer firstElem = points->Lower();
1394   Standard_Integer lastElem  = points->Upper();
1395   Standard_Integer i;
1396   Standard_Integer nbPntDropped = 0;
1397   Standard_Integer lastValid = firstElem; // indice of last undropped point
1398
1399   // will store 0 when the point is to be removed, 1 otherwise
1400   TColStd_Array1OfInteger tmpParam(firstElem, lastElem);
1401   for (i = firstElem; i<=lastElem ; i++) tmpParam.SetValue(i,1);
1402   Standard_Real DistMin2 = RealLast();
1403   gp_Pnt Prev = points->Value (lastValid);
1404   gp_Pnt Curr;
1405   for (i = firstElem + 1; i <= lastElem ; i ++) {
1406     Curr = points->Value(i);
1407     Standard_Real CurDist2 = Prev.SquareDistance(Curr);
1408     if (CurDist2 < gp::Resolution()) {  // test 0
1409       nbPntDropped ++;
1410       if ( i == lastElem ) tmpParam.SetValue(lastValid, 0); // last point kept
1411       else tmpParam.SetValue(i, 0);    // current dropped, lastValid unchanged
1412     } else {
1413       if (CurDist2 < DistMin2)
1414         DistMin2 = CurDist2;
1415       // lastValid becomes the current (i.e. i)
1416       lastValid = i;
1417       Prev = Curr;
1418     }
1419   }
1420   if (DistMin2 < RealLast())
1421     preci = 0.9 * Sqrt (DistMin2); // preci est la distance min entre les points on la reduit un peu
1422   if (nbPntDropped == 0)
1423     return;
1424
1425 #ifdef OCCT_DEBUG
1426   cout << "Warning : removing 3d points for interpolation" << endl;
1427 #endif
1428   // Build new HArrays
1429   Standard_Integer newLast = lastElem - nbPntDropped;
1430   if ((newLast - firstElem + 1)  < 2) {
1431 #ifdef OCCT_DEBUG
1432     cout << "Too many degenerated points for 3D interpolation" << endl;
1433 #endif
1434     return;
1435   }
1436   Handle(TColgp_HArray1OfPnt) newPnts =
1437     new TColgp_HArray1OfPnt(firstElem, newLast);
1438   Handle(TColStd_HArray1OfReal) newParams =
1439     new TColStd_HArray1OfReal(firstElem, newLast);
1440   Standard_Integer newCurr = 1;
1441   for (i = firstElem; i<= lastElem ; i++) {
1442     if (tmpParam.Value(i) == 1) {
1443       newPnts->SetValue(newCurr, points->Value(i));
1444       newParams->SetValue(newCurr, params->Value(i));
1445       newCurr ++;
1446     }
1447   }
1448   points = newPnts;
1449   params = newParams;
1450   // on la reduit un peu
1451 }
1452
1453 //=======================================================================
1454 //function : CheckPoints2d
1455 //purpose  :
1456 //=======================================================================
1457
1458  void ShapeConstruct_ProjectCurveOnSurface::CheckPoints2d(Handle(TColgp_HArray1OfPnt2d)& points,
1459                                                           Handle(TColStd_HArray1OfReal)& params,
1460                                                           Standard_Real& preci) const
1461 {
1462   Standard_Integer firstElem = points->Lower();
1463   Standard_Integer lastElem  = points->Upper();
1464   Standard_Integer i;
1465   Standard_Integer nbPntDropped = 0;
1466   Standard_Integer lastValid = firstElem; // indice of last undropped point
1467
1468   // will store 0 when the point is to be removed, 1 otherwise
1469   TColStd_Array1OfInteger tmpParam(firstElem, lastElem);
1470   for (i = firstElem; i<=lastElem ; i++) {
1471     tmpParam.SetValue(i,1);
1472   }
1473   Standard_Real DistMin2 = RealLast();
1474   gp_Pnt2d Prev = points->Value(lastValid);
1475   gp_Pnt2d Curr;
1476   for (i = firstElem + 1; i<=lastElem ; i++) {
1477     Curr = points->Value(i);
1478     Standard_Real CurDist2 = Prev.SquareDistance(Curr);
1479     if (CurDist2 < gp::Resolution()) {  // test 0
1480       nbPntDropped ++;
1481       if ( i == lastElem ) tmpParam.SetValue(lastValid, 0); // last point kept
1482       else tmpParam.SetValue(i, 0);    // current dropped, lastValid unchanged
1483     } else {
1484       if (CurDist2 < DistMin2)
1485         DistMin2 = CurDist2;
1486       // lastValid becomes the current (i.e. i)
1487       lastValid = i;
1488       Prev = Curr;
1489     }
1490   }
1491   if (DistMin2 < RealLast())
1492     preci = 0.9 * Sqrt (DistMin2);
1493   if (nbPntDropped == 0)
1494     return;
1495
1496 #ifdef OCCT_DEBUG
1497   cout << "Warning : removing 2d points for interpolation" << endl;
1498 #endif
1499   // Build new HArrays
1500   Standard_Integer newLast = lastElem - nbPntDropped;
1501   if ((newLast - firstElem + 1)  < 2) {
1502 #ifdef OCCT_DEBUG
1503     cout << "Too many degenerated points for 2D interpolation" << endl;
1504 #endif
1505     //pdn 12.02.99 S4135 Creating pcurve with minimal length.
1506     tmpParam.SetValue(firstElem,1);
1507     tmpParam.SetValue(lastElem,1);
1508     gp_XY  lastPnt = points->Value(lastElem).XY();
1509     lastPnt.Add(gp_XY(preci,preci));
1510     points->SetValue(lastElem,lastPnt);
1511     newLast = firstElem+1;
1512     //return;
1513   }
1514   Handle(TColgp_HArray1OfPnt2d) newPnts =
1515     new TColgp_HArray1OfPnt2d(firstElem, newLast);
1516   Handle(TColStd_HArray1OfReal) newParams =
1517     new TColStd_HArray1OfReal(firstElem, newLast);
1518   Standard_Integer newCurr = 1;
1519   for (i = firstElem; i <= lastElem ; i++) {
1520     if (tmpParam.Value(i) == 1) {
1521 #ifdef OCCT_DEBUG
1522       cout << "Point " << i << " : " << points->Value(i).X() << " " << points->Value(i).Y() << " at param " <<  params->Value(i) << endl;
1523 #endif
1524       newPnts->SetValue(newCurr, points->Value(i));
1525       newParams->SetValue(newCurr, params->Value(i));
1526       newCurr ++;
1527     }
1528     else {
1529 #ifdef OCCT_DEBUG
1530       cout << "Removed " << i << " : " << points->Value(i).X() << " " << points->Value(i).Y() << " at param " <<  params->Value(i) << endl;
1531 #endif
1532     }
1533   }
1534   points = newPnts;
1535   params = newParams;
1536 }
1537
1538 //=======================================================================
1539 //function : IsAnIsoparametric
1540 //purpose  :
1541 //=======================================================================
1542 //:S4030: modified for optimization
1543 //:p9 abv 11 Mar 99: PRO7226 #489490: find nearest boundary instead of first one
1544
1545  Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::IsAnIsoparametric(const Standard_Integer nbrPnt,
1546                                                                           const TColgp_Array1OfPnt& points,
1547                                                                           const TColStd_Array1OfReal& params,
1548                                                                           Standard_Boolean& isoTypeU,
1549                                                                           Standard_Boolean& p1OnIso,
1550                                                                           gp_Pnt2d& valueP1,
1551                                                                           Standard_Boolean& p2OnIso,
1552                                                                           gp_Pnt2d& valueP2,
1553                                                                           Standard_Boolean& isoPar2d3d,
1554                                                                           Handle(Geom_Curve)& cIso,
1555                                                                           Standard_Real& t1,
1556                                                                           Standard_Real& t2,
1557                                                                           TColStd_Array1OfReal& pout) const
1558 {
1559   try {    // RAJOUT
1560     OCC_CATCH_SIGNALS
1561
1562   Standard_Real prec = Precision::Confusion();//myPreci;
1563
1564   Standard_Boolean isoParam = Standard_False;
1565   isoPar2d3d = Standard_False;
1566
1567   Standard_Real U1, U2, V1, V2;
1568   mySurf->Bounds(U1, U2, V1, V2);
1569
1570   if ( mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
1571     Handle(Geom_RectangularTrimmedSurface) sTrim =
1572       Handle(Geom_RectangularTrimmedSurface)::DownCast(mySurf->Surface());
1573     sTrim->Bounds(U1, U2, V1, V2);
1574   }
1575
1576   gp_Pnt pt;
1577   Standard_Integer mpt[2]; mpt[0] = mpt[1] = 0;
1578   Standard_Real t, tpar[2] = { 0.0, 0.0 }, isoValue=0.;
1579   Standard_Real mindist2;
1580   Standard_Real mind2[2];
1581   mindist2 = mind2[0] = mind2[1] = 4*prec*prec;
1582
1583   p1OnIso = Standard_False;
1584   p2OnIso = Standard_False;
1585   const Bnd_Box* aBox = 0;
1586
1587   for (Standard_Integer j=1; (j<=4) /*&& !isoParam*/; j++) {
1588     Standard_Real isoVal=0.;
1589     Standard_Boolean isoU=Standard_False; //szv#4:S4163:12Mar99 `isoU` must be Standard_Boolean
1590     Handle(Geom_Curve) cI;
1591     Standard_Real tt1, tt2;
1592
1593     if      (j == 1 ) {
1594       if (Precision::IsInfinite(U1)) continue;
1595       cI = mySurf->UIso(U1);
1596       isoU = Standard_True;
1597       isoVal = U1;
1598       aBox = & mySurf->GetBoxUF();
1599     }
1600     else if (j == 2) {
1601       if (Precision::IsInfinite(U2)) continue;
1602       cI = mySurf->UIso(U2);
1603       isoU = Standard_True;
1604       isoVal = U2;
1605       aBox = & mySurf->GetBoxUL();
1606     }
1607     else if (j == 3) {
1608       if (Precision::IsInfinite(V1)) continue;
1609       cI = mySurf->VIso(V1);
1610       isoU = Standard_False;
1611       isoVal = V1;
1612       aBox = & mySurf->GetBoxVF();
1613     }
1614     else if (j == 4) {
1615       if (Precision::IsInfinite(V2)) continue;
1616       cI = mySurf->VIso(V2);
1617       isoU = Standard_False;
1618       isoVal = V2;
1619       aBox = & mySurf->GetBoxVL();
1620     }
1621     if(cI.IsNull())
1622       continue;
1623
1624     if (isoU)  {  tt1 = V1;  tt2 = V2;  }
1625     else       {  tt1 = U1;  tt2 = U2;  }
1626
1627     gp_Pnt ext1, ext2;
1628     cI->D0(tt1, ext1);
1629     cI->D0(tt2, ext2);
1630
1631 // PATCH CKY 9-JUL-1998 : protection contre singularite
1632     gp_Pnt extmi;
1633     cI->D0( (tt1+tt2)/2,extmi);
1634     if (ext1.IsEqual(ext2,prec) && ext1.IsEqual(extmi,prec)) continue;
1635
1636     Standard_Boolean PtEQext1 = Standard_False;
1637     Standard_Boolean PtEQext2 = Standard_False;
1638
1639     Standard_Real currd2[2], tp[2] = {0, 0};
1640     Standard_Integer mp[2];
1641
1642     for (Standard_Integer i=0; i<2; i++) {
1643       mp[i] = 0;
1644       Standard_Integer k = (i == 0 ? 1 : nbrPnt);
1645
1646       // si ext1 == ext2 => valueP1 == valueP2 => vect null plus tard
1647       currd2[i] = points(k).SquareDistance ( ext1 );
1648       if ( currd2[i] <= prec*prec && !PtEQext1) {
1649         mp[i] = 1;
1650         tp[i] = tt1;
1651         PtEQext1 = Standard_True;
1652         continue;
1653       }
1654
1655       currd2[i] = points(k).SquareDistance ( ext2 );
1656       if ( currd2[i] <= prec*prec && !PtEQext2) {
1657         mp[i] = 2;
1658         tp[i] = tt2;
1659         PtEQext2 = Standard_True;
1660         continue;
1661       }
1662
1663       // On evite de projecter sur un iso degenere
1664       // on doit egalement le faire pour l apex du cone
1665       if (mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) && !isoU) {
1666         continue;
1667       }
1668
1669       if(aBox->IsOut(points(k))) continue;
1670
1671       Standard_Real Cf = cI->FirstParameter();
1672       Standard_Real Cl = cI->LastParameter();
1673       if (Precision::IsInfinite(Cf))  Cf = -1000;
1674       if (Precision::IsInfinite(Cl))  Cl = +1000;
1675
1676       ShapeAnalysis_Curve sac;
1677       Standard_Real dist = sac.Project (cI,points(k),prec,pt,t,Cf,Cl);
1678       currd2[i] = dist * dist;
1679       if ((dist <= prec) && (t>= Cf) && (t<=Cl)) {
1680         mp[i]  = 3;
1681         tp[i] = t;
1682       }
1683     }
1684
1685     //:e7 abv 21 Apr 98: ProSTEP TR8, r0501_pe #56679:
1686     // avoid possible null-length curves
1687     if ( mp[0] >0 && mp[1] >0 &&
1688          Abs ( tp[0] - tp[1] ) < Precision::PConfusion() ) continue;
1689
1690
1691     if (mp[0] > 0 &&
1692         ( ! p1OnIso || currd2[0] < mind2[0] ) ) {
1693       p1OnIso = Standard_True;
1694       mind2[0] = currd2[0]; // LP2.stp #105899: FLT_INVALID_OPERATION on Windows 7 VC 9 Release mode on the whole file
1695       if (isoU) valueP1.SetCoord(isoVal, tp[0]);
1696       else      valueP1.SetCoord(tp[0], isoVal);
1697     }
1698
1699     if (mp[1] > 0 &&
1700         ( ! p2OnIso || currd2[1] < mind2[1] ) ) {
1701       p2OnIso = Standard_True;
1702       mind2[1] = currd2[1];
1703       if (isoU) valueP2.SetCoord(isoVal, tp[1]);
1704       else      valueP2.SetCoord(tp[1], isoVal);
1705     }
1706
1707     if ( mp[0] <=0 || mp[1] <=0 ) continue;
1708
1709     Standard_Real md2 = currd2[0] + currd2[1];
1710     if ( mindist2 <= md2 ) continue;
1711
1712     mindist2 = md2;
1713     mpt[0] = mp[0];
1714     mpt[1] = mp[1];
1715     tpar[0] = tp[0];
1716     tpar[1] = tp[1];
1717     isoTypeU = isoU;
1718     isoValue = isoVal;
1719     cIso = cI;
1720     t1 = tt1;
1721     t2 = tt2;
1722   }
1723
1724   // probablely it concerns an isoparametrics
1725   if ( mpt[0] >0 && mpt[1] >0 ) {
1726
1727     p1OnIso = p2OnIso = Standard_True;
1728     if (isoTypeU) {
1729       valueP1.SetCoord(isoValue, tpar[0]);
1730       valueP2.SetCoord(isoValue, tpar[1]);
1731     }
1732     else {
1733       valueP1.SetCoord(tpar[0], isoValue);
1734       valueP2.SetCoord(tpar[1], isoValue);
1735     }
1736
1737     if ( mpt[0] != 3 && mpt[1] != 3 ) {
1738       isoPar2d3d = Standard_True;
1739       for (Standard_Integer i=2; i < nbrPnt && isoPar2d3d; i++){
1740         if (tpar[1] > tpar[0])  t = params(i);
1741         else                    t = t1+t2-params(i);
1742         cIso->D0(t, pt);
1743         if (!points(i).IsEqual(pt, prec)) isoPar2d3d = Standard_False;
1744       }
1745     }
1746
1747     if (isoPar2d3d) isoParam = Standard_True;
1748     else {
1749       Standard_Real prevParam = tpar[0];
1750       Standard_Real Cf, Cl;
1751       Standard_Boolean isoByDistance = Standard_True;
1752       Cf = cIso->FirstParameter();
1753       Cl = cIso->LastParameter();
1754       if (Precision::IsInfinite(Cf))  Cf = -1000;
1755       if (Precision::IsInfinite(Cl))  Cl = +1000;
1756
1757       ShapeAnalysis_Curve sac;
1758       for (Standard_Integer i=2; i < nbrPnt && isoByDistance; i++) {
1759         Standard_Real dist = sac.NextProject (prevParam,cIso,points(i),
1760                                               prec,pt,t,Cf,Cl,
1761                                               Standard_False); //:j8 abv 10.12.98: TR10 r0501_db.stp #9423: avoid adjusting to ends
1762         prevParam = t;
1763         pout(i)=t;
1764         if( (dist > prec) || (t < Cf) || (t > Cl) )
1765           isoByDistance = Standard_False;
1766       }
1767       if (isoByDistance) isoParam = Standard_True;
1768     }
1769   }
1770 /*  if (!isoParam) {    CKY 29-mai-1997 : garder tout ce qu on peut ?
1771     p1OnIso = Standard_False;
1772     p2OnIso = Standard_False;
1773   }  */
1774   return isoParam;
1775   }  // RAJOUT
1776   catch(Standard_Failure) {
1777 //  pb : on affiche ce qu on peut
1778 #ifdef OCCT_DEBUG
1779     for (Standard_Integer numpnt = 1; numpnt <= nbrPnt; numpnt ++) {
1780       cout<<"["<<numpnt<<"]param="<<params(numpnt)<<" point=("<<
1781         points(numpnt).X()<<"  "<<points(numpnt).Y()<<"  "<<points(numpnt).Z()<<")"<<endl;
1782     }
1783 #endif
1784 #ifdef OCCT_DEBUG //:s5
1785     cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::IsAnIsoparametric(): Exception: ";
1786     Standard_Failure::Caught()->Print(cout); cout << endl;
1787 #endif
1788     return Standard_False;
1789   }
1790 }
1791
1792 /* S4135 : BestExtremum is commented after IsAnIsoparametric works with Precision::Confusion()
1793 //=======================================================================
1794 //function : BestExtremum
1795 //purpose  : auxiliaire prenant le meilleur extremum si ISO car doute possible
1796 //=======================================================================
1797
1798  gp_Pnt2d ShapeConstruct_ProjectCurveOnSurface::BestExtremum(const gp_Pnt2d& P2iso,const gp_Pnt& P3ext,const gp_Pnt& P3next) const
1799 {
1800 //  P2iso a ete calcule depuis P3ext sur une iso externe de la surface
1801 //  En principe bon mais circularite possible ... et IsU/VClosed faillible
1802 //    (si baillement 1e-4 ou 1e-5, on est dedans !). DONC
1803 //  1/ on privilegie l iso mais a tout hasard on verifie si Surf meilleur
1804 //  2/ si iso, attention a la circularite (cas limite)
1805
1806 //  NB : si isoParam, on suppose que P2iso est bon (car il y en a 2). A voir...
1807
1808 //  D abord, calcul p2ext depuis la surface. choix surface/iso
1809   return P2iso;
1810   Standard_Real prec = Precision::Confusion();//myPreci;
1811   gp_Pnt2d P2cal = mySurf->ValueOfUV(P3ext, prec);
1812   gp_Pnt   P3cal = mySurf->Value (P2cal);
1813   Standard_Real dcal = P3ext.Distance (P3cal);
1814   Standard_Real dnxt = P3ext.Distance (P3next);
1815   if (dcal > dnxt) return P2iso;    // en fait protection sur BUG (PRO8468)
1816
1817 //  On choisit entre P2iso et P2cal, le plus proche de P2next ... !!!
1818   gp_Pnt2d P2next = mySurf->ValueOfUV(P3next, prec);
1819   if (P2next.Distance(P2cal) < P2next.Distance(P2iso)) return P2cal;
1820   return P2iso;
1821 }
1822 */