src/ShapeConstruct/ShapeConstruct_ProjectCurveOnSurface.cxx

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