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