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 <ElCLib.hxx>
  62 #include <NCollection_Sequence.hxx>
  63 #include <Precision.hxx>
  64 #include <ProjLib_CompProjectedCurve.hxx>
  65 #include <ProjLib_HCompProjectedCurve.hxx>
  66 #include <ProjLib_ProjectedCurve.hxx>
  67 #include <ShapeAnalysis.hxx>
  68 #include <ShapeAnalysis_Curve.hxx>
  69 #include <ShapeAnalysis_Surface.hxx>
  70 #include <ShapeConstruct_ProjectCurveOnSurface.hxx>
  71 #include <ShapeExtend.hxx>
  72 #include <Standard_ErrorHandler.hxx>
  73 #include <Standard_Failure.hxx>
  74 #include <Standard_Type.hxx>
  75 #include <TColgp_Array1OfPnt.hxx>
  76 #include <TColStd_Array1OfInteger.hxx>
  77 #include <IntRes2d_Domain.hxx>
  78 #include <IntCurve_IntConicConic.hxx>
  79
  80 #include <algorithm>
  81 IMPLEMENT_STANDARD_RTTIEXT(ShapeConstruct_ProjectCurveOnSurface,Standard_Transient)
  82
  83 #define NCONTROL 23
  84
  85
  86 static void AdjustSecondPointToFirstPoint(const gp_Pnt2d& theFirstPoint,
  87                                           gp_Pnt2d& theSecondPoint,
  88                                           const Handle(Geom_Surface)& theSurf)
  89 {
  90   if (theSurf->IsUPeriodic())
  91   {
  92     Standard_Real UPeriod = theSurf->UPeriod();
  93     Standard_Real NewU = ElCLib::InPeriod(theSecondPoint.X(),
  94                                           theFirstPoint.X() - UPeriod/2,
  95                                           theFirstPoint.X() + UPeriod/2);
  96     theSecondPoint.SetX(NewU);
  97   }
  98   if (theSurf->IsVPeriodic())
  99   {
 100     Standard_Real VPeriod = theSurf->VPeriod();
 101     Standard_Real NewV = ElCLib::InPeriod(theSecondPoint.Y(),
 102                                           theFirstPoint.Y() - VPeriod/2,
 103                                           theFirstPoint.Y() + VPeriod/2);
 104     theSecondPoint.SetY(NewV);
 105   }
 106 }
 107
 108
 109 //=======================================================================
 110 //function : ShapeConstruct_ProjectCurveOnSurface
 111 //purpose  :
 112 //=======================================================================
 113
 114 ShapeConstruct_ProjectCurveOnSurface::ShapeConstruct_ProjectCurveOnSurface()
 115 {
 116   myPreci = Precision::Confusion();
 117   myBuild = Standard_False;
 118   myAdjustOverDegen = 1; //:c0 //szv#4:S4163:12Mar99 was boolean
 119   myNbCashe = 0; //:q9
 120 }
 121
 122 //=======================================================================
 123 //function : Init
 124 //purpose  :
 125 //=======================================================================
 126
 127  void ShapeConstruct_ProjectCurveOnSurface::Init(const Handle(Geom_Surface)& surf,const Standard_Real preci)
 128 {
 129   Init (new ShapeAnalysis_Surface (surf), preci);
 130 }
 131
 132 //=======================================================================
 133 //function : Init
 134 //purpose  :
 135 //=======================================================================
 136
 137  void ShapeConstruct_ProjectCurveOnSurface::Init(const Handle(ShapeAnalysis_Surface)& surf,const Standard_Real preci)
 138 {
 139   SetSurface (surf);
 140   SetPrecision (preci);
 141 }
 142
 143 //=======================================================================
 144 //function : SetSurface
 145 //purpose  :
 146 //=======================================================================
 147
 148  void ShapeConstruct_ProjectCurveOnSurface::SetSurface(const Handle(Geom_Surface)& surf)
 149 {
 150   SetSurface (new ShapeAnalysis_Surface (surf));
 151 }
 152
 153 //=======================================================================
 154 //function : SetSurface
 155 //purpose  :
 156 //=======================================================================
 157
 158  void ShapeConstruct_ProjectCurveOnSurface::SetSurface(const Handle(ShapeAnalysis_Surface)& surf)
 159 {
 160   if ( mySurf == surf ) return;
 161   mySurf = surf;
 162   myNbCashe = 0; //:q9
 163 }
 164
 165 //=======================================================================
 166 //function : SetPrecision
 167 //purpose  :
 168 //=======================================================================
 169
 170  void ShapeConstruct_ProjectCurveOnSurface::SetPrecision(const Standard_Real preci)
 171 {
 172   myPreci = preci;
 173 }
 174
 175 //=======================================================================
 176 //function : BuildCurveMode
 177 //purpose  :
 178 //=======================================================================
 179
 180  Standard_Boolean& ShapeConstruct_ProjectCurveOnSurface::BuildCurveMode()
 181 {
 182   return myBuild;
 183 }
 184
 185 //=======================================================================
 186 //function : AdjustOverDegenMode
 187 //purpose  :
 188 //=======================================================================
 189 //:c0
 190
 191 //szv#4:S4163:12Mar99 was Boolean
 192  Standard_Integer& ShapeConstruct_ProjectCurveOnSurface::AdjustOverDegenMode()
 193  {
 194    return myAdjustOverDegen;
 195  }
 196
 197 //=======================================================================
 198 //function : Status
 199 //purpose  :
 200 //=======================================================================
 201
 202  Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::Status (const ShapeExtend_Status theStatus) const
 203 {
 204   return ShapeExtend::DecodeStatus (myStatus, theStatus);
 205 }
 206
 207 //=======================================================================
 208 //function : Perform
 209 //purpose  :
 210 //=======================================================================
 211 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::Perform (Handle(Geom_Curve)& c3d,
 212                                                                 const Standard_Real First,
 213                                                                 const Standard_Real Last,
 214                                                                 Handle(Geom2d_Curve)& c2d,
 215                                                                 const Standard_Real TolFirst,
 216                                                                 const Standard_Real TolLast)
 217 {
 218   myStatus = ShapeExtend::EncodeStatus (ShapeExtend_OK);
 219   //Standard_Boolean OK = Standard_True; //szv#4:S4163:12Mar99 not needed
 220
 221   if (mySurf.IsNull()) {
 222     c2d.Nullify();
 223     myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_FAIL1);
 224     return Standard_False;
 225   }
 226 //  Projection Analytique
 227   Handle(Geom_Curve) crv3dtrim = c3d;
 228   if ( ! c3d->IsKind(STANDARD_TYPE(Geom_BoundedCurve)) )
 229     crv3dtrim = new Geom_TrimmedCurve ( c3d, First, Last );
 230   c2d = ProjectAnalytic ( crv3dtrim );
 231   if (!c2d.IsNull()) {
 232     myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE1);
 233     return Standard_True;
 234   }
 235
 236 //  Projection par approximation
 237
 238   // discretize the 3d curve
 239
 240   Standard_Integer nbrPnt;
 241
 242 //   $$$$    :92 abv 28 Jan 98   see PRO10107, big BSplineCurve C0
 243   Standard_Integer nbPini = NCONTROL;  // as in BRepCheck_Edge (RLN/Nijni)
 244  // 20; // number of points for interpolation, should be "parametric dependent"
 245
 246   //:92 abv 28 Jan 98: if curve is BSpline with many intervals,
 247   // increase number of points to provide at least Degree()+1 points per interval
 248   Handle(Geom_BSplineCurve) bspl;
 249   if ( c3d->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)) ) {
 250     Handle(Geom_TrimmedCurve) ctrim = Handle(Geom_TrimmedCurve)::DownCast(c3d);
 251     bspl = Handle(Geom_BSplineCurve)::DownCast ( ctrim->BasisCurve() );
 252   }
 253   else bspl = Handle(Geom_BSplineCurve)::DownCast ( c3d );
 254   if ( ! bspl.IsNull() ) {
 255     Standard_Integer nint = 0;
 256     for ( Standard_Integer i=1; i < bspl->NbKnots(); i++ )
 257     {
 258       if ( bspl->Knot(i+1) > First && bspl->Knot(i) < Last )
 259         nint++;
 260     }
 261     Standard_Integer minPnt = nint * ( bspl->Degree() + 1 );
 262     while ( nbPini < minPnt ) nbPini += NCONTROL - 1;
 263 #ifdef OCCT_DEBUG
 264     if ( nbPini > NCONTROL )
 265       std::cout << "Warning: number of points for projecting is " << nbPini << std::endl;
 266 #endif
 267   }
 268
 269 //    $$$$    end :92 (big BSplineCurve C0)
 270
 271   // this number should be "parametric dependent"
 272   TColgp_SequenceOfPnt points;
 273   TColStd_SequenceOfReal params;
 274   NCollection_Sequence<Standard_Real> aKnotCoeffs;
 275   gp_Pnt p3d;
 276   Standard_Integer iPnt;
 277
 278   // In case of bspline compute parametrization speed on each
 279   // knot interval inside [aFirstParam, aLastParam].
 280   // If quotient = (MaxSpeed / MinSpeed) >= aMaxQuotientCoeff then
 281   // use PerformByProjLib algorithm.
 282   if(!bspl.IsNull())
 283   {
 284     Standard_Real aFirstParam = First; // First parameter of current interval.
 285     Standard_Real aLastParam = Last; // Last parameter of current interval.
 286
 287     // First index computation.
 288     Standard_Integer anIdx = 1;
 289     for(; anIdx <= bspl->NbKnots() && aFirstParam < Last; anIdx++)
 290     {
 291       if(bspl->Knot(anIdx) > First)
 292       {
 293         break;
 294       }
 295     }
 296
 297     GeomAdaptor_Curve aC3DAdaptor(c3d);
 298     Standard_Real aMinParSpeed = Precision::Infinite(); // Minimal parameterization speed.
 299     for(; anIdx <= bspl->NbKnots() && aFirstParam < Last; anIdx++)
 300     {
 301       // Fill current knot interval.
 302       aLastParam = Min(Last, bspl->Knot(anIdx));
 303       Standard_Integer aNbIntPnts = NCONTROL;
 304       // Number of inner points is adapted according to the length of the interval
 305       // to avoid a lot of calculations on small range of parameters.
 306       if (anIdx > 1)
 307       {
 308         const Standard_Real aLenThres = 1.e-2;
 309         const Standard_Real aLenRatio =
 310             (aLastParam - aFirstParam) / (bspl->Knot(anIdx) - bspl->Knot(anIdx - 1));
 311         if (aLenRatio < aLenThres)
 312         {
 313           aNbIntPnts = Standard_Integer(aLenRatio / aLenThres * aNbIntPnts);
 314           if (aNbIntPnts < 2)
 315             aNbIntPnts = 2;
 316         }
 317       }
 318       Standard_Real aStep = (aLastParam - aFirstParam) / (aNbIntPnts - 1);
 319       Standard_Integer anIntIdx;
 320       gp_Pnt p3d1, p3d2;
 321       // Start filling from first point.
 322       aC3DAdaptor.D0(aFirstParam, p3d1);
 323
 324       Standard_Real aLength3d = 0.0;
 325       for(anIntIdx = 1; anIntIdx < aNbIntPnts; anIntIdx++)
 326       {
 327         Standard_Real aParam = aFirstParam + aStep * anIntIdx;
 328         aC3DAdaptor.D0 (aParam, p3d2);
 329         const Standard_Real aDist = p3d2.Distance(p3d1);
 330
 331         aLength3d += aDist;
 332         p3d1 = p3d2;
 333
 334         aMinParSpeed = Min(aMinParSpeed, aDist / aStep);
 335       }
 336       const Standard_Real aCoeff = aLength3d / (aLastParam - aFirstParam);
 337       if (Abs(aCoeff) > gp::Resolution())
 338         aKnotCoeffs.Append(aCoeff);
 339       aFirstParam = aLastParam;
 340     }
 341
 342     Standard_Real anEvenlyCoeff = 0;
 343     if (aKnotCoeffs.Size() > 0)
 344     {
 345       anEvenlyCoeff = *std::max_element(aKnotCoeffs.begin(), aKnotCoeffs.end()) /
 346                       *std::min_element(aKnotCoeffs.begin(), aKnotCoeffs.end());
 347     }
 348
 349     const Standard_Real aMaxQuotientCoeff = 1500.0;
 350     if (anEvenlyCoeff > aMaxQuotientCoeff &&
 351         aMinParSpeed > Precision::Confusion() )
 352     {
 353       PerformByProjLib(c3d, First, Last, c2d);
 354       // PerformByProjLib fail detection:
 355       if (!c2d.IsNull())
 356       {
 357         return Status (ShapeExtend_DONE);
 358       }
 359     }
 360   }
 361
 362   Standard_Real deltaT, t;
 363   deltaT = (Last - First) / (nbPini-1);
 364   nbrPnt = nbPini;
 365   for (iPnt = 1; iPnt <= nbPini; iPnt ++)
 366   {
 367     if      (iPnt == 1)      t = First;
 368     else if (iPnt == nbPini) t = Last;
 369     else                     t = First + (iPnt - 1) * deltaT;
 370
 371     c3d->D0 (t, p3d);
 372     points.Append(p3d);
 373     params.Append(t);
 374   }
 375
 376   //  CALCUL par approximation
 377   TColgp_SequenceOfPnt2d pnt2d;
 378   ApproxPCurve (nbrPnt,c3d,TolFirst,TolLast,
 379                 points,params,pnt2d,c2d); //szv#4:S4163:12Mar99 OK not needed
 380   nbPini = points.Length();
 381   if (!c2d.IsNull()) {
 382     myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE2);
 383     return Standard_True;
 384   }// cas particulier d iso
 385
 386 //  INTERPOLATION du resultat
 387
 388   if ( myBuild ) {
 389     Handle(TColgp_HArray1OfPnt) thePnts = new TColgp_HArray1OfPnt  (1, nbPini);
 390     Handle(TColStd_HArray1OfReal) theParams = new TColStd_HArray1OfReal(1, nbPini);
 391     for (iPnt = 1; iPnt <= nbPini ; iPnt ++) {
 392       thePnts->SetValue(iPnt, points(iPnt));
 393       theParams->SetValue(iPnt, params(iPnt));
 394     }
 395
 396     Handle(Geom_Curve) newc3d = InterpolateCurve3d (nbPini,thePnts,theParams, c3d);
 397     if ( newc3d.IsNull() ) myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_FAIL2);
 398     else {
 399       myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE3);
 400       c3d = newc3d;
 401     }
 402   }
 403
 404   Handle(TColgp_HArray1OfPnt2d) thePnts2d = new TColgp_HArray1OfPnt2d(1, nbPini);
 405   Handle(TColStd_HArray1OfReal) theParams2d = new TColStd_HArray1OfReal(1, nbPini);
 406   for (iPnt = 1; iPnt <= nbPini ; iPnt ++) {
 407     theParams2d->SetValue(iPnt, params(iPnt));
 408     thePnts2d->SetValue(iPnt, pnt2d(iPnt));
 409   }
 410
 411   c2d = InterpolatePCurve (nbPini, thePnts2d, theParams2d, c3d);
 412 //  c2d = ApproximatePCurve (nbPini, thePnts2d, theParams2d, c3d);
 413 //  Faut-il aussi reprendre la C3D ?
 414   myStatus |= ShapeExtend::EncodeStatus (c2d.IsNull() ? ShapeExtend_FAIL1 : ShapeExtend_DONE2);
 415   return Status (ShapeExtend_DONE);
 416 }
 417
 418 //=======================================================================
 419 //function : PerformByProjLib
 420 //purpose  :
 421 //=======================================================================
 422
 423 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformByProjLib(Handle(Geom_Curve)& c3d,
 424                                                                         const Standard_Real First,
 425                                                                         const Standard_Real Last,
 426                                                                         Handle(Geom2d_Curve)& c2d,
 427                                                                         const GeomAbs_Shape /*continuity*/,
 428                                                                         const Standard_Integer /*maxdeg */,
 429                                                                         const Standard_Integer /*nbinterval */)
 430 {
 431   //Standard_Boolean OK = Standard_True; //szv#4:S4163:12Mar99 unused
 432   c2d.Nullify();
 433   if (mySurf.IsNull()) {
 434     myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL1);
 435     return Standard_False;
 436   }
 437
 438   try
 439   {
 440     OCC_CATCH_SIGNALS
 441     Handle(GeomAdaptor_HSurface) GAS = mySurf->Adaptor3d();
 442     Handle(GeomAdaptor_HCurve) GAC = new GeomAdaptor_HCurve (c3d,First,Last);
 443     ProjLib_ProjectedCurve Projector(GAS, GAC);
 444
 445     switch (Projector.GetType())
 446     {
 447     case GeomAbs_Line :
 448       c2d = new Geom2d_Line(Projector.Line());
 449       break;
 450     case GeomAbs_Circle :
 451       c2d = new Geom2d_Circle(Projector.Circle());
 452       break;
 453     case GeomAbs_Ellipse :
 454       c2d = new Geom2d_Ellipse(Projector.Ellipse());
 455       break;
 456     case GeomAbs_Parabola :
 457       c2d = new Geom2d_Parabola(Projector.Parabola());
 458       break;
 459     case GeomAbs_Hyperbola :
 460       c2d = new Geom2d_Hyperbola(Projector.Hyperbola());
 461       break;
 462     case GeomAbs_BSplineCurve :
 463       c2d = Projector.BSpline();
 464       break;
 465     default:
 466       // Not possible, handling added to avoid gcc warning.
 467       break;
 468     }
 469
 470     if(c2d.IsNull())
 471     {
 472       myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL2);
 473       return Standard_False;
 474     }
 475     else
 476     {
 477       myStatus = ShapeExtend::EncodeStatus (ShapeExtend_DONE1);
 478       return Standard_True;
 479     }
 480
 481   }
 482   catch(Standard_Failure const& anException) {
 483 #ifdef OCCT_DEBUG
 484     std::cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::PerformByProjLib(): Exception: ";
 485     anException.Print(std::cout); std::cout << std::endl;
 486 #endif
 487     (void)anException;
 488     myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL3);
 489     c2d.Nullify();
 490   }
 491   return Standard_False;
 492 }
 493
 494 //=======================================================================
 495 //function : ProjectAnalytic
 496 //purpose  :
 497 //=======================================================================
 498
 499  Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::ProjectAnalytic(const Handle(Geom_Curve)& c3d) const
 500 {
 501   Handle(Geom2d_Curve) result;
 502
 503   //:k1 abv 16 Dec 98: limit analytic cases by Plane surfaces only
 504   // This is necessary for K4L since it fails on other surfaces
 505   // when general method GeomProjLib::Curve2d() is used
 506   // Projection is done as in BRep_Tool and BRepCheck_Edge
 507   Handle(Geom_Surface) surf = mySurf->Surface();
 508   Handle(Geom_Plane) Plane = Handle(Geom_Plane)::DownCast ( surf );
 509   if ( Plane.IsNull() ) {
 510     Handle(Geom_RectangularTrimmedSurface) RTS =
 511       Handle(Geom_RectangularTrimmedSurface)::DownCast ( surf );
 512     if ( ! RTS.IsNull() ) Plane = Handle(Geom_Plane)::DownCast ( RTS->BasisSurface() );
 513     else {
 514       Handle(Geom_OffsetSurface) OS =
 515         Handle(Geom_OffsetSurface)::DownCast ( surf );
 516       if ( ! OS.IsNull() )
 517         Plane = Handle(Geom_Plane)::DownCast ( OS->BasisSurface() );
 518     }
 519   }
 520   if ( ! Plane.IsNull() ) {
 521     Handle(Geom_Curve) ProjOnPlane =
 522       GeomProjLib::ProjectOnPlane (c3d, Plane,
 523                                    Plane->Position().Direction(), Standard_True);
 524     Handle(GeomAdaptor_HCurve) HC = new GeomAdaptor_HCurve ( ProjOnPlane );
 525     ProjLib_ProjectedCurve Proj ( mySurf->Adaptor3d(), HC );
 526
 527     result = Geom2dAdaptor::MakeCurve(Proj);
 528     if ( result.IsNull() ) return result;
 529     if ( result->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve)) ) {
 530       Handle(Geom2d_TrimmedCurve) TC = Handle(Geom2d_TrimmedCurve)::DownCast ( result );
 531       result = TC->BasisCurve();
 532     }
 533
 534     return result;
 535   }
 536
 537   return result;
 538 }
 539
 540  //! Fix possible period jump and handle walking period parameter.
 541  static Standard_Boolean fixPeriodictyTroubles(gp_Pnt2d *thePnt, // pointer to gp_Pnt2d[4] beginning
 542                                                Standard_Integer theIdx, // Index of objective coord: 1 ~ X, 2 ~ Y
 543                                                Standard_Real thePeriod, // Period on objective coord
 544                                                Standard_Integer theSavedPoint, // Point number to choose period
 545                                                Standard_Real theSavedParam) // Param from cashe to choose period
 546 {
 547   Standard_Real aSavedParam;
 548   Standard_Integer aSavedPoint;
 549   Standard_Real aMinParam = 0.0, aMaxParam = thePeriod;
 550   if (theSavedPoint < 0) {
 551     // normalize to first period by default
 552     aSavedParam = 0.5 * thePeriod;
 553     aSavedPoint = 0;
 554   }
 555   else {
 556     aSavedParam =  theSavedParam;
 557     aSavedPoint = theSavedPoint;
 558     while (aMinParam > aSavedParam) {
 559       aMinParam -= thePeriod;
 560       aMaxParam -= thePeriod;
 561     }
 562     while (aMaxParam < aSavedParam) {
 563       aMinParam += thePeriod;
 564       aMaxParam += thePeriod;
 565     }
 566   }
 567
 568   Standard_Real aFixIsoParam = aMinParam;
 569   Standard_Boolean isIsoLine = Standard_False;
 570   if (aMaxParam - aSavedParam < Precision::PConfusion() ||
 571       aSavedParam - aMinParam < Precision::PConfusion()) {
 572     aFixIsoParam = aSavedParam;
 573     isIsoLine = Standard_True;
 574   }
 575   // normalize all coordinates to [aMinParam, aMaxParam)
 576   for (Standard_Integer i = 0; i < 4; i++) {
 577     Standard_Real aParam = thePnt[i].Coord(theIdx);
 578     Standard_Real aShift = ShapeAnalysis::AdjustToPeriod(aParam, aMinParam, aMaxParam);
 579     aParam += aShift;
 580     // Walk over period coord -> not walking on another isoline in parameter space.
 581     if (isIsoLine) {
 582       if (aMaxParam - aParam < Precision::PConfusion() || aParam - aMinParam < Precision::PConfusion())
 583         aParam = aFixIsoParam;
 584     }
 585     else {
 586       if (aMaxParam - aParam < Precision::PConfusion())
 587         aParam = aMaxParam;
 588       if (aParam - aMinParam < Precision::PConfusion())
 589           aParam = aMinParam;
 590     }
 591
 592     thePnt[i].SetCoord(theIdx, aParam);
 593   }
 594
 595   // find possible period jump and increasing or decreasing coordinates vector
 596   Standard_Boolean isJump = Standard_False;
 597   Standard_Real aPrevDiff = 0.0;
 598   Standard_Real aSumDiff = 1.0;
 599   for (Standard_Integer i = 0; i < 3; i++) {
 600     Standard_Real aDiff = thePnt[i + 1].Coord(theIdx) - thePnt[i].Coord(theIdx);
 601     if (aDiff < -Precision::PConfusion()) {
 602       aSumDiff *= -1.0;
 603     }
 604     //if first derivative changes its sign then period jump may exists in this place
 605     if (aDiff * aPrevDiff < -Precision::PConfusion()) {
 606       isJump = Standard_True;
 607     }
 608     aPrevDiff = aDiff;
 609   }
 610
 611   if (!isJump)
 612     return Standard_False;
 613
 614   if (aSumDiff > 0) { // decreasing sequence (parameters decrease twice(--) and one period jump(+))
 615     for (Standard_Integer i = aSavedPoint; i > 0; i--)
 616       if (thePnt[i].Coord(theIdx) > thePnt[i - 1].Coord(theIdx)) {
 617         thePnt[i - 1].SetCoord(theIdx, thePnt[i - 1].Coord(theIdx) + thePeriod);
 618       }
 619     for (Standard_Integer i = aSavedPoint; i < 3; i++)
 620       if (thePnt[i].Coord(theIdx) < thePnt[i + 1].Coord(theIdx)) {
 621         thePnt[i + 1].SetCoord(theIdx, thePnt[i + 1].Coord(theIdx) - thePeriod);
 622       }
 623   }
 624   else {// increasing sequence (parameters increase twice(++) and one period jump(-))
 625     for (Standard_Integer i = aSavedPoint; i > 0; i--)
 626       if (thePnt[i].Coord(theIdx) < thePnt[i - 1].Coord(theIdx)) {
 627         thePnt[i - 1].SetCoord(theIdx, thePnt[i - 1].Coord(theIdx) - thePeriod);
 628       }
 629     for (Standard_Integer i = aSavedPoint; i < 3; i++)
 630       if (thePnt[i].Coord(theIdx) > thePnt[i + 1].Coord(theIdx)) {
 631         thePnt[i + 1].SetCoord(theIdx, thePnt[i + 1].Coord(theIdx) + thePeriod);
 632       }
 633   }
 634
 635   // Do not return false, because for nonlinear 2d curves vector of parameters
 636   // may change its first derivative and shifted parameters will be broken for this case.
 637   return Standard_True;
 638  }
 639
 640 //=======================================================================
 641 //function : getLine
 642 //purpose  :
 643 //=======================================================================
 644
 645  Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::getLine(
 646    const TColgp_SequenceOfPnt& thepoints,
 647    const TColStd_SequenceOfReal& theparams,
 648    TColgp_SequenceOfPnt2d& thePnt2ds,
 649    Standard_Real theTol,
 650    Standard_Boolean &isRecompute,
 651    Standard_Boolean &isFromCashe) const
 652  {
 653    Standard_Integer nb =  thepoints.Length();
 654    gp_Pnt aP[4];
 655    aP[0] = thepoints(1);
 656    aP[1] = thepoints(2);
 657    aP[2] = thepoints(nb - 1);
 658    aP[3] = thepoints(nb);
 659    gp_Pnt2d aP2d[4];
 660    Standard_Integer i = 0;
 661
 662    Standard_Real aTol2 = theTol * theTol;
 663    Standard_Boolean isPeriodicU = mySurf->Surface()->IsUPeriodic();
 664    Standard_Boolean isPeriodicV = mySurf->Surface()->IsVPeriodic();
 665
 666    // Workaround:
 667    // Protection against bad "tolerance" shapes.
 668    if (aTol2 > 1.0)
 669    {
 670      theTol = Precision::Confusion();
 671      aTol2 = theTol * theTol;
 672    }
 673    if (aTol2 < Precision::SquareConfusion())
 674      aTol2 = Precision::SquareConfusion();
 675    Standard_Real anOldTol2 = aTol2;
 676    // auxiliary variables to choose period for connection with previous 2dcurve (if exist)
 677    Standard_Integer aSavedPointNum = -1;
 678    gp_Pnt2d aSavedPoint;
 679
 680    // project first and last points
 681    for( ; i < 4; i +=3)
 682    {
 683      Standard_Integer j;
 684      for (j = 0; j < myNbCashe; j++)
 685      {
 686        if ( myCashe3d[j].SquareDistance (aP[i] ) < aTol2)
 687        {
 688          aP2d[i] = mySurf->NextValueOfUV (myCashe2d[j], aP[i], theTol, theTol);
 689          aSavedPointNum = i;
 690          aSavedPoint = myCashe2d[j];
 691          if (i == 0)
 692            isFromCashe = Standard_True;
 693          break;
 694        }
 695      }
 696
 697      if (j >= myNbCashe)
 698        aP2d[i] = mySurf->ValueOfUV(aP[i], theTol);
 699
 700      Standard_Real aDist = mySurf->Gap();
 701      Standard_Real aCurDist = aDist * aDist;
 702      if( aTol2 < aDist * aDist)
 703        aTol2 = aCurDist;
 704    }
 705
 706    if ( isPeriodicU || isPeriodicV )
 707    {
 708      // Compute second and last but one c2d points.
 709      for(i = 1; i < 3; i++)
 710      {
 711        Standard_Integer j;
 712        for (j = 0; j < myNbCashe; j++)
 713        {
 714          if ( myCashe3d[j].SquareDistance (aP[i] ) < aTol2)
 715          {
 716            aP2d[i] = mySurf->NextValueOfUV (myCashe2d[j], aP[i], theTol, theTol);
 717            aSavedPointNum = i;
 718            aSavedPoint = myCashe2d[j];
 719            break;
 720          }
 721        }
 722
 723        if (j >= myNbCashe)
 724          aP2d[i] = mySurf->ValueOfUV(aP[i], theTol);
 725
 726        Standard_Real aDist = mySurf->Gap();
 727        Standard_Real aCurDist = aDist * aDist;
 728        if( aTol2 < aDist * aDist)
 729          aTol2 = aCurDist;
 730      }
 731
 732      if (isPeriodicU)
 733      {
 734        isRecompute = fixPeriodictyTroubles(&aP2d[0], 1 /* X Coord */, mySurf->Surface()->UPeriod(), aSavedPointNum, aSavedPoint.X());
 735      }
 736
 737      if (isPeriodicV)
 738      {
 739        isRecompute = fixPeriodictyTroubles(&aP2d[0], 2 /* Y Coord */, mySurf->Surface()->VPeriod(), aSavedPointNum, aSavedPoint.Y());
 740      }
 741    }
 742
 743    if (isRecompute && mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_SphericalSurface))) {
 744      // Do not try to make line, because in this case may be very special case when 3d curve
 745      // go over the pole of, e.g., sphere, and partly lies along seam
 746      // (see ApproxPCurve() for more information).
 747      return 0;
 748    }
 749
 750    thePnt2ds.SetValue(1, aP2d[0]);
 751    thePnt2ds.SetValue(nb, aP2d[3]);
 752
 753    // Restore old tolerance in 2d space to avoid big gap cases.
 754    aTol2 = anOldTol2;
 755    // Check that straight line in 2d with parameterisation as in 3d will fit
 756    // fit 3d curve at all points.
 757    Standard_Real dPar = theparams(nb) - theparams(1);
 758    if ( Abs(dPar) < Precision::PConfusion() )
 759      return 0;
 760    gp_Vec2d aVec0 (aP2d[0], aP2d[3]);
 761    gp_Vec2d aVec = aVec0 / dPar;
 762    Handle(Geom_Surface) aSurf = mySurf->Surface();
 763    Standard_Boolean isNormalCheck = aSurf->IsCNu(1) && aSurf->IsCNv(1);
 764    if (isNormalCheck) {
 765      for(i = 1; i <= nb; i++)
 766      {
 767        gp_XY aCurPoint = aP2d[0].XY() + aVec.XY() * (theparams(i) - theparams(1));
 768        gp_Pnt aCurP;
 769        gp_Vec aNormalVec, aDu, aDv;
 770        aSurf->D1(aCurPoint.X(), aCurPoint.Y(), aCurP, aDu, aDv);
 771        aNormalVec = aDu ^ aDv;
 772        if (aNormalVec.SquareMagnitude() < Precision::SquareConfusion()) {
 773          isNormalCheck = Standard_False;
 774          break;
 775        }
 776        gp_Lin aNormalLine(aCurP, gp_Dir(aNormalVec));
 777        Standard_Real aDist = aNormalLine.Distance(thepoints(i));
 778        if (aDist > theTol)
 779          return 0;
 780      }
 781    }
 782    if (!isNormalCheck) {
 783      Standard_Real aFirstPointDist = mySurf->Surface()->Value(aP2d[0].X(), aP2d[0].Y()).
 784                                      SquareDistance(thepoints(1));
 785      aTol2 = Max(aTol2, aTol2 * 2 * aFirstPointDist);
 786      for(i = 2; i <  nb; i++)
 787      {
 788        gp_XY aCurPoint = aP2d[0].XY() + aVec.XY() * (theparams(i) - theparams(1));
 789        gp_Pnt aCurP;
 790        aSurf->D0(aCurPoint.X(), aCurPoint.Y(), aCurP);
 791        Standard_Real aDist1 = aCurP.SquareDistance(thepoints(i));
 792
 793        if(Abs (aFirstPointDist - aDist1) > aTol2)
 794          return 0;
 795      }
 796    }
 797
 798    // check if pcurve can be represented by Geom2d_Line (parameterised by length)
 799    Standard_Real aLLength = aVec0.Magnitude();
 800    if ( Abs (aLLength - dPar) <= Precision::PConfusion() )
 801    {
 802      gp_XY aDirL = aVec0.XY() / aLLength;
 803      gp_Pnt2d aPL (aP2d[0].XY() - theparams(1) * aDirL);
 804      return new Geom2d_Line (aPL, gp_Dir2d(aDirL));
 805    }
 806
 807    // create straight bspline
 808    TColgp_Array1OfPnt2d aPoles(1, 2);
 809    aPoles(1) = aP2d[0];
 810    aPoles(2) = aP2d[3];
 811
 812    TColStd_Array1OfReal aKnots(1,2);
 813    aKnots(1) = theparams(1);
 814    aKnots(2) = theparams(theparams.Length());
 815
 816    TColStd_Array1OfInteger aMults(1,2);
 817    aMults(1) = 2;
 818    aMults(2) = 2;
 819    Standard_Integer aDegree = 1;
 820    Handle(Geom2d_BSplineCurve) abspl2d =
 821      new Geom2d_BSplineCurve (aPoles, aKnots, aMults, aDegree);
 822    return abspl2d;
 823  }
 824
 825 //=======================================================================
 826 //function : ApproxPCurve
 827 //purpose  :
 828 //=======================================================================
 829
 830   Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::ApproxPCurve(const Standard_Integer nbrPnt,
 831                                                                       const Handle(Geom_Curve)& c3d,
 832                                                                       const Standard_Real TolFirst,
 833                                                                       const Standard_Real TolLast,
 834                                                                       TColgp_SequenceOfPnt& points,
 835                                                                       TColStd_SequenceOfReal& params,
 836                                                                       TColgp_SequenceOfPnt2d& pnt2d,
 837                                                                       Handle(Geom2d_Curve)& c2d)
 838 {
 839   // for performance, first try to handle typical case when pcurve is straight
 840   Standard_Boolean isRecompute = Standard_False;
 841   Standard_Boolean isFromCasheLine = Standard_False;
 842   for (Standard_Integer iseq = 1; iseq <= nbrPnt; iseq++)
 843   {
 844     gp_Pnt2d aP2d(0.,0.);
 845     pnt2d.Append(aP2d);
 846   }
 847   c2d = getLine(points, params, pnt2d, myPreci, isRecompute, isFromCasheLine);
 848   if(!c2d.IsNull())
 849   {
 850     // fill cashe
 851     Standard_Boolean ChangeCycle = Standard_False;
 852     if(myNbCashe>0 && myCashe3d[0].Distance(points(1)) > myCashe3d[0].Distance(points(nbrPnt)) &&
 853        myCashe3d[0].Distance(points(nbrPnt))<Precision::Confusion())
 854       ChangeCycle = Standard_True;
 855     myNbCashe = 2;
 856     if(ChangeCycle) {
 857       myCashe3d[0] = points(1);
 858       myCashe3d[1] = points(nbrPnt);
 859       myCashe2d[0] = pnt2d(1);
 860       myCashe2d[1] = pnt2d(nbrPnt);
 861     }
 862     else {
 863       myCashe3d[1] = points(1);
 864       myCashe3d[0] = points(nbrPnt);
 865       myCashe2d[1] = pnt2d(1);
 866       myCashe2d[0] = pnt2d(nbrPnt);
 867     }
 868     return Standard_True;
 869   }
 870     Standard_Boolean isDone = Standard_True;
 871   // test if the curve 3d is a boundary of the surface
 872   // (only for Bezier or BSpline surface)
 873
 874   Standard_Boolean isoParam, isoPar2d3d, isoTypeU, p1OnIso, p2OnIso, isoclosed;
 875   gp_Pnt2d valueP1, valueP2;
 876   Handle(Geom_Curve) cIso;
 877   Standard_Real t1, t2;
 878
 879   Handle(Standard_Type) sType = mySurf->Surface()->DynamicType();
 880   Standard_Boolean isAnalytic = Standard_True;
 881   if (sType == STANDARD_TYPE(Geom_BezierSurface) || sType == STANDARD_TYPE(Geom_BSplineSurface)) isAnalytic = Standard_False;
 882   Standard_Real uf, ul, vf, vl;
 883   mySurf->Surface()->Bounds(uf, ul, vf, vl);
 884   isoclosed = Standard_False;
 885   TColStd_Array1OfReal pout(1, nbrPnt);
 886
 887   isoParam = IsAnIsoparametric(nbrPnt, points, params,
 888                                isoTypeU, p1OnIso, valueP1, p2OnIso, valueP2,
 889                                isoPar2d3d, cIso, t1, t2, pout);
 890
 891   // projection of the points on surfaces
 892
 893   gp_Pnt p3d;
 894   gp_Pnt2d p2d;
 895   Standard_Real isoValue=0., isoPar1=0., isoPar2=0., tPar=0., tdeb,tfin;
 896   Standard_Real Cf, Cl, parf, parl; //szv#4:S4163:12Mar99 dist not needed
 897
 898   //  Le calcul part-il dans le bon sens, c-a-d deb et fin dans le bon ordre ?
 899   //  Si uclosed et iso en V, attention isoPar1 ET/OU 2 peut toucher la fermeture
 900   if(isoParam){
 901     if(isoTypeU){
 902       isoValue = valueP1.X();
 903       isoPar1 = valueP1.Y();
 904       isoPar2 = valueP2.Y();
 905       isoclosed = mySurf->IsVClosed(myPreci);//#78 rln 12.03.99 S4135
 906       parf = vf;  parl = vl;
 907     }
 908     else {
 909       isoValue = valueP1.Y();
 910       isoPar1 = valueP1.X();
 911       isoPar2 = valueP2.X();
 912       isoclosed = mySurf->IsUClosed(myPreci);//#78 rln 12.03.99 S4135
 913       parf = uf;  parl = ul;
 914     }
 915     if (!isoPar2d3d && !isAnalytic) {
 916       Cf = cIso->FirstParameter();
 917       Cl = cIso->LastParameter();
 918       if (Precision::IsInfinite(Cf))    Cf = -1000;
 919       if (Precision::IsInfinite(Cl))    Cl = +1000;
 920       //pdn S4030 optimizing and fix isopar case on PRO41323
 921       tdeb = pout(2);
 922       //    dist = ShapeAnalysis_Curve().Project (cIso,points(2),myPreci,pt,tdeb,Cf,Cl);
 923       //  Chacun des par1 ou par2 est-il sur un bord. Attention first/last : recaler
 924       if (isoclosed && (isoPar1 == parf || isoPar1 == parl)) {
 925         if (Abs(tdeb-parf) < Abs(tdeb-parl)) isoPar1 = parf;
 926         else isoPar1 = parl;
 927         if (isoTypeU) valueP1.SetY (isoPar1);
 928         else          valueP1.SetX (isoPar1);
 929       }
 930       if (isoclosed && (isoPar2 == parf || isoPar2 == parl)) {
 931         //pdn S4030 optimizing and fix isopar case on PRO41323
 932         tfin = pout(nbrPnt-1);
 933         //dist =  ShapeAnalysis_Curve().Project (cIso,points(nbrPnt-1),myPreci,pt,tfin,Cf,Cl);
 934         if (Abs(tfin-parf) < Abs(tfin-parl)) isoPar2 = parf;
 935         else isoPar2 = parl;
 936         if (isoTypeU) valueP2.SetY (isoPar2);
 937         else          valueP2.SetX (isoPar2);
 938       }
 939
 940       //  Interversion Par1/Par2 (ne veut que si les 2 sont sur les bords ...)
 941       //  Est-ce encore necessaire apres ce qui vient d etre fait ?
 942
 943       // PTV 05.02.02 fix for translation face from 12_hp_mouse (PARASOLID) face 24008
 944       // if curve is periodic do not change the points
 945       // skl change "if" for pout(nbrPnt-1) 19.11.2003
 946       if (!isoclosed) {
 947         if( (Abs(tdeb-isoPar1)>Abs(tdeb-isoPar2)) &&
 948             (Abs(pout(nbrPnt-1)-isoPar2)>Abs(pout(nbrPnt-1)-isoPar1)) ) {
 949           gp_Pnt2d valueTmp = valueP1;
 950           valueP1 = valueP2;  valueP2 = valueTmp;
 951           if (isoTypeU) {
 952             isoValue = valueP1.X();
 953             isoPar1 = valueP1.Y();
 954             isoPar2 = valueP2.Y();
 955           }
 956           else {
 957             isoValue = valueP1.Y();
 958             isoPar1 = valueP1.X();
 959             isoPar2 = valueP2.X();
 960           }
 961           //  Fin calcul sens de courbe iso
 962         }
 963       } // end of fix check 05.02.02
 964     }
 965   }
 966
 967   //  Si pas isoParam, on a quand meme du p1OnIso/p2OnIso possible ... !!!
 968   //  (utile pour detromper bug de projection). Mais detromper aussi circularite
 969   //else {
 970     //if (p1OnIso) valueP1 =
 971     //BestExtremum (valueP1,points(1),points(2));
 972     //if (p2OnIso) valueP2 =
 973     //BestExtremum (valueP2,points(nbrPnt),points(nbrPnt-1));
 974     //}
 975
 976   Standard_Real gap = myPreci; //:q1
 977   Standard_Boolean ChangeCycle = Standard_False; //skl for OCC3430
 978   // auxiliaruy variables to shift 2dcurve, according to previous
 979   Standard_Boolean isFromCashe = Standard_False;
 980   gp_Pnt2d aSavedPoint;
 981   if( myNbCashe>0 && myCashe3d[0].Distance(points(1))>myCashe3d[0].Distance(points(nbrPnt)) )
 982     //if(myCashe3d[0].Distance(points(nbrPnt))<myPreci)
 983     if(myCashe3d[0].Distance(points(nbrPnt))<Precision::Confusion())
 984       ChangeCycle = Standard_True;
 985   //for( i = 1; i <= nbrPnt; i ++) {
 986   for(Standard_Integer ii=1; ii<=nbrPnt; ii++) {
 987     const Standard_Integer aPntIndex = ChangeCycle ? (nbrPnt - ii + 1) : ii;
 988     p3d = points (aPntIndex);
 989     if (isoParam) {
 990
 991       if (isoPar2d3d) {
 992         if (isoPar2 > isoPar1) tPar = params (aPntIndex);
 993         else                   tPar = t1 + t2 - params(aPntIndex);
 994       } else if (!isAnalytic) {
 995         // projection to iso
 996         if      (aPntIndex == 1)      tPar = isoPar1;
 997         else if (aPntIndex == nbrPnt) tPar = isoPar2;
 998         else {
 999           tPar = pout(aPntIndex);
1000           //:S4030  ShapeAnalysis_Curve().Project (cIso,p3d,myPreci,pt,tPar,Cf,Cl); //szv#4:S4163:12Mar99 `dist=` not needed
1001         }
1002       }
1003
1004       if (!isoPar2d3d && isAnalytic) {
1005         if      (aPntIndex == 1)      p2d = valueP1;
1006         else if (aPntIndex == nbrPnt) p2d = valueP2;
1007         else {
1008           p2d = mySurf->NextValueOfUV(p2d,p3d, myPreci, //%12 pdn 15.02.99 optimizing
1009                                       Precision::Confusion()+1000*gap); //:q1
1010           gap = mySurf->Gap();
1011         }
1012       } else {
1013         if(isoTypeU)  {  p2d.SetX(isoValue);  p2d.SetY(tPar);     }
1014         else          {  p2d.SetX(tPar);      p2d.SetY(isoValue); }
1015       }
1016     }
1017
1018     else {
1019       if     ( (aPntIndex == 1)      && p1OnIso)  p2d = valueP1;
1020       else if( (aPntIndex == nbrPnt) && p2OnIso)  p2d = valueP2;
1021       else  {// general case (not an iso)  mais attention aux singularites !
1022         // first and last points are already computed by getLine()
1023         if (aPntIndex == 1 || aPntIndex == nbrPnt)
1024         {
1025           if (!isRecompute)
1026           {
1027             p2d = pnt2d (aPntIndex);
1028             gap = mySurf->Gap();
1029             if (aPntIndex == 1) {
1030               isFromCashe = isFromCasheLine;
1031               aSavedPoint = p2d;
1032             }
1033             continue;
1034           }
1035           else
1036           {
1037             //:q9 abv 23 Mar 99: use cashe as 1st approach
1038             Standard_Integer j; // svv #1
1039             for (j = 0; j < myNbCashe; ++j)
1040             {
1041               if ( myCashe3d[j].SquareDistance ( p3d ) < myPreci*myPreci )
1042               {
1043                 p2d = mySurf->NextValueOfUV (myCashe2d[j], p3d, myPreci, Precision::Confusion()+gap);
1044                 if (aPntIndex == 1)
1045                 {
1046                   isFromCashe = Standard_True;
1047                   aSavedPoint = myCashe2d[j];
1048                 }
1049                 break;
1050               }
1051             }
1052             if (j >= myNbCashe)
1053             {
1054               p2d = mySurf->ValueOfUV(p3d, myPreci);
1055             }
1056           }
1057         }
1058         else {
1059           p2d = mySurf->NextValueOfUV (p2d, p3d, myPreci, //:S4030: optimizing
1060                                        Precision::Confusion()+1000*gap); //:q1
1061         }
1062         gap = mySurf->Gap();
1063       }
1064     }
1065     pnt2d (aPntIndex) = p2d;
1066     if ( ii > 1 ) {
1067       if(ChangeCycle)
1068         p2d.SetXY ( 2. * p2d.XY() - pnt2d(aPntIndex + 1).XY() );
1069       else
1070         p2d.SetXY ( 2. * p2d.XY() - pnt2d(aPntIndex - 1).XY() );
1071     }
1072   }
1073
1074   //pdn %12 11.02.99 PRO9234 entity 15402
1075   if (!isoPar2d3d) {
1076     mySurf->ProjectDegenerated(nbrPnt,points,pnt2d,myPreci,Standard_True);
1077     mySurf->ProjectDegenerated(nbrPnt,points,pnt2d,myPreci,Standard_False);
1078   }
1079
1080   //Check the extremities of 3d curve for coinciding with singularities of surf
1081   //Standard_Integer NbSing = mySurf->NbSingularities(Precision::Confusion());
1082   gp_Pnt PointFirst = points.First(), PointLast = points.Last();
1083   Standard_Real aTolFirst = (TolFirst == -1)? Precision::Confusion() : TolFirst;
1084   Standard_Real aTolLast  = (TolLast == -1)?  Precision::Confusion() : TolLast;
1085   for (Standard_Integer i = 1; ; i++)
1086   {
1087     Standard_Real aPreci, aFirstPar, aLastPar;
1088     gp_Pnt aP3d;
1089     gp_Pnt2d aFirstP2d, aLastP2d;
1090     Standard_Boolean IsUiso;
1091     if (!mySurf->Singularity(i, aPreci, aP3d, aFirstP2d, aLastP2d, aFirstPar, aLastPar, IsUiso))
1092       break;
1093     if (aPreci <= Precision::Confusion() &&
1094         PointFirst.Distance(aP3d) <= aTolFirst)
1095     {
1096       CorrectExtremity(c3d, params, pnt2d,
1097                        Standard_True, //first point
1098                        aFirstP2d,
1099                        IsUiso);
1100     }
1101     if (aPreci <= Precision::Confusion() &&
1102         PointLast.Distance(aP3d) <= aTolLast)
1103     {
1104       CorrectExtremity(c3d, params, pnt2d,
1105                        Standard_False, //last point
1106                        aFirstP2d,
1107                        IsUiso);
1108     }
1109   }
1110
1111   //  attention aux singularites ... (hors cas iso qui les traite deja)
1112   //  if (!isoParam) {
1113     //    p2d = pnt2d (1);
1114     //    if (mySurf->ProjectDegenerated (points(1),myPreci,pnt2d (2),p2d))
1115     //      pnt2d (1) = p2d;
1116     //    p2d = pnt2d (nbrPnt);
1117     //    if (mySurf->ProjectDegenerated (points(nbrPnt),myPreci,pnt2d (nbrPnt-1),p2d))
1118     //      pnt2d (nbrPnt) = p2d;
1119     //  }
1120
1121   // Si la surface est UCLosed et VClosed, on recadre les points
1122   // algo un peu complique, on retarde l implementation
1123   Standard_Real Up = ul - uf;
1124   Standard_Real Vp = vl - vf;
1125   Standard_Real dist2d;
1126   const Standard_Real TolOnUPeriod = Precision::Confusion() * Up;
1127   const Standard_Real TolOnVPeriod = Precision::Confusion() * Vp;
1128 #ifdef OCCT_DEBUG
1129   if (mySurf->IsUClosed(myPreci) && mySurf->IsVClosed(myPreci)) {//#78 rln 12.03.99 S4135
1130     std::cout << "WARNING : Recadrage incertain sur U & VClosed" << std::endl;
1131   }
1132 #endif
1133   // Si la surface est UCLosed, on recadre les points
1134   if (mySurf->IsUClosed(myPreci)) {//#78 rln 12.03.99 S4135
1135     // Premier point dans le domain [uf, ul]
1136     Standard_Real prevX, firstX = pnt2d (1).X();
1137     if (!isFromCashe) {
1138       // do not shift 2dcurve, if it connects to previous
1139       while (firstX < uf)  {  firstX += Up;   pnt2d (1).SetX(firstX);  }
1140       while (firstX > ul)  {  firstX -= Up;   pnt2d (1).SetX(firstX);  }
1141     }
1142     // shift first point, according to cashe
1143     if (mySurf->Surface()->IsUPeriodic() && isFromCashe) {
1144       Standard_Real aMinParam = uf, aMaxParam = ul;
1145       while (aMinParam > aSavedPoint.X()) {
1146         aMinParam -= Up;
1147         aMaxParam -= Up;
1148       }
1149       while (aMaxParam < aSavedPoint.X()) {
1150         aMinParam += Up;
1151         aMaxParam += Up;
1152       }
1153       Standard_Real aShift = ShapeAnalysis::AdjustToPeriod(firstX, aMinParam, aMaxParam);
1154       firstX += aShift;
1155       pnt2d(1).SetX(firstX);
1156     }
1157     prevX = firstX;
1158
1159     //:97 abv 1 Feb 98: treat case when curve is whole out of surface bounds
1160     Standard_Real minX = firstX, maxX = firstX;
1161     Standard_Boolean ToAdjust = Standard_False;
1162
1163     // On decalle toujours le suivant
1164     for (Standard_Integer aPntIter = 2; aPntIter <= pnt2d.Length(); ++aPntIter)
1165     {
1166       //      dist2d = pnt2d (aPntIter - 1).Distance(pnt2d (aPntIter));
1167       Standard_Real CurX = pnt2d (aPntIter).X();
1168       dist2d = Abs (CurX - prevX);
1169       if (dist2d > ( Up / 2) )
1170       {
1171         InsertAdditionalPointOrAdjust(ToAdjust, 1, Up, TolOnUPeriod,
1172                                       CurX, prevX,
1173                                       c3d,
1174                                       aPntIter,
1175                                       points, params, pnt2d);
1176       }
1177       prevX = CurX;
1178       if ( minX > CurX ) minX = CurX;      //:97
1179       else if ( maxX < CurX ) maxX = CurX; //:97
1180     }
1181
1182     //:97
1183     if (!isFromCashe) {
1184       // do not shift 2dcurve, if it connects to previous
1185       Standard_Real midX = 0.5 * ( minX + maxX );
1186       Standard_Real shiftX=0.;
1187       if ( midX > ul ) shiftX = -Up;
1188       else if ( midX < uf ) shiftX = Up;
1189       if ( shiftX != 0. )
1190         for (Standard_Integer aPntIter = 1; aPntIter <= pnt2d.Length(); ++aPntIter)
1191           pnt2d (aPntIter).SetX ( pnt2d (aPntIter).X() + shiftX );
1192       }
1193   }
1194   // Si la surface est VCLosed, on recadre les points
1195   // Same code as UClosed : optimisation souhaitable !!
1196   // CKY : d abord un code IDENTIQUE A UClosed; PUIS le special Seam ...
1197   // Si la surface est UCLosed, on recadre les points
1198   //
1199   //#69 rln 01.03.99 S4135 bm2_sd_t4-A.stp entity 30
1200   //#78 rln 12.03.99 S4135
1201   if (mySurf->IsVClosed(myPreci) || mySurf->Surface()->IsKind (STANDARD_TYPE (Geom_SphericalSurface))) {
1202     // Premier point dans le domain [vf, vl]
1203     Standard_Real prevY, firstY = pnt2d (1).Y();
1204     if (!isFromCashe) {
1205       // do not shift 2dcurve, if it connects to previous
1206       while (firstY < vf)  {  firstY += Vp;  pnt2d (1).SetY(firstY);  }
1207       while (firstY > vl)  {  firstY -= Vp;  pnt2d (1).SetY(firstY);  }
1208     }
1209     // shift first point, according to cashe
1210     if (mySurf->Surface()->IsVPeriodic() && isFromCashe) {
1211       Standard_Real aMinParam = vf, aMaxParam = vl;
1212       while (aMinParam > aSavedPoint.Y()) {
1213         aMinParam -= Vp;
1214         aMaxParam -= Vp;
1215       }
1216       while (aMaxParam < aSavedPoint.Y()) {
1217         aMinParam += Vp;
1218         aMaxParam += Vp;
1219       }
1220       Standard_Real aShift = ShapeAnalysis::AdjustToPeriod(firstY, aMinParam, aMaxParam);
1221       firstY += aShift;
1222       pnt2d(1).SetY(firstY);
1223     }
1224     prevY = firstY;
1225
1226     //:97 abv 1 Feb 98: treat case when curve is whole out of surface bounds
1227     Standard_Real minY = firstY, maxY = firstY;
1228     Standard_Boolean ToAdjust = Standard_False;
1229
1230     // On decalle toujours le suivant
1231     for (Standard_Integer aPntIter = 2; aPntIter <= pnt2d.Length(); ++aPntIter)
1232     {
1233       //      dist2d = pnt2d (i-1).Distance(pnt2d (i));
1234       Standard_Real CurY = pnt2d (aPntIter).Y();
1235       dist2d = Abs (CurY - prevY);
1236       if (dist2d > ( Vp / 2) )
1237       {
1238         InsertAdditionalPointOrAdjust(ToAdjust, 2, Vp, TolOnVPeriod,
1239                                       CurY, prevY,
1240                                       c3d,
1241                                       aPntIter,
1242                                       points, params, pnt2d);
1243       }
1244       prevY = CurY;
1245       if ( minY > CurY ) minY = CurY;      //:97
1246       else if ( maxY < CurY ) maxY = CurY; //:97
1247     }
1248
1249     //:97
1250     if (!isFromCashe) {
1251       // do not shift 2dcurve, if it connects to previous
1252       Standard_Real midY = 0.5 * ( minY + maxY );
1253       Standard_Real shiftY=0.;
1254       if ( midY > vl ) shiftY = -Vp;
1255       else if ( midY < vf ) shiftY = Vp;
1256       if ( shiftY != 0. )
1257         for (Standard_Integer aPntIter = 1; aPntIter <= pnt2d.Length(); ++aPntIter)
1258           pnt2d(aPntIter).SetY ( pnt2d(aPntIter).Y() + shiftY );
1259       }
1260   }
1261
1262   //#69 rln 01.03.99 S4135 bm2_sd_t4-A.stp entity 30
1263   //#78 rln 12.03.99 S4135
1264   if (mySurf->IsVClosed(myPreci) || mySurf->Surface()->IsKind (STANDARD_TYPE (Geom_SphericalSurface))) {
1265     for (Standard_Integer aPntIter = 2; aPntIter <= pnt2d.Length(); ++aPntIter) {
1266       //#1 rln 11/02/98 ca_exhaust.stp entity #9869 dist2d = pnt2d (i-1).Distance(pnt2d (i));
1267       dist2d = Abs (pnt2d (aPntIter).Y() - pnt2d (aPntIter - 1).Y());
1268       if (dist2d > ( Vp / 2) ) {
1269         // ATTENTION : il faut regarder ou le decalage se fait.
1270         // si plusieurs points sont decalles, il faut plusieurs passes
1271         // pour obtenir un resultat correct.
1272         // NOT YET IMPLEMENTED
1273
1274         // one of those point is incorrectly placed
1275         // i.e on the wrong side of the "seam"
1276         // on prend le point le plus pres des bords vf ou vl
1277         Standard_Boolean prevOnFirst = Standard_False;
1278         Standard_Boolean prevOnLast  = Standard_False;
1279         Standard_Boolean currOnFirst = Standard_False;
1280         Standard_Boolean currOnLast  = Standard_False;
1281
1282         //  .X ?  plutot .Y ,  non ?
1283         Standard_Real distPrevVF = Abs(pnt2d (aPntIter - 1).Y() - vf);
1284         Standard_Real distPrevVL = Abs(pnt2d (aPntIter - 1).Y() - vl);
1285         Standard_Real distCurrVF = Abs(pnt2d (aPntIter).Y() - vf);
1286         Standard_Real distCurrVL = Abs(pnt2d (aPntIter).Y() - vl);
1287
1288         Standard_Real theMin = distPrevVF;
1289         prevOnFirst = Standard_True;
1290         if (distPrevVL < theMin) {
1291           theMin = distPrevVL;
1292           prevOnFirst = Standard_False;
1293           prevOnLast  = Standard_True;
1294         }
1295         if (distCurrVF < theMin) {
1296           theMin = distCurrVF;
1297           prevOnFirst = Standard_False;
1298           prevOnLast  = Standard_False;
1299           currOnFirst = Standard_True;
1300         }
1301         if (distCurrVL < theMin) {
1302           theMin = distCurrVL;
1303           prevOnFirst = Standard_False;
1304           prevOnLast  = Standard_False;
1305           currOnFirst = Standard_False;
1306           currOnLast  = Standard_True;
1307         }
1308         //  Modifs RLN/Nijni  3-DEC-1997
1309         if (prevOnFirst) {
1310           // on decalle le point (aPntIter - 1) en V Last
1311           gp_Pnt2d newPrev(pnt2d (aPntIter - 1).X(), vf); // instead of  vl RLN/Nijni
1312           pnt2d (aPntIter - 1) = newPrev;
1313         }
1314         else if (prevOnLast) {
1315           // on decalle le point (aPntIter - 1) en V first
1316           gp_Pnt2d newPrev(pnt2d (aPntIter - 1).X(), vl); // instead of  vf RLN/Nijni
1317           pnt2d (aPntIter - 1) = newPrev;
1318         }
1319         else if (currOnFirst) {
1320           // on decalle le point (aPntIter) en V Last
1321           gp_Pnt2d newCurr(pnt2d (aPntIter).X(),vf);  // instead of vl  RLN/Nijni
1322           pnt2d (aPntIter) = newCurr;
1323         }
1324         else if (currOnLast) {
1325           // on decalle le point (aPntIter) en V First
1326           gp_Pnt2d newCurr(pnt2d (aPntIter).X(), vl); // instead of vf  RLN/Nijni
1327           pnt2d (aPntIter) = newCurr;
1328         }
1329         // on verifie
1330 #ifdef OCCT_DEBUG
1331         dist2d = pnt2d (aPntIter - 1).Distance(pnt2d (aPntIter));
1332         if (dist2d > ( Vp / 2) ) {
1333           std::cout << "Echec dans le recadrage" << std::endl;
1334         }
1335 #endif
1336       }
1337     }
1338   }
1339
1340   //:c0 abv 20 Feb 98: treat very special case when 3d curve
1341   // go over the pole of, e.g., sphere, and partly lies along seam.
1342   // 2d representation of such a curve should consist of 3 parts - one on
1343   // regular part of surface (interior), one part along degenerated boundary
1344   // and one along seam.
1345   // Since it cannot be adjusted later by arranging pcurves (curve is single),
1346   // to fix it it is nesessary to have a possibility of adjusting seam
1347   // part of such curve either to left or right boundary of surface.
1348   // Test is performed only if flag AdjustOverDegen is not -1.
1349   // If AdjustOverDegen is True, seam part of curve is adjusted to
1350   // the left, and if False - to the right parametric boundary
1351   // If treated case is detected, flag DONE4 is set to status
1352   // NOTE: currently, precision is Precision::PConfusion() since it
1353   // is enough on encountered example
1354   // (ug_turbine-A.stp from ProSTEP Benchmark #3, entities ##2470 & 5680)
1355   // (r1001_ac.stp from Test Rally #10, face #35027 and others)
1356   if ( myAdjustOverDegen != -1 ) {
1357     if ( mySurf->IsUClosed(myPreci) ) {//#78 rln 12.03.99 S4135
1358       mySurf->IsDegenerated ( gp_Pnt(0,0,0), myPreci );  // pour calculer les dgnr
1359       if ( mySurf->NbSingularities(myPreci) > 0 ) { //rln S4135
1360         // 1st, find gap point (degenerated pole)
1361         Standard_Real PrevX=0.;
1362         Standard_Integer OnBound=0, PrevOnBound=0;
1363         Standard_Integer ind; // svv #1
1364         Standard_Boolean start = Standard_True;
1365         for ( ind=1; ind <= pnt2d.Length(); ind++ ) {
1366           Standard_Real CurX = pnt2d(ind).X();
1367           // abv 16 Mar 00: trj3_s1-ug.stp #697: ignore points in singularity
1368           if ( mySurf->IsDegenerated ( points(ind), Precision::Confusion() ) )
1369             continue;
1370           OnBound = ( Abs ( Abs ( CurX - 0.5 * ( ul + uf ) ) - Up/2 ) <=
1371                      Precision::PConfusion() );
1372           if ( ! start &&  Abs ( Abs ( CurX - PrevX ) - Up/2 ) <= 0.01*Up )
1373             break;
1374           start = Standard_False;
1375           PrevX = CurX;
1376           PrevOnBound = OnBound;
1377         }
1378         // if found, adjust seam part
1379         if ( ind <= pnt2d.Length() ) {
1380           PrevX = ( myAdjustOverDegen ? uf : ul );
1381           Standard_Real dU = Up/2 + Precision::PConfusion();
1382           if ( PrevOnBound ) {
1383             pnt2d(ind-1).SetX ( PrevX );
1384             for ( Standard_Integer j=ind-2; j >0; j-- ) {
1385               Standard_Real CurX = pnt2d(j).X();
1386               while ( CurX < PrevX - dU ) pnt2d(j).SetX ( CurX += Up );
1387               while ( CurX > PrevX + dU ) pnt2d(j).SetX ( CurX -= Up );
1388             }
1389           }
1390           else if ( OnBound ) {
1391             pnt2d(ind).SetX ( PrevX );
1392             for ( Standard_Integer j=ind+1; j <= pnt2d.Length(); j++ ) {
1393               Standard_Real CurX = pnt2d(j).X();
1394               while ( CurX < PrevX - dU ) pnt2d(j).SetX ( CurX += Up );
1395               while ( CurX > PrevX + dU ) pnt2d(j).SetX ( CurX -= Up );
1396             }
1397           }
1398           myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE4);
1399         }
1400       }
1401     }
1402     else if ( mySurf->IsVClosed(myPreci) ) {//#78 rln 12.03.99 S4135
1403       mySurf->IsDegenerated ( gp_Pnt(0,0,0), myPreci );  // pour calculer les dgnr
1404       if ( mySurf->NbSingularities(myPreci) > 0 ) { //rln S4135
1405         // 1st, find gap point (degenerated pole)
1406         Standard_Real PrevY=0.;
1407         Standard_Integer OnBound=0, PrevOnBound=0;
1408         Standard_Integer ind; // svv #1
1409         Standard_Boolean start = Standard_True;
1410         for ( ind=1; ind <= pnt2d.Length(); ind++ ) {
1411           Standard_Real CurY = pnt2d(ind).Y();
1412           // abv 16 Mar 00: trj3_s1-ug.stp #697: ignore points in singularity
1413           if ( mySurf->IsDegenerated ( points(ind), Precision::Confusion() ) )
1414             continue;
1415           OnBound = ( Abs ( Abs ( CurY - 0.5 * ( vl + vf ) ) - Vp/2 ) <=
1416                      Precision::PConfusion() );
1417           if ( ! start &&  Abs ( Abs ( CurY - PrevY ) - Vp/2 ) <= 0.01*Vp )
1418             break;
1419           start = Standard_False;
1420           PrevY = CurY;
1421           PrevOnBound = OnBound;
1422         }
1423         // if found, adjust seam part
1424         if ( ind <= pnt2d.Length() ) {
1425           PrevY = ( myAdjustOverDegen ? vf : vl );
1426           Standard_Real dV = Vp/2 + Precision::PConfusion();
1427           if ( PrevOnBound ) {
1428             pnt2d(ind-1).SetY ( PrevY );
1429             for ( Standard_Integer j=ind-2; j >0; j-- ) {
1430               Standard_Real CurY = pnt2d(j).Y();
1431               while ( CurY < PrevY - dV ) pnt2d(j).SetY ( CurY += Vp );
1432               while ( CurY > PrevY + dV ) pnt2d(j).SetY ( CurY -= Vp );
1433             }
1434           }
1435           else if ( OnBound ) {
1436             pnt2d(ind).SetY ( PrevY );
1437             for ( Standard_Integer j=ind+1; j <= pnt2d.Length(); j++ ) {
1438               Standard_Real CurY = pnt2d(j).Y();
1439               while ( CurY < PrevY - dV ) pnt2d(j).SetY ( CurY += Vp );
1440               while ( CurY > PrevY + dV ) pnt2d(j).SetY ( CurY -= Vp );
1441             }
1442           }
1443           myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE4);
1444         }
1445       }
1446     }
1447   }
1448
1449   //:q9: fill cashe
1450   myNbCashe = 2;
1451   if(ChangeCycle) {  // msv 10.08.04: avoid using of uninitialised field
1452   //if(myCashe3d[0].Distance(points(1))>Precision::Confusion() &&
1453   //   myCashe3d[1].Distance(points(1))>Precision::Confusion()) {
1454     myCashe3d[0] = points(1);
1455     myCashe3d[1] = points.Last();
1456     myCashe2d[0] = pnt2d(1);
1457     myCashe2d[1] = pnt2d.Last();
1458   }
1459   else {
1460     myCashe3d[1] = points(1);
1461     myCashe3d[0] = points.Last();
1462     myCashe2d[1] = pnt2d(1);
1463     myCashe2d[0] = pnt2d.Last();
1464   }
1465   return isDone;
1466 }
1467
1468 //=======================================================================
1469 //function : ApproximatePCurve
1470 //purpose  :
1471 //=======================================================================
1472
1473 Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(const Standard_Integer /*nbrPnt*/,
1474                                                                              Handle(TColgp_HArray1OfPnt2d)& points2d,
1475                                                                              Handle(TColStd_HArray1OfReal)& params,
1476                                                                              const Handle(Geom_Curve)& /*orig*/) const
1477 {
1478 //  Standard_Real resol = Min(mySurf->Adaptor3d()->VResolution(myPreci), mySurf->Adaptor3d()->UResolution(myPreci));
1479   Standard_Real theTolerance2d = myPreci; // (100*nbrPnt);//resol;
1480   Handle(Geom2d_Curve) C2d;
1481   try {
1482     OCC_CATCH_SIGNALS
1483     CheckPoints2d (points2d, params, theTolerance2d);
1484     Standard_Integer numberPnt = points2d->Length();
1485
1486     TColgp_Array1OfPnt points3d(1,numberPnt);
1487     gp_Pnt2d pnt2d;
1488     gp_Pnt pnt;
1489     Standard_Integer i; // svv #1
1490     for( i = 1; i <= numberPnt; i++) {
1491       pnt2d = points2d->Value(i);
1492       pnt.SetCoord(pnt2d.X(),pnt2d.Y(),0);
1493       points3d(i) = pnt;
1494     }
1495
1496     GeomAPI_PointsToBSpline appr(points3d, params->Array1(), 1, 10, GeomAbs_C1, theTolerance2d);
1497     Handle(Geom_BSplineCurve) crv3d = appr.Curve();
1498     Standard_Integer NbPoles = crv3d->NbPoles();
1499     TColgp_Array1OfPnt poles3d (1, NbPoles);
1500     TColgp_Array1OfPnt2d poles2d (1, NbPoles);
1501     crv3d->Poles(poles3d);
1502     for( i = 1; i <= NbPoles; i++) {
1503       pnt2d.SetCoord(poles3d(i).X(),poles3d(i).Y());
1504       poles2d(i) = pnt2d;
1505     }
1506     TColStd_Array1OfReal weights (1,NbPoles);
1507     TColStd_Array1OfInteger multiplicities (1,crv3d->NbKnots());
1508     TColStd_Array1OfReal knots(1,crv3d->NbKnots());
1509     crv3d->Knots(knots);
1510     crv3d->Weights(weights);
1511     crv3d->Multiplicities(multiplicities);
1512     C2d = new Geom2d_BSplineCurve  ( poles2d, weights, knots, multiplicities, crv3d->Degree(), crv3d->IsPeriodic());
1513     return C2d;
1514   }
1515   catch(Standard_Failure const& anException) {
1516 #ifdef OCCT_DEBUG
1517  //:s5
1518     //    debug ...
1519     Standard_Integer nbp = params->Length();
1520     Standard_Integer nb2 = points2d->Length();
1521     std::cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(): Exception: ";
1522     anException.Print(std::cout);
1523     std::cout<<"Pb Geom2dAPI_Approximate, tol2d="<<theTolerance2d<<" NbParams="<<nbp<<" NbPnts="<<nb2<<std::endl;
1524 //     if (nb2 > nbp) nb2 = nbp;
1525 //     Standard_Real rbp,rb2; rbp = nbp; rb2 = nb2;
1526 //     //    dbl.AddString ("NbP2d/NbParams puis  X Y Param -> mini");
1527 //     dbl.AddReals (rb2,rbp);
1528 //     for (Standard_Integer i = 1; i <= nb2; i ++) {
1529 //       gp_XYZ quoi (points2d->Value(i).X(),points2d->Value(i).Y(),params->Value(i) );
1530 //       dbl.AddXYZ (quoi);
1531 //     }
1532 #endif
1533      (void)anException;
1534      C2d.Nullify();
1535   }
1536   return C2d;
1537 }
1538
1539 //=======================================================================
1540 //function : InterpolatePCurve
1541 //purpose  :
1542 //=======================================================================
1543
1544  Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolatePCurve(const Standard_Integer nbrPnt,
1545                                                                               Handle(TColgp_HArray1OfPnt2d)& points2d,
1546                                                                               Handle(TColStd_HArray1OfReal)& params,
1547                                                                               const Handle(Geom_Curve)& /*orig*/) const
1548 {
1549   Handle(Geom2d_Curve) C2d;    // NULL si echec
1550   Standard_Real theTolerance2d = myPreci / (100 * nbrPnt);
1551   try {
1552     OCC_CATCH_SIGNALS
1553     // on verifie d abord s il n y a pas de points confondus
1554     // si besoin on retouche les valeurs ...
1555     CheckPoints2d (points2d, params, theTolerance2d);
1556     Geom2dAPI_Interpolate myInterPol2d (points2d, params,
1557                                         Standard_False, theTolerance2d);
1558     myInterPol2d.Perform();
1559     if (myInterPol2d.IsDone()) C2d = myInterPol2d.Curve();
1560   }
1561   catch(Standard_Failure const& anException) {
1562 #ifdef OCCT_DEBUG
1563 //:s5
1564 // //    debug ...
1565     Standard_Integer nbp = params->Length();
1566     Standard_Integer nb2 = points2d->Length();
1567     std::cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::InterpolatePCurve(): Exception: ";
1568     anException.Print(std::cout);
1569     std::cout<<"Pb Geom2dAPI_Interpolate, tol2d="<<theTolerance2d<<" NbParams="<<nbp<<" NbPnts="<<nb2<<std::endl;
1570 //     if (nb2 > nbp) nb2 = nbp;
1571 //     Standard_Real rbp,rb2; rbp = nbp; rb2 = nb2;
1572 // //    dbl.AddString ("NbP2d/NbParams puis  X Y Param -> mini");
1573 //     dbl.AddReals (rb2,rbp);
1574 //     for (Standard_Integer i = 1; i <= nb2; i ++) {
1575 //       gp_XYZ quoi (points2d->Value(i).X(),points2d->Value(i).Y(),params->Value(i) );
1576 //       dbl.AddXYZ (quoi);
1577 //     }
1578 #endif
1579     (void)anException;
1580     C2d.Nullify();
1581   }
1582   return C2d;
1583 }
1584
1585 //=======================================================================
1586 //function : InterpolateCurve3d
1587 //purpose  :
1588 //=======================================================================
1589
1590 Handle(Geom_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(const Standard_Integer,
1591                                                                             Handle(TColgp_HArray1OfPnt)& points,
1592                                                                             Handle(TColStd_HArray1OfReal)& params,
1593                                                                             const Handle(Geom_Curve)& /*orig*/) const
1594 {
1595   Handle(Geom_Curve) C3d;    // NULL si echec
1596   try {
1597     OCC_CATCH_SIGNALS
1598     Standard_Real Tol = myPreci;
1599     CheckPoints(points, params, Tol);
1600     GeomAPI_Interpolate myInterPol(points, params, Standard_False, Tol);
1601     myInterPol.Perform();
1602     if (myInterPol.IsDone()) C3d = myInterPol.Curve();
1603   }
1604   catch(Standard_Failure const& anException) {
1605 #ifdef OCCT_DEBUG
1606  //:s5
1607     std::cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(): Exception: ";
1608     anException.Print(std::cout); std::cout << std::endl;
1609 #endif
1610     (void)anException;
1611     C3d.Nullify();
1612   }
1613   return C3d;
1614 }
1615
1616 //============================================================================================
1617 //function : CorrectExtremity
1618 //purpose  : corrects first or last 2d point of future curve
1619 //           in the case when it coincids with a singularity of surface
1620 //============================================================================================
1621
1622 void ShapeConstruct_ProjectCurveOnSurface::CorrectExtremity(const Handle(Geom_Curve)& theC3d,
1623                                                             const TColStd_SequenceOfReal& theParams,
1624                                                             TColgp_SequenceOfPnt2d& thePnt2d,
1625                                                             const Standard_Boolean theIsFirstPoint,
1626                                                             const gp_Pnt2d& thePointOnIsoLine,
1627                                                             const Standard_Boolean theIsUiso)
1628 {
1629   Standard_Integer NbPnt = thePnt2d.Length();
1630   Standard_Integer IndCoord = (theIsUiso)? 2 : 1;
1631   Standard_Real SingularityCoord = thePointOnIsoLine.Coord(3-IndCoord);
1632   gp_Pnt2d EndPoint = (theIsFirstPoint)? thePnt2d(1) : thePnt2d(NbPnt);
1633   Standard_Real FinishCoord = EndPoint.Coord(3-IndCoord); //the constant coord of isoline
1634
1635   gp_Dir2d aDir = (theIsUiso)? gp::DY2d() : gp::DX2d();
1636   gp_Lin2d anIsoLine(EndPoint, aDir);
1637   IntRes2d_Domain Dom1, Dom2;
1638
1639   Standard_Boolean IsPeriodic = (theIsUiso)?
1640     mySurf->Surface()->IsVPeriodic() : mySurf->Surface()->IsUPeriodic();
1641
1642   gp_Pnt2d FirstPointOfLine, SecondPointOfLine;
1643   Standard_Real FinishParam, FirstParam, SecondParam;
1644
1645   if (theIsFirstPoint)
1646   {
1647     FirstPointOfLine = thePnt2d(3);
1648     SecondPointOfLine = thePnt2d(2);
1649     FinishParam = theParams(1);
1650     FirstParam = theParams(3);
1651     SecondParam = theParams(2);
1652   }
1653   else //last point
1654   {
1655     FirstPointOfLine = thePnt2d(NbPnt-2);
1656     SecondPointOfLine = thePnt2d(NbPnt-1);
1657     FinishParam = theParams(NbPnt);
1658     FirstParam = theParams(NbPnt-2);
1659     SecondParam = theParams(NbPnt-1);
1660   }
1661
1662   if (SingularityCoord > FinishCoord &&
1663       SecondPointOfLine.Coord(3-IndCoord) > FinishCoord)
1664     return; //the curve passes through the singularity, do nothing
1665   if (SingularityCoord < FinishCoord &&
1666       SecondPointOfLine.Coord(3-IndCoord) < FinishCoord)
1667     return; //the curve passes through the singularity, do nothing
1668   //Check correctness of <EndPoint>
1669   {
1670     const Standard_Real aPrevDist = Abs(SecondPointOfLine.Coord(IndCoord) - FirstPointOfLine.Coord(IndCoord));
1671     const Standard_Real aCurDist  = Abs(EndPoint.Coord(IndCoord) - SecondPointOfLine.Coord(IndCoord));
1672     if (aCurDist <= 2 * aPrevDist)
1673       return;
1674   }
1675
1676   gp_Pnt2d FinishPoint = (theIsUiso)? gp_Pnt2d(FinishCoord, SecondPointOfLine.Y()) :
1677     gp_Pnt2d(SecondPointOfLine.X(), FinishCoord); //first approximation of <FinishPoint>
1678
1679   for (;;)
1680   {
1681     if (Abs(SecondPointOfLine.Coord(3-IndCoord) - FinishCoord) <= 2*Precision::PConfusion())
1682       break;
1683
1684     gp_Vec2d aVec(FirstPointOfLine, SecondPointOfLine);
1685     Standard_Real aSqMagnitude = aVec.SquareMagnitude();
1686     if (aSqMagnitude <= 1.e-32)
1687       break;
1688     aDir.SetCoord(aVec.X(), aVec.Y());
1689
1690     gp_Lin2d aLine(FirstPointOfLine, aDir);
1691     IntCurve_IntConicConic Intersector(anIsoLine, Dom1,
1692                                        aLine,     Dom2,
1693                                        1.e-10, 1.e-10);
1694     if (Intersector.IsDone() && !Intersector.IsEmpty())
1695     {
1696       IntRes2d_IntersectionPoint IntPoint = Intersector.Point(1);
1697       FinishPoint = IntPoint.Value();
1698     }
1699     else
1700       FinishPoint = (theIsUiso)? gp_Pnt2d(FinishCoord, SecondPointOfLine.Y()) :
1701         gp_Pnt2d(SecondPointOfLine.X(), FinishCoord);
1702
1703     gp_Pnt2d PrevPoint = FirstPointOfLine;
1704     FirstPointOfLine = SecondPointOfLine;
1705     FirstParam = SecondParam;
1706     SecondParam = (FirstParam + FinishParam)/2;
1707     if (Abs(SecondParam - FirstParam) <= 2*Precision::PConfusion())
1708       break;
1709     gp_Pnt aP3d;
1710     theC3d->D0(SecondParam, aP3d);
1711     SecondPointOfLine = mySurf->NextValueOfUV(FirstPointOfLine, aP3d,
1712                                               myPreci, Precision::Confusion());
1713     if (IsPeriodic)
1714       AdjustSecondPointToFirstPoint(FirstPointOfLine, SecondPointOfLine, mySurf->Surface());
1715
1716     //Check <SecondPointOfLine> to be enough close to <FirstPointOfLine>
1717     //because when a projected point is too close to singularity,
1718     //the non-constant coordinate becomes random.
1719     const Standard_Real aPrevDist = Abs(FirstPointOfLine.Coord(IndCoord) - PrevPoint.Coord(IndCoord));
1720     const Standard_Real aCurDist  = Abs(SecondPointOfLine.Coord(IndCoord) - FirstPointOfLine.Coord(IndCoord));
1721     if (aCurDist > 2 * aPrevDist)
1722       break;
1723   }
1724
1725   if (theIsFirstPoint)
1726     thePnt2d(1) = FinishPoint;
1727   else
1728     thePnt2d(NbPnt) = FinishPoint;
1729 }
1730
1731 //============================================================================================
1732 //function : InsertAdditionalPointOrAdjust
1733 //purpose  : If the current point is too far from the previous point
1734 //           (more than half-period of surface), it can happen in two cases:
1735 //           1. Real current step on corresponding coordinate is small, all we need is adjust;
1736 //           2. Current step on corresponding coordinate is really bigger than half-period of
1737 //           surface in this parametric direction, so we must add additional point to exclude
1738 //           such big intervals between points in 2d space.
1739 //============================================================================================
1740
1741 void ShapeConstruct_ProjectCurveOnSurface::
1742 InsertAdditionalPointOrAdjust(Standard_Boolean& ToAdjust,
1743                               const Standard_Integer theIndCoord,
1744                               const Standard_Real Period,
1745                               const Standard_Real TolOnPeriod,
1746                               Standard_Real& CurCoord,
1747                               const Standard_Real prevCoord,
1748                               const Handle(Geom_Curve)& c3d,
1749                               Standard_Integer& theIndex,
1750                               TColgp_SequenceOfPnt& points,
1751                               TColStd_SequenceOfReal& params,
1752                               TColgp_SequenceOfPnt2d& pnt2d)
1753 {
1754   Standard_Real CorrectedCurCoord = ElCLib::InPeriod(CurCoord,
1755                                                      prevCoord - Period/2,
1756                                                      prevCoord + Period/2);
1757   if (!ToAdjust)
1758   {
1759     Standard_Real CurPar = params(theIndex);
1760     Standard_Real PrevPar = params(theIndex-1);
1761     Standard_Real MidPar = (PrevPar + CurPar)/2;
1762     gp_Pnt MidP3d;
1763     c3d->D0(MidPar, MidP3d);
1764     gp_Pnt2d MidP2d = mySurf->ValueOfUV(MidP3d, myPreci);
1765     Standard_Real MidCoord = MidP2d.Coord(theIndCoord);
1766     MidCoord = ElCLib::InPeriod(MidCoord, prevCoord - Period/2, prevCoord + Period/2);
1767     Standard_Real FirstCoord = prevCoord, LastCoord = CorrectedCurCoord;
1768     if (LastCoord < FirstCoord)
1769     {Standard_Real tmp = FirstCoord; FirstCoord = LastCoord; LastCoord = tmp;}
1770     if (LastCoord - FirstCoord <= TolOnPeriod)
1771       ToAdjust = Standard_True;
1772     else if (FirstCoord <= MidCoord && MidCoord <= LastCoord)
1773       ToAdjust = Standard_True;
1774     else //add mid point
1775     {
1776       //Standard_Real RefU = prevX;
1777       Standard_Boolean Success = Standard_True;
1778       Standard_Real FirstT = PrevPar; //params(i-1)
1779       Standard_Real LastT  = CurPar;  //params(i)
1780       MidCoord = MidP2d.Coord(theIndCoord);
1781       while (Abs(MidCoord - prevCoord) >= Period/2 - TolOnPeriod ||
1782              Abs(CurCoord - MidCoord) >= Period/2 - TolOnPeriod)
1783       {
1784         if (MidPar - FirstT <= Precision::PConfusion() ||
1785             LastT - MidPar <= Precision::PConfusion())
1786         {
1787           Success = Standard_False;
1788           break; //wrong choice
1789         }
1790         if (Abs(MidCoord - prevCoord) >= Period/2 - TolOnPeriod)
1791           LastT = (FirstT + LastT)/2;
1792         else
1793           FirstT = (FirstT + LastT)/2;
1794         MidPar = (FirstT + LastT)/2;
1795         c3d->D0(MidPar, MidP3d);
1796         MidP2d = mySurf->ValueOfUV(MidP3d, myPreci);
1797         MidCoord = MidP2d.Coord(theIndCoord);
1798       }
1799       if (Success)
1800       {
1801         points.InsertBefore(theIndex, MidP3d);
1802         params.InsertBefore(theIndex, MidPar);
1803         pnt2d.InsertBefore(theIndex, MidP2d);
1804         theIndex++;
1805       }
1806       else
1807         ToAdjust = Standard_True;
1808     } //add mid point
1809   } //if (!ToAdjust)
1810   if (ToAdjust)
1811   {
1812     CurCoord = CorrectedCurCoord;
1813     pnt2d(theIndex).SetCoord (theIndCoord, CurCoord);
1814   }
1815 }
1816
1817 //=======================================================================
1818 //function : CheckPoints
1819 //purpose  :
1820 //=======================================================================
1821
1822  void ShapeConstruct_ProjectCurveOnSurface::CheckPoints(Handle(TColgp_HArray1OfPnt)& points,Handle(TColStd_HArray1OfReal)& params,Standard_Real& preci) const
1823 {
1824   Standard_Integer firstElem = points->Lower();
1825   Standard_Integer lastElem  = points->Upper();
1826   Standard_Integer i;
1827   Standard_Integer nbPntDropped = 0;
1828   Standard_Integer lastValid = firstElem; // indice of last undropped point
1829
1830   // will store 0 when the point is to be removed, 1 otherwise
1831   TColStd_Array1OfInteger tmpParam(firstElem, lastElem);
1832   for (i = firstElem; i<=lastElem ; i++)
1833     tmpParam.SetValue(i,1);
1834   Standard_Real DistMin2 = RealLast();
1835   gp_Pnt Prev = points->Value (lastValid);
1836   gp_Pnt Curr;
1837   for (i = firstElem + 1; i <= lastElem ; i ++) {
1838     Curr = points->Value(i);
1839     Standard_Real CurDist2 = Prev.SquareDistance(Curr);
1840     if (CurDist2 < gp::Resolution()) {  // test 0
1841       nbPntDropped ++;
1842       if ( i == lastElem ) tmpParam.SetValue(lastValid, 0); // last point kept
1843       else tmpParam.SetValue(i, 0);    // current dropped, lastValid unchanged
1844     } else {
1845       if (CurDist2 < DistMin2)
1846         DistMin2 = CurDist2;
1847       // lastValid becomes the current (i.e. i)
1848       lastValid = i;
1849       Prev = Curr;
1850     }
1851   }
1852   if (DistMin2 < RealLast())
1853     preci = 0.9 * Sqrt (DistMin2); // preci est la distance min entre les points on la reduit un peu
1854   if (nbPntDropped == 0)
1855     return;
1856
1857 #ifdef OCCT_DEBUG
1858   std::cout << "Warning : removing 3d points for interpolation" << std::endl;
1859 #endif
1860   // Build new HArrays
1861   Standard_Integer newLast = lastElem - nbPntDropped;
1862   if ((newLast - firstElem + 1)  < 2) {
1863 #ifdef OCCT_DEBUG
1864     std::cout << "Too many degenerated points for 3D interpolation" << std::endl;
1865 #endif
1866     return;
1867   }
1868   Handle(TColgp_HArray1OfPnt) newPnts =
1869     new TColgp_HArray1OfPnt(firstElem, newLast);
1870   Handle(TColStd_HArray1OfReal) newParams =
1871     new TColStd_HArray1OfReal(firstElem, newLast);
1872   Standard_Integer newCurr = 1;
1873   for (i = firstElem; i<= lastElem ; i++) {
1874     if (tmpParam.Value(i) == 1) {
1875       newPnts->SetValue(newCurr, points->Value(i));
1876       newParams->SetValue(newCurr, params->Value(i));
1877       newCurr ++;
1878     }
1879   }
1880   points = newPnts;
1881   params = newParams;
1882   // on la reduit un peu
1883 }
1884
1885 //=======================================================================
1886 //function : CheckPoints2d
1887 //purpose  :
1888 //=======================================================================
1889
1890  void ShapeConstruct_ProjectCurveOnSurface::CheckPoints2d(Handle(TColgp_HArray1OfPnt2d)& points,
1891                                                           Handle(TColStd_HArray1OfReal)& params,
1892                                                           Standard_Real& preci) const
1893 {
1894   Standard_Integer firstElem = points->Lower();
1895   Standard_Integer lastElem  = points->Upper();
1896   Standard_Integer i;
1897   Standard_Integer nbPntDropped = 0;
1898   Standard_Integer lastValid = firstElem; // indice of last undropped point
1899
1900   // will store 0 when the point is to be removed, 1 otherwise
1901   TColStd_Array1OfInteger tmpParam(firstElem, lastElem);
1902   for (i = firstElem; i<=lastElem ; i++) {
1903     tmpParam.SetValue(i,1);
1904   }
1905   Standard_Real DistMin2 = RealLast();
1906   gp_Pnt2d Prev = points->Value(lastValid);
1907   gp_Pnt2d Curr;
1908   for (i = firstElem + 1; i<=lastElem ; i++) {
1909     Curr = points->Value(i);
1910     Standard_Real CurDist2 = Prev.SquareDistance(Curr);
1911     if (CurDist2 < gp::Resolution()) {  // test 0
1912       nbPntDropped ++;
1913       if ( i == lastElem ) tmpParam.SetValue(lastValid, 0); // last point kept
1914       else tmpParam.SetValue(i, 0);    // current dropped, lastValid unchanged
1915     } else {
1916       if (CurDist2 < DistMin2)
1917         DistMin2 = CurDist2;
1918       // lastValid becomes the current (i.e. i)
1919       lastValid = i;
1920       Prev = Curr;
1921     }
1922   }
1923   if (DistMin2 < RealLast())
1924     preci = 0.9 * Sqrt (DistMin2);
1925   if (nbPntDropped == 0)
1926     return;
1927
1928 #ifdef OCCT_DEBUG
1929   std::cout << "Warning : removing 2d points for interpolation" << std::endl;
1930 #endif
1931   // Build new HArrays
1932   Standard_Integer newLast = lastElem - nbPntDropped;
1933   if ((newLast - firstElem + 1)  < 2) {
1934 #ifdef OCCT_DEBUG
1935     std::cout << "Too many degenerated points for 2D interpolation" << std::endl;
1936 #endif
1937     //pdn 12.02.99 S4135 Creating pcurve with minimal length.
1938     tmpParam.SetValue(firstElem,1);
1939     tmpParam.SetValue(lastElem,1);
1940     gp_XY  lastPnt = points->Value(lastElem).XY();
1941     lastPnt.Add(gp_XY(preci,preci));
1942     points->SetValue(lastElem,lastPnt);
1943     newLast = firstElem+1;
1944     //return;
1945   }
1946   Handle(TColgp_HArray1OfPnt2d) newPnts =
1947     new TColgp_HArray1OfPnt2d(firstElem, newLast);
1948   Handle(TColStd_HArray1OfReal) newParams =
1949     new TColStd_HArray1OfReal(firstElem, newLast);
1950   Standard_Integer newCurr = 1;
1951   for (i = firstElem; i <= lastElem ; i++) {
1952     if (tmpParam.Value(i) == 1) {
1953 #ifdef OCCT_DEBUG
1954       std::cout << "Point " << i << " : " << points->Value(i).X() << " " << points->Value(i).Y() << " at param " <<  params->Value(i) << std::endl;
1955 #endif
1956       newPnts->SetValue(newCurr, points->Value(i));
1957       newParams->SetValue(newCurr, params->Value(i));
1958       newCurr ++;
1959     }
1960     else {
1961 #ifdef OCCT_DEBUG
1962       std::cout << "Removed " << i << " : " << points->Value(i).X() << " " << points->Value(i).Y() << " at param " <<  params->Value(i) << std::endl;
1963 #endif
1964     }
1965   }
1966   points = newPnts;
1967   params = newParams;
1968 }
1969
1970 //=======================================================================
1971 //function : IsAnIsoparametric
1972 //purpose  :
1973 //=======================================================================
1974 //:S4030: modified for optimization
1975 //:p9 abv 11 Mar 99: PRO7226 #489490: find nearest boundary instead of first one
1976
1977  Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::IsAnIsoparametric(const Standard_Integer nbrPnt,
1978                                                                           const TColgp_SequenceOfPnt& points,
1979                                                                           const TColStd_SequenceOfReal& params,
1980                                                                           Standard_Boolean& isoTypeU,
1981                                                                           Standard_Boolean& p1OnIso,
1982                                                                           gp_Pnt2d& valueP1,
1983                                                                           Standard_Boolean& p2OnIso,
1984                                                                           gp_Pnt2d& valueP2,
1985                                                                           Standard_Boolean& isoPar2d3d,
1986                                                                           Handle(Geom_Curve)& cIso,
1987                                                                           Standard_Real& t1,
1988                                                                           Standard_Real& t2,
1989                                                                           TColStd_Array1OfReal& pout) const
1990 {
1991   try {    // RAJOUT
1992     OCC_CATCH_SIGNALS
1993
1994   Standard_Real prec = Precision::Confusion();//myPreci;
1995
1996   Standard_Boolean isoParam = Standard_False;
1997   isoPar2d3d = Standard_False;
1998
1999   Standard_Real U1, U2, V1, V2;
2000   mySurf->Bounds(U1, U2, V1, V2);
2001
2002   if ( mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
2003     Handle(Geom_RectangularTrimmedSurface) sTrim =
2004       Handle(Geom_RectangularTrimmedSurface)::DownCast(mySurf->Surface());
2005     sTrim->Bounds(U1, U2, V1, V2);
2006   }
2007
2008   gp_Pnt pt;
2009   Standard_Integer mpt[2]; mpt[0] = mpt[1] = 0;
2010   Standard_Real t, tpar[2] = { 0.0, 0.0 }, isoValue=0.;
2011   Standard_Real mindist2;
2012   Standard_Real mind2[2];
2013   mindist2 = mind2[0] = mind2[1] = 4*prec*prec;
2014
2015   p1OnIso = Standard_False;
2016   p2OnIso = Standard_False;
2017   const Bnd_Box* aBox = 0;
2018
2019   for (Standard_Integer j=1; (j<=4) /*&& !isoParam*/; j++) {
2020     Standard_Real isoVal=0.;
2021     Standard_Boolean isoU=Standard_False; //szv#4:S4163:12Mar99 `isoU` must be Standard_Boolean
2022     Handle(Geom_Curve) cI;
2023     Standard_Real tt1, tt2;
2024
2025     if      (j == 1 ) {
2026       if (Precision::IsInfinite(U1)) continue;
2027       cI = mySurf->UIso(U1);
2028       isoU = Standard_True;
2029       isoVal = U1;
2030       aBox = & mySurf->GetBoxUF();
2031     }
2032     else if (j == 2) {
2033       if (Precision::IsInfinite(U2)) continue;
2034       cI = mySurf->UIso(U2);
2035       isoU = Standard_True;
2036       isoVal = U2;
2037       aBox = & mySurf->GetBoxUL();
2038     }
2039     else if (j == 3) {
2040       if (Precision::IsInfinite(V1)) continue;
2041       cI = mySurf->VIso(V1);
2042       isoU = Standard_False;
2043       isoVal = V1;
2044       aBox = & mySurf->GetBoxVF();
2045     }
2046     else if (j == 4) {
2047       if (Precision::IsInfinite(V2)) continue;
2048       cI = mySurf->VIso(V2);
2049       isoU = Standard_False;
2050       isoVal = V2;
2051       aBox = & mySurf->GetBoxVL();
2052     }
2053     if(cI.IsNull())
2054       continue;
2055
2056     if (isoU)  {  tt1 = V1;  tt2 = V2;  }
2057     else       {  tt1 = U1;  tt2 = U2;  }
2058
2059     gp_Pnt ext1, ext2;
2060     cI->D0(tt1, ext1);
2061     cI->D0(tt2, ext2);
2062
2063 // PATCH CKY 9-JUL-1998 : protection contre singularite
2064     gp_Pnt extmi;
2065     cI->D0( (tt1+tt2)/2,extmi);
2066     if (ext1.IsEqual(ext2,prec) && ext1.IsEqual(extmi,prec)) continue;
2067
2068     Standard_Boolean PtEQext1 = Standard_False;
2069     Standard_Boolean PtEQext2 = Standard_False;
2070
2071     Standard_Real currd2[2], tp[2] = {0, 0};
2072     Standard_Integer mp[2];
2073
2074     for (Standard_Integer i=0; i<2; i++) {
2075       mp[i] = 0;
2076       Standard_Integer k = (i == 0 ? 1 : nbrPnt);
2077
2078       // si ext1 == ext2 => valueP1 == valueP2 => vect null plus tard
2079       currd2[i] = points(k).SquareDistance ( ext1 );
2080       if ( currd2[i] <= prec*prec && !PtEQext1) {
2081         mp[i] = 1;
2082         tp[i] = tt1;
2083         PtEQext1 = Standard_True;
2084         continue;
2085       }
2086
2087       currd2[i] = points(k).SquareDistance ( ext2 );
2088       if ( currd2[i] <= prec*prec && !PtEQext2) {
2089         mp[i] = 2;
2090         tp[i] = tt2;
2091         PtEQext2 = Standard_True;
2092         continue;
2093       }
2094
2095       // On evite de projecter sur un iso degenere
2096       // on doit egalement le faire pour l apex du cone
2097       if (mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) && !isoU) {
2098         continue;
2099       }
2100
2101       if(aBox->IsOut(points(k))) continue;
2102
2103       Standard_Real Cf = cI->FirstParameter();
2104       Standard_Real Cl = cI->LastParameter();
2105       if (Precision::IsInfinite(Cf))  Cf = -1000;
2106       if (Precision::IsInfinite(Cl))  Cl = +1000;
2107
2108       ShapeAnalysis_Curve sac;
2109       Standard_Real dist = sac.Project (cI,points(k),prec,pt,t,Cf,Cl);
2110       currd2[i] = dist * dist;
2111       if ((dist <= prec) && (t>= Cf) && (t<=Cl)) {
2112         mp[i]  = 3;
2113         tp[i] = t;
2114       }
2115     }
2116
2117     //:e7 abv 21 Apr 98: ProSTEP TR8, r0501_pe #56679:
2118     // avoid possible null-length curves
2119     if ( mp[0] >0 && mp[1] >0 &&
2120          Abs ( tp[0] - tp[1] ) < Precision::PConfusion() ) continue;
2121
2122
2123     if (mp[0] > 0 &&
2124         ( ! p1OnIso || currd2[0] < mind2[0] ) ) {
2125       p1OnIso = Standard_True;
2126       mind2[0] = currd2[0]; // LP2.stp #105899: FLT_INVALID_OPERATION on Windows 7 VC 9 Release mode on the whole file
2127       if (isoU) valueP1.SetCoord(isoVal, tp[0]);
2128       else      valueP1.SetCoord(tp[0], isoVal);
2129     }
2130
2131     if (mp[1] > 0 &&
2132         ( ! p2OnIso || currd2[1] < mind2[1] ) ) {
2133       p2OnIso = Standard_True;
2134       mind2[1] = currd2[1];
2135       if (isoU) valueP2.SetCoord(isoVal, tp[1]);
2136       else      valueP2.SetCoord(tp[1], isoVal);
2137     }
2138
2139     if ( mp[0] <=0 || mp[1] <=0 ) continue;
2140
2141     Standard_Real md2 = currd2[0] + currd2[1];
2142     if ( mindist2 <= md2 ) continue;
2143
2144     mindist2 = md2;
2145     mpt[0] = mp[0];
2146     mpt[1] = mp[1];
2147     tpar[0] = tp[0];
2148     tpar[1] = tp[1];
2149     isoTypeU = isoU;
2150     isoValue = isoVal;
2151     cIso = cI;
2152     t1 = tt1;
2153     t2 = tt2;
2154   }
2155
2156   // probablely it concerns an isoparametrics
2157   if ( mpt[0] >0 && mpt[1] >0 ) {
2158
2159     p1OnIso = p2OnIso = Standard_True;
2160     if (isoTypeU) {
2161       valueP1.SetCoord(isoValue, tpar[0]);
2162       valueP2.SetCoord(isoValue, tpar[1]);
2163     }
2164     else {
2165       valueP1.SetCoord(tpar[0], isoValue);
2166       valueP2.SetCoord(tpar[1], isoValue);
2167     }
2168
2169     if ( mpt[0] != 3 && mpt[1] != 3 ) {
2170       isoPar2d3d = Standard_True;
2171       for (Standard_Integer i=2; i < nbrPnt && isoPar2d3d; i++){
2172         if (tpar[1] > tpar[0])  t = params(i);
2173         else                    t = t1+t2-params(i);
2174         cIso->D0(t, pt);
2175         if (!points(i).IsEqual(pt, prec)) isoPar2d3d = Standard_False;
2176       }
2177     }
2178
2179     if (isoPar2d3d) isoParam = Standard_True;
2180     else {
2181       Standard_Real prevParam = tpar[0];
2182       Standard_Real Cf, Cl;
2183       Standard_Boolean isoByDistance = Standard_True;
2184       Cf = cIso->FirstParameter();
2185       Cl = cIso->LastParameter();
2186       if (Precision::IsInfinite(Cf))  Cf = -1000;
2187       if (Precision::IsInfinite(Cl))  Cl = +1000;
2188
2189       ShapeAnalysis_Curve sac;
2190       for (Standard_Integer i=2; i < nbrPnt && isoByDistance; i++) {
2191         Standard_Real dist = sac.NextProject (prevParam,cIso,points(i),
2192                                               prec,pt,t,Cf,Cl,
2193                                               Standard_False); //:j8 abv 10.12.98: TR10 r0501_db.stp #9423: avoid adjusting to ends
2194         prevParam = t;
2195         pout(i)=t;
2196         if( (dist > prec) || (t < Cf) || (t > Cl) )
2197           isoByDistance = Standard_False;
2198       }
2199       if (isoByDistance) isoParam = Standard_True;
2200     }
2201   }
2202 /*  if (!isoParam) {    CKY 29-mai-1997 : garder tout ce qu on peut ?
2203     p1OnIso = Standard_False;
2204     p2OnIso = Standard_False;
2205   }  */
2206   return isoParam;
2207   }  // RAJOUT
2208   catch(Standard_Failure const& anException) {
2209 #ifdef OCCT_DEBUG
2210 //  pb : on affiche ce qu on peut
2211     for (Standard_Integer numpnt = 1; numpnt <= nbrPnt; numpnt ++) {
2212       std::cout<<"["<<numpnt<<"]param="<<params(numpnt)<<" point=("<<
2213         points(numpnt).X()<<"  "<<points(numpnt).Y()<<"  "<<points(numpnt).Z()<<")"<<std::endl;
2214     }
2215     std::cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::IsAnIsoparametric(): Exception: ";
2216     anException.Print(std::cout); std::cout << std::endl;
2217 #endif
2218     (void)anException;
2219     return Standard_False;
2220   }
2221 }
2222
2223 /* S4135 : BestExtremum is commented after IsAnIsoparametric works with Precision::Confusion()
2224 //=======================================================================
2225 //function : BestExtremum
2226 //purpose  : auxiliaire prenant le meilleur extremum si ISO car doute possible
2227 //=======================================================================
2228
2229  gp_Pnt2d ShapeConstruct_ProjectCurveOnSurface::BestExtremum(const gp_Pnt2d& P2iso,const gp_Pnt& P3ext,const gp_Pnt& P3next) const
2230 {
2231 //  P2iso a ete calcule depuis P3ext sur une iso externe de la surface
2232 //  En principe bon mais circularite possible ... et IsU/VClosed faillible
2233 //    (si baillement 1e-4 ou 1e-5, on est dedans !). DONC
2234 //  1/ on privilegie l iso mais a tout hasard on verifie si Surf meilleur
2235 //  2/ si iso, attention a la circularite (cas limite)
2236
2237 //  NB : si isoParam, on suppose que P2iso est bon (car il y en a 2). A voir...
2238
2239 //  D abord, calcul p2ext depuis la surface. choix surface/iso
2240   return P2iso;
2241   Standard_Real prec = Precision::Confusion();//myPreci;
2242   gp_Pnt2d P2cal = mySurf->ValueOfUV(P3ext, prec);
2243   gp_Pnt   P3cal = mySurf->Value (P2cal);
2244   Standard_Real dcal = P3ext.Distance (P3cal);
2245   Standard_Real dnxt = P3ext.Distance (P3next);
2246   if (dcal > dnxt) return P2iso;    // en fait protection sur BUG (PRO8468)
2247
2248 //  On choisit entre P2iso et P2cal, le plus proche de P2next ... !!!
2249   gp_Pnt2d P2next = mySurf->ValueOfUV(P3next, prec);
2250   if (P2next.Distance(P2cal) < P2next.Distance(P2iso)) return P2cal;
2251   return P2iso;
2252 }
2253 */