0028643: Coding rules - eliminate GCC compiler warnings -Wmisleading-indentation
[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 const& anException) {
475 #ifdef OCCT_DEBUG
476     cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::PerformByProjLib(): Exception: ";
477     anException.Print(cout); cout << endl;
478 #endif
479     (void)anException;
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      {
708        if ( myCashe3d[j].SquareDistance (aP[i] ) < aTol2)
709        {
710          aP2d[i] = mySurf->NextValueOfUV (myCashe2d[j], aP[i], theTol, theTol);
711          aSavedPointNum = i;
712          aSavedPoint = myCashe2d[j];
713          if (i == 0)
714            isFromCashe = Standard_True;
715          break;
716        }
717      }
718
719      if (j >= myNbCashe)
720      {
721        aP2d[i] = mySurf->ValueOfUV(aP[i], theTol);
722      }
723
724      Standard_Real aDist = mySurf->Gap();
725      Standard_Real aCurDist = aDist * aDist;
726      if (aTol2 < aDist * aDist)
727      {
728        aTol2 = aCurDist;
729      }
730    }
731
732    if ( isPeriodicU || isPeriodicV )
733    {
734      // Compute second and last but one c2d points.
735      for(i = 1; i < 3; i++)
736      {
737        Standard_Integer j;
738        for (j = 0; j < myNbCashe; ++j)
739        {
740          if ( myCashe3d[j].SquareDistance (aP[i] ) < aTol2)
741          {
742            aP2d[i] = mySurf->NextValueOfUV (myCashe2d[j], aP[i], theTol, theTol);
743            aSavedPointNum = i;
744            aSavedPoint = myCashe2d[j];
745            break;
746          }
747        }
748
749        if (j >= myNbCashe)
750        {
751          aP2d[i] = mySurf->ValueOfUV(aP[i], theTol);
752        }
753
754        Standard_Real aDist = mySurf->Gap();
755        Standard_Real aCurDist = aDist * aDist;
756        if (aTol2 < aDist * aDist)
757        {
758          aTol2 = aCurDist;
759        }
760      }
761
762      if (isPeriodicU)
763      {
764        isRecompute = fixPeriodictyTroubles(&aP2d[0], 1 /* X Coord */, mySurf->Surface()->UPeriod(), aSavedPointNum, aSavedPoint.X());
765      }
766
767      if (isPeriodicV)
768      {
769        isRecompute = fixPeriodictyTroubles(&aP2d[0], 2 /* Y Coord */, mySurf->Surface()->VPeriod(), aSavedPointNum, aSavedPoint.Y());
770      }
771    }
772
773    if (isRecompute && mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_SphericalSurface))) {
774      // Do not try to make line, because in this case may be very special case when 3d curve
775      // go over the pole of, e.g., sphere, and partly lies along seam
776      // (see ApproxPCurve() for more information).
777      return 0;
778    }
779
780    thePnt2ds.SetValue(1, aP2d[0]);
781    thePnt2ds.SetValue(nb, aP2d[3]);
782
783    // Restore old tolerance in 2d space to avoid big gap cases.
784    aTol2 = anOldTol2;
785    // Check that straight line in 2d with parameterisation as in 3d will fit 
786    // fit 3d curve at all points.
787    Standard_Real dPar = theparams(nb) - theparams(1);
788    if ( Abs(dPar) < Precision::PConfusion() )
789      return 0;
790    gp_Vec2d aVec0 (aP2d[0], aP2d[3]);
791    gp_Vec2d aVec = aVec0 / dPar;
792    Handle(Geom_Surface) aSurf = mySurf->Surface();
793    Standard_Boolean isNormalCheck = aSurf->IsCNu(1) && aSurf->IsCNv(1);
794    if (isNormalCheck) {
795      for(i = 1; i <= nb; i++)
796      {
797        gp_XY aCurPoint = aP2d[0].XY() + aVec.XY() * (theparams(i) - theparams(1));
798        gp_Pnt aCurP;
799        gp_Vec aNormalVec, aDu, aDv;
800        aSurf->D1(aCurPoint.X(), aCurPoint.Y(), aCurP, aDu, aDv);
801        aNormalVec = aDu ^ aDv;
802        if (aNormalVec.SquareMagnitude() < Precision::SquareConfusion()) {
803          isNormalCheck = Standard_False;
804          break;
805        }
806        gp_Lin aNormalLine(aCurP, gp_Dir(aNormalVec));
807        Standard_Real aDist = aNormalLine.Distance(thepoints(i));
808        if (aDist > theTol)
809          return 0;
810      }
811    }
812    if (!isNormalCheck) {
813      Standard_Real aFirstPointDist = mySurf->Surface()->Value(aP2d[0].X(), aP2d[0].Y()). 
814                                      SquareDistance(thepoints(1));
815      aTol2 = Max(aTol2, aTol2 * 2 * aFirstPointDist);
816      for(i = 2; i <  nb; i++)
817      {
818        gp_XY aCurPoint = aP2d[0].XY() + aVec.XY() * (theparams(i) - theparams(1));
819        gp_Pnt aCurP;
820        aSurf->D0(aCurPoint.X(), aCurPoint.Y(), aCurP);
821        Standard_Real aDist1 = aCurP.SquareDistance(thepoints(i));
822      
823        if(Abs (aFirstPointDist - aDist1) > aTol2)
824          return 0;
825      }
826    }
827
828    // check if pcurve can be represented by Geom2d_Line (parameterised by length)
829    Standard_Real aLLength = aVec0.Magnitude();
830    if ( Abs (aLLength - dPar) <= Precision::PConfusion() )
831    {
832      gp_XY aDirL = aVec0.XY() / aLLength;
833      gp_Pnt2d aPL (aP2d[0].XY() - theparams(1) * aDirL);
834      return new Geom2d_Line (aPL, gp_Dir2d(aDirL));
835    }
836
837    // create straight bspline
838    TColgp_Array1OfPnt2d aPoles(1, 2);
839    aPoles(1) = aP2d[0];
840    aPoles(2) = aP2d[3];
841
842    TColStd_Array1OfReal aKnots(1,2);
843    aKnots(1) = theparams(1);
844    aKnots(2) = theparams(theparams.Length());
845
846    TColStd_Array1OfInteger aMults(1,2);
847    aMults(1) = 2;
848    aMults(2) = 2;
849    Standard_Integer aDegree = 1;
850    Handle(Geom2d_BSplineCurve) abspl2d =
851      new Geom2d_BSplineCurve (aPoles, aKnots, aMults, aDegree);
852    return abspl2d;
853  }
854
855 //=======================================================================
856 //function : ApproxPCurve
857 //purpose  : 
858 //=======================================================================
859
860   Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::ApproxPCurve(const Standard_Integer nbrPnt,
861                                                                      const TColgp_Array1OfPnt& points,
862                                                                      const TColStd_Array1OfReal& params,
863                                                                      TColgp_Array1OfPnt2d& pnt2d,
864                                                                      Handle(Geom2d_Curve)& c2d) 
865 {
866   // for performance, first try to handle typical case when pcurve is straight
867   Standard_Boolean isRecompute = Standard_False;
868   Standard_Boolean isFromCasheLine = Standard_False;
869   c2d = getLine(points, params, pnt2d, myPreci, isRecompute, isFromCasheLine);
870   if(!c2d.IsNull())
871   {
872     // fill cashe
873     Standard_Boolean ChangeCycle = Standard_False;
874     if(myNbCashe>0 && myCashe3d[0].Distance(points(1)) > myCashe3d[0].Distance(points(nbrPnt)) &&
875        myCashe3d[0].Distance(points(nbrPnt))<Precision::Confusion())
876       ChangeCycle = Standard_True;
877     myNbCashe = 2;
878     if(ChangeCycle) {
879       myCashe3d[0] = points(1);
880       myCashe3d[1] = points(nbrPnt);
881       myCashe2d[0] = pnt2d(1);
882       myCashe2d[1] = pnt2d(nbrPnt);
883     }
884     else {
885       myCashe3d[1] = points(1);
886       myCashe3d[0] = points(nbrPnt);
887       myCashe2d[1] = pnt2d(1);
888       myCashe2d[0] = pnt2d(nbrPnt);
889     }
890     return Standard_True;
891   }
892     Standard_Boolean isDone = Standard_True;
893   // test if the curve 3d is a boundary of the surface 
894   // (only for Bezier or BSpline surface)
895   
896   Standard_Boolean isoParam, isoPar2d3d, isoTypeU, p1OnIso, p2OnIso, isoclosed;
897   gp_Pnt2d valueP1, valueP2;
898   Handle(Geom_Curve) cIso;
899   Standard_Real t1, t2;
900   
901   Handle(Standard_Type) sType = mySurf->Surface()->DynamicType();
902   Standard_Boolean isAnalytic = Standard_True;
903   if (sType == STANDARD_TYPE(Geom_BezierSurface) || sType == STANDARD_TYPE(Geom_BSplineSurface)) isAnalytic = Standard_False; 
904   Standard_Real uf, ul, vf, vl;
905   mySurf->Surface()->Bounds(uf, ul, vf, vl);
906   isoclosed = Standard_False;
907   TColStd_Array1OfReal pout(1, nbrPnt);
908   
909   isoParam = IsAnIsoparametric(nbrPnt, points, params, 
910                                isoTypeU, p1OnIso, valueP1, p2OnIso, valueP2,
911                                isoPar2d3d, cIso, t1, t2, pout);
912   
913   // projection of the points on surfaces
914   
915   gp_Pnt p3d;
916   gp_Pnt2d p2d;
917   Standard_Integer i;
918   Standard_Real isoValue=0., isoPar1=0., isoPar2=0., tPar=0., tdeb,tfin;
919   Standard_Real Cf, Cl, parf, parl; //szv#4:S4163:12Mar99 dist not needed
920   
921   //  Le calcul part-il dans le bon sens, c-a-d deb et fin dans le bon ordre ?
922   //  Si uclosed et iso en V, attention isoPar1 ET/OU 2 peut toucher la fermeture
923   if(isoParam){
924     if(isoTypeU){
925       isoValue = valueP1.X();
926       isoPar1 = valueP1.Y();
927       isoPar2 = valueP2.Y();
928       isoclosed = mySurf->IsVClosed(myPreci);//#78 rln 12.03.99 S4135
929       parf = vf;  parl = vl;
930     }
931     else { 
932       isoValue = valueP1.Y();
933       isoPar1 = valueP1.X();
934       isoPar2 = valueP2.X();
935       isoclosed = mySurf->IsUClosed(myPreci);//#78 rln 12.03.99 S4135
936       parf = uf;  parl = ul;
937     }
938     if (!isoPar2d3d && !isAnalytic) {
939       Cf = cIso->FirstParameter();
940       Cl = cIso->LastParameter();
941       if (Precision::IsInfinite(Cf))    Cf = -1000;
942       if (Precision::IsInfinite(Cl))    Cl = +1000;
943       //pdn S4030 optimizing and fix isopar case on PRO41323
944       tdeb = pout(2);
945       //    dist = ShapeAnalysis_Curve().Project (cIso,points(2),myPreci,pt,tdeb,Cf,Cl);
946       //  Chacun des par1 ou par2 est-il sur un bord. Attention first/last : recaler
947       if (isoclosed && (isoPar1 == parf || isoPar1 == parl)) {
948         if (Abs(tdeb-parf) < Abs(tdeb-parl)) isoPar1 = parf;
949         else isoPar1 = parl;
950         if (isoTypeU) valueP1.SetY (isoPar1);
951         else          valueP1.SetX (isoPar1);
952       }
953       if (isoclosed && (isoPar2 == parf || isoPar2 == parl)) {
954         //pdn S4030 optimizing and fix isopar case on PRO41323
955         tfin = pout(nbrPnt-1);
956         //dist =  ShapeAnalysis_Curve().Project (cIso,points(nbrPnt-1),myPreci,pt,tfin,Cf,Cl);
957         if (Abs(tfin-parf) < Abs(tfin-parl)) isoPar2 = parf;
958         else isoPar2 = parl;
959         if (isoTypeU) valueP2.SetY (isoPar2);
960         else          valueP2.SetX (isoPar2);
961       }
962       
963       //  Interversion Par1/Par2 (ne veut que si les 2 sont sur les bords ...)
964       //  Est-ce encore necessaire apres ce qui vient d etre fait ?
965       
966       // PTV 05.02.02 fix for translation face from 12_hp_mouse (PARASOLID) face 24008
967       // if curve is periodic do not change the points
968       // skl change "if" for pout(nbrPnt-1) 19.11.2003
969       if (!isoclosed) {
970         if( (Abs(tdeb-isoPar1)>Abs(tdeb-isoPar2)) &&
971             (Abs(pout(nbrPnt-1)-isoPar2)>Abs(pout(nbrPnt-1)-isoPar1)) ) {
972           gp_Pnt2d valueTmp = valueP1;
973           valueP1 = valueP2;  valueP2 = valueTmp;
974           if (isoTypeU) {
975             isoValue = valueP1.X();
976             isoPar1 = valueP1.Y();
977             isoPar2 = valueP2.Y();
978           }
979           else { 
980             isoValue = valueP1.Y();
981             isoPar1 = valueP1.X();
982             isoPar2 = valueP2.X();
983           }
984           //  Fin calcul sens de courbe iso
985         }
986       } // end of fix check 05.02.02
987     }
988   }
989   
990   //  Si pas isoParam, on a quand meme du p1OnIso/p2OnIso possible ... !!!
991   //  (utile pour detromper bug de projection). Mais detromper aussi circularite
992   //else {
993     //if (p1OnIso) valueP1 =
994     //BestExtremum (valueP1,points(1),points(2));
995     //if (p2OnIso) valueP2 =
996     //BestExtremum (valueP2,points(nbrPnt),points(nbrPnt-1));
997     //}
998   
999   Standard_Real gap = myPreci; //:q1
1000   Standard_Boolean ChangeCycle = Standard_False; //skl for OCC3430
1001   // auxiliaruy variables to shift 2dcurve, according to previous
1002   Standard_Boolean isFromCashe = Standard_False;
1003   gp_Pnt2d aSavedPoint;
1004   if( myNbCashe>0 && myCashe3d[0].Distance(points(1))>myCashe3d[0].Distance(points(nbrPnt)) )
1005     //if(myCashe3d[0].Distance(points(nbrPnt))<myPreci)
1006     if(myCashe3d[0].Distance(points(nbrPnt))<Precision::Confusion())
1007       ChangeCycle = Standard_True;
1008   //for( i = 1; i <= nbrPnt; i ++) {
1009   for(Standard_Integer ii=1; ii<=nbrPnt; ii++) {
1010     if(ChangeCycle) //skl for OCC3430
1011       i=nbrPnt-ii+1;
1012     else
1013       i=ii;
1014     p3d = points(i);
1015     if (isoParam) {
1016       
1017       if (isoPar2d3d) {
1018         if (isoPar2 > isoPar1) tPar = params(i);
1019         else                   tPar = t1 + t2 - params(i);
1020       } else if (!isAnalytic) {
1021         // projection to iso
1022         if      (i==1)      tPar = isoPar1;
1023         else if (i==nbrPnt) tPar = isoPar2;
1024         else {
1025           tPar = pout(i);
1026           //:S4030  ShapeAnalysis_Curve().Project (cIso,p3d,myPreci,pt,tPar,Cf,Cl); //szv#4:S4163:12Mar99 `dist=` not needed
1027         }
1028       }
1029       
1030       if (!isoPar2d3d && isAnalytic) {
1031         if      (i == 1)      p2d = valueP1;
1032         else if (i == nbrPnt) p2d = valueP2;
1033         else {
1034           p2d = mySurf->NextValueOfUV(p2d,p3d, myPreci, //%12 pdn 15.02.99 optimizing
1035                                       Precision::Confusion()+1000*gap); //:q1
1036           gap = mySurf->Gap();
1037         }
1038       } else {
1039         if(isoTypeU)  {  p2d.SetX(isoValue);  p2d.SetY(tPar);     }
1040         else          {  p2d.SetX(tPar);      p2d.SetY(isoValue); }
1041       }
1042     }
1043     
1044     else {
1045       if     ( (i == 1)      && p1OnIso)  p2d = valueP1;
1046       else if( (i == nbrPnt) && p2OnIso)  p2d = valueP2;
1047       else  {// general case (not an iso)  mais attention aux singularites !
1048         // first and last points are already computed by getLine()
1049         if ( (i == 1 || i == nbrPnt))
1050         {
1051           if (!isRecompute)
1052           {
1053             p2d = pnt2d(i);
1054             gap = mySurf->Gap();
1055             if (i == 1) {
1056               isFromCashe = isFromCasheLine;
1057               aSavedPoint = p2d;
1058             }
1059             continue;
1060           }
1061           else
1062           {
1063             //:q9 abv 23 Mar 99: use cashe as 1st approach
1064             Standard_Integer j; // svv #1
1065             for (j = 0; j < myNbCashe; ++j)
1066             {
1067               if ( myCashe3d[j].SquareDistance ( p3d ) < myPreci*myPreci )
1068               {
1069                 p2d = mySurf->NextValueOfUV (myCashe2d[j], p3d, myPreci, Precision::Confusion()+gap);
1070                 if (i == 1)
1071                 {
1072                   isFromCashe = Standard_True;
1073                   aSavedPoint = myCashe2d[j];
1074                 }
1075                 break;
1076               }
1077             }
1078             if (j >= myNbCashe)
1079             {
1080               p2d = mySurf->ValueOfUV(p3d, myPreci);
1081             }
1082           }
1083         }
1084         else {
1085           p2d = mySurf->NextValueOfUV (p2d, p3d, myPreci, //:S4030: optimizing
1086                                        Precision::Confusion()+1000*gap); //:q1
1087         }
1088         gap = mySurf->Gap();
1089       }
1090     }
1091     pnt2d (i) = p2d;
1092     if ( ii > 1 ) {
1093       if(ChangeCycle)
1094         p2d.SetXY ( 2. * p2d.XY() - pnt2d(i+1).XY() );
1095       else
1096         p2d.SetXY ( 2. * p2d.XY() - pnt2d(i-1).XY() );
1097     }
1098   }
1099
1100   //pdn %12 11.02.99 PRO9234 entity 15402
1101   if (!isoPar2d3d) {
1102     mySurf->ProjectDegenerated(nbrPnt,points,pnt2d,myPreci,Standard_True);
1103     mySurf->ProjectDegenerated(nbrPnt,points,pnt2d,myPreci,Standard_False);
1104   }
1105   
1106   //  attention aux singularites ... (hors cas iso qui les traite deja)
1107   //  if (!isoParam) {
1108     //    p2d = pnt2d (1);
1109     //    if (mySurf->ProjectDegenerated (points(1),myPreci,pnt2d (2),p2d))
1110     //      pnt2d (1) = p2d;
1111     //    p2d = pnt2d (nbrPnt);
1112     //    if (mySurf->ProjectDegenerated (points(nbrPnt),myPreci,pnt2d (nbrPnt-1),p2d))
1113     //      pnt2d (nbrPnt) = p2d;
1114     //  }
1115   
1116   // Si la surface est UCLosed et VClosed, on recadre les points
1117   // algo un peu complique, on retarde l implementation
1118   Standard_Real Up = ul - uf;
1119   Standard_Real Vp = vl - vf;
1120   Standard_Real dist2d;
1121 #ifdef OCCT_DEBUG
1122   if (mySurf->IsUClosed(myPreci) && mySurf->IsVClosed(myPreci)) {//#78 rln 12.03.99 S4135
1123     cout << "WARNING : Recadrage incertain sur U & VClosed" << endl;
1124   }
1125 #endif
1126   // Si la surface est UCLosed, on recadre les points
1127   if (mySurf->IsUClosed(myPreci)) {//#78 rln 12.03.99 S4135
1128     // Premier point dans le domain [uf, ul]
1129     Standard_Real prevX, firstX = pnt2d (1).X();
1130     if (!isFromCashe) {
1131       // do not shift 2dcurve, if it connects to previous
1132       while (firstX < uf)  {  firstX += Up;   pnt2d (1).SetX(firstX);  }
1133       while (firstX > ul)  {  firstX -= Up;   pnt2d (1).SetX(firstX);  }
1134     }
1135     // shift first point, according to cashe
1136     if (mySurf->Surface()->IsUPeriodic() && isFromCashe) {
1137       Standard_Real aMinParam = uf, aMaxParam = ul;
1138       while (aMinParam > aSavedPoint.X()) {
1139         aMinParam -= Up;
1140         aMaxParam -= Up;
1141       }
1142       while (aMaxParam < aSavedPoint.X()) {
1143         aMinParam += Up;
1144         aMaxParam += Up;
1145       }
1146       Standard_Real aShift = ShapeAnalysis::AdjustToPeriod(firstX, aMinParam, aMaxParam);
1147       firstX += aShift;
1148       pnt2d(1).SetX(firstX);
1149     }
1150     prevX = firstX;
1151     
1152     //:97 abv 1 Feb 98: treat case when curve is whole out of surface bounds
1153     Standard_Real minX = firstX, maxX = firstX;
1154     
1155     // On decalle toujours le suivant
1156     for (i = 2; i <= nbrPnt; i++) {
1157       //      dist2d = pnt2d (i-1).Distance(pnt2d (i));
1158       Standard_Real CurX = pnt2d (i).X();
1159       dist2d = Abs (CurX - prevX);
1160       if (dist2d > ( Up / 2) ) {
1161         if        (CurX > prevX + Up/2) {
1162           while (CurX > prevX + Up/2) {  CurX -= Up;  pnt2d (i).SetX (CurX);  }
1163         } else if (CurX < prevX - Up/2) {
1164           while (CurX < prevX - Up/2) {  CurX += Up;  pnt2d (i).SetX (CurX);  }
1165         }
1166         
1167       }
1168       prevX = CurX;
1169       if ( minX > CurX ) minX = CurX;      //:97
1170       else if ( maxX < CurX ) maxX = CurX; //:97
1171     }
1172     
1173     //:97
1174     if (!isFromCashe) {
1175       // do not shift 2dcurve, if it connects to previous
1176       Standard_Real midX = 0.5 * ( minX + maxX );
1177       Standard_Real shiftX=0.;
1178       if ( midX > ul ) shiftX = -Up;
1179       else if ( midX < uf ) shiftX = Up;
1180       if ( shiftX != 0. ) 
1181         for ( i=1; i <= nbrPnt; i++ ) pnt2d(i).SetX ( pnt2d(i).X() + shiftX );
1182       }
1183   }
1184   // Si la surface est VCLosed, on recadre les points
1185   // Same code as UClosed : optimisation souhaitable !!
1186   // CKY : d abord un code IDENTIQUE A UClosed; PUIS le special Seam ...
1187   // Si la surface est UCLosed, on recadre les points
1188   //
1189   //#69 rln 01.03.99 S4135 bm2_sd_t4-A.stp entity 30
1190   //#78 rln 12.03.99 S4135
1191   if (mySurf->IsVClosed(myPreci) || mySurf->Surface()->IsKind (STANDARD_TYPE (Geom_SphericalSurface))) {
1192     // Premier point dans le domain [vf, vl]
1193     Standard_Real prevY, firstY = pnt2d (1).Y();
1194     if (!isFromCashe) {
1195       // do not shift 2dcurve, if it connects to previous
1196       while (firstY < vf)  {  firstY += Vp;  pnt2d (1).SetY(firstY);  }
1197       while (firstY > vl)  {  firstY -= Vp;  pnt2d (1).SetY(firstY);  }
1198     }
1199     // shift first point, according to cashe
1200     if (mySurf->Surface()->IsVPeriodic() && isFromCashe) {
1201       Standard_Real aMinParam = vf, aMaxParam = vl;
1202       while (aMinParam > aSavedPoint.Y()) {
1203         aMinParam -= Vp;
1204         aMaxParam -= Vp;
1205       }
1206       while (aMaxParam < aSavedPoint.Y()) {
1207         aMinParam += Vp;
1208         aMaxParam += Vp;
1209       }
1210       Standard_Real aShift = ShapeAnalysis::AdjustToPeriod(firstY, aMinParam, aMaxParam);
1211       firstY += aShift;
1212       pnt2d(1).SetY(firstY);
1213     }
1214     prevY = firstY;
1215     
1216     //:97 abv 1 Feb 98: treat case when curve is whole out of surface bounds
1217     Standard_Real minY = firstY, maxY = firstY;
1218     
1219     // On decalle toujours le suivant
1220     for (i = 2; i <= nbrPnt; i ++) {
1221       //      dist2d = pnt2d (i-1).Distance(pnt2d (i));
1222       Standard_Real CurY = pnt2d (i).Y();
1223       dist2d = Abs (CurY - prevY);
1224       if (dist2d > ( Vp / 2) ) {
1225         if        (CurY > prevY + Vp/2) {
1226           while (CurY > prevY + Vp/2) {  CurY -= Vp;  pnt2d (i).SetY (CurY);  }
1227         } else if (CurY < prevY - Vp/2) {
1228           while (CurY < prevY - Vp/2) {  CurY += Vp;  pnt2d (i).SetY (CurY);  }
1229         }
1230       }
1231       prevY = CurY;
1232       if ( minY > CurY ) minY = CurY;      //:97
1233       else if ( maxY < CurY ) maxY = CurY; //:97
1234     }
1235     
1236     //:97
1237     if (!isFromCashe) {
1238       // do not shift 2dcurve, if it connects to previous
1239       Standard_Real midY = 0.5 * ( minY + maxY );
1240       Standard_Real shiftY=0.;
1241       if ( midY > vl ) shiftY = -Vp;
1242       else if ( midY < vf ) shiftY = Vp;
1243       if ( shiftY != 0. ) 
1244         for ( i=1; i <= nbrPnt; i++ ) pnt2d(i).SetY ( pnt2d(i).Y() + shiftY );
1245       }
1246   }
1247   
1248   //#69 rln 01.03.99 S4135 bm2_sd_t4-A.stp entity 30
1249   //#78 rln 12.03.99 S4135
1250   if (mySurf->IsVClosed(myPreci) || mySurf->Surface()->IsKind (STANDARD_TYPE (Geom_SphericalSurface))) {
1251     for (i = 2; i <= nbrPnt; i++) {
1252       //#1 rln 11/02/98 ca_exhaust.stp entity #9869 dist2d = pnt2d (i-1).Distance(pnt2d (i));
1253       dist2d = Abs (pnt2d(i).Y() - pnt2d(i - 1).Y());
1254       if (dist2d > ( Vp / 2) ) {
1255         // ATTENTION : il faut regarder ou le decalage se fait.
1256         // si plusieurs points sont decalles, il faut plusieurs passes
1257         // pour obtenir un resultat correct.
1258         // NOT YET IMPLEMENTED
1259         
1260         // one of those point is incorrectly placed
1261         // i.e on the wrong side of the "seam"
1262         // on prend le point le plus pres des bords vf ou vl
1263         Standard_Boolean prevOnFirst = Standard_False;
1264         Standard_Boolean prevOnLast  = Standard_False;
1265         Standard_Boolean currOnFirst = Standard_False;
1266         Standard_Boolean currOnLast  = Standard_False;
1267         
1268         //  .X ?  plutot .Y ,  non ?
1269         Standard_Real distPrevVF = Abs(pnt2d (i-1).Y() - vf);
1270         Standard_Real distPrevVL = Abs(pnt2d (i-1).Y() - vl);
1271         Standard_Real distCurrVF = Abs(pnt2d (i).Y() - vf);
1272         Standard_Real distCurrVL = Abs(pnt2d (i).Y() - vl);
1273         
1274         Standard_Real theMin = distPrevVF;
1275         prevOnFirst = Standard_True;
1276         if (distPrevVL < theMin) {
1277           theMin = distPrevVL;
1278           prevOnFirst = Standard_False;
1279           prevOnLast  = Standard_True;
1280         }
1281         if (distCurrVF < theMin) {
1282           theMin = distCurrVF;
1283           prevOnFirst = Standard_False;
1284           prevOnLast  = Standard_False;
1285           currOnFirst = Standard_True;
1286         }
1287         if (distCurrVL < theMin) {
1288           theMin = distCurrVL;
1289           prevOnFirst = Standard_False;
1290           prevOnLast  = Standard_False;
1291           currOnFirst = Standard_False;
1292           currOnLast  = Standard_True;
1293         }
1294         //  Modifs RLN/Nijni  3-DEC-1997
1295         if (prevOnFirst) {
1296           // on decalle le point (i-1) en V Last
1297           gp_Pnt2d newPrev(pnt2d (i-1).X(), vf); // instead of  vl RLN/Nijni
1298           pnt2d (i-1) = newPrev;
1299         }
1300         else if (prevOnLast) {
1301           // on decalle le point (i-1) en V first
1302           gp_Pnt2d newPrev(pnt2d (i-1).X(), vl); // instead of  vf RLN/Nijni
1303           pnt2d (i-1) = newPrev;
1304         }
1305         else if (currOnFirst) {
1306           // on decalle le point (i) en V Last
1307           gp_Pnt2d newCurr(pnt2d (i).X(),vf);  // instead of vl  RLN/Nijni
1308           pnt2d (i) = newCurr;
1309         }
1310         else if (currOnLast) {
1311           // on decalle le point (i) en V First
1312           gp_Pnt2d newCurr(pnt2d (i).X(), vl); // instead of vf  RLN/Nijni
1313           pnt2d (i) = newCurr;
1314         }
1315         // on verifie
1316 #ifdef OCCT_DEBUG
1317         dist2d = pnt2d (i-1).Distance(pnt2d (i));
1318         if (dist2d > ( Vp / 2) ) {
1319           cout << "Echec dans le recadrage" << endl;
1320         }
1321 #endif
1322       }
1323     }
1324   }
1325   
1326   //:c0 abv 20 Feb 98: treat very special case when 3d curve
1327   // go over the pole of, e.g., sphere, and partly lies along seam.
1328   // 2d representation of such a curve should consist of 3 parts - one on
1329   // regular part of surface (interior), one part along degenerated boundary
1330   // and one along seam.
1331   // Since it cannot be adjusted later by arranging pcurves (curve is single),
1332   // to fix it it is nesessary to have a possibility of adjusting seam
1333   // part of such curve either to left or right boundary of surface.
1334   // Test is performed only if flag AdjustOverDegen is not -1.
1335   // If AdjustOverDegen is True, seam part of curve is adjusted to
1336   // the left, and if False - to the right parametric boundary 
1337   // If treated case is detected, flag DONE4 is set to status
1338   // NOTE: currently, precision is Precision::PConfusion() since it 
1339   // is enough on encountered example
1340   // (ug_turbine-A.stp from ProSTEP Benchmark #3, entities ##2470 & 5680)
1341   // (r1001_ac.stp from Test Rally #10, face #35027 and others)
1342   if ( myAdjustOverDegen != -1 ) {
1343     if ( mySurf->IsUClosed(myPreci) ) {//#78 rln 12.03.99 S4135
1344       mySurf->IsDegenerated ( gp_Pnt(0,0,0), myPreci );  // pour calculer les dgnr
1345       if ( mySurf->NbSingularities(myPreci) > 0 ) { //rln S4135
1346         // 1st, find gap point (degenerated pole)
1347         Standard_Real PrevX=0.;
1348         Standard_Integer OnBound=0, PrevOnBound=0;
1349         Standard_Integer ind; // svv #1
1350         Standard_Boolean start = Standard_True;
1351         for ( ind=1; ind <= nbrPnt; ind++ ) {
1352           Standard_Real CurX = pnt2d(ind).X();
1353           // abv 16 Mar 00: trj3_s1-ug.stp #697: ignore points in singularity
1354           if ( mySurf->IsDegenerated ( points(ind), Precision::Confusion() ) )
1355             continue;
1356           OnBound = ( Abs ( Abs ( CurX - 0.5 * ( ul + uf ) ) - Up/2 ) <=
1357                      Precision::PConfusion() );
1358           if ( ! start &&  Abs ( Abs ( CurX - PrevX ) - Up/2 ) <= 0.01*Up ) 
1359             break;
1360           start = Standard_False;
1361           PrevX = CurX;
1362           PrevOnBound = OnBound;
1363         }
1364         // if found, adjust seam part
1365         if ( ind <= nbrPnt ) {
1366           PrevX = ( myAdjustOverDegen ? uf : ul );
1367           Standard_Real dU = Up/2 + Precision::PConfusion();
1368           if ( PrevOnBound ) {
1369             pnt2d(ind-1).SetX ( PrevX );
1370             for ( Standard_Integer j=ind-2; j >0; j-- ) {
1371               Standard_Real CurX = pnt2d(j).X();
1372               while ( CurX < PrevX - dU ) pnt2d(j).SetX ( CurX += Up );
1373               while ( CurX > PrevX + dU ) pnt2d(j).SetX ( CurX -= Up );
1374             }
1375           }
1376           else if ( OnBound ) {
1377             pnt2d(ind).SetX ( PrevX );
1378             for ( Standard_Integer j=ind+1; j <= nbrPnt; j++ ) {
1379               Standard_Real CurX = pnt2d(j).X();
1380               while ( CurX < PrevX - dU ) pnt2d(j).SetX ( CurX += Up );
1381               while ( CurX > PrevX + dU ) pnt2d(j).SetX ( CurX -= Up );
1382             }
1383           }
1384           myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE4);
1385         }
1386       }
1387     }
1388     else if ( mySurf->IsVClosed(myPreci) ) {//#78 rln 12.03.99 S4135
1389       mySurf->IsDegenerated ( gp_Pnt(0,0,0), myPreci );  // pour calculer les dgnr
1390       if ( mySurf->NbSingularities(myPreci) > 0 ) { //rln S4135
1391         // 1st, find gap point (degenerated pole)
1392         Standard_Real PrevY=0.;
1393         Standard_Integer OnBound=0, PrevOnBound=0;
1394         Standard_Integer ind; // svv #1
1395         Standard_Boolean start = Standard_True;
1396         for ( ind=1; ind <= nbrPnt; ind++ ) {
1397           Standard_Real CurY = pnt2d(ind).Y();
1398           // abv 16 Mar 00: trj3_s1-ug.stp #697: ignore points in singularity
1399           if ( mySurf->IsDegenerated ( points(ind), Precision::Confusion() ) )
1400             continue;
1401           OnBound = ( Abs ( Abs ( CurY - 0.5 * ( vl + vf ) ) - Vp/2 ) <=
1402                      Precision::PConfusion() );
1403           if ( ! start &&  Abs ( Abs ( CurY - PrevY ) - Vp/2 ) <= 0.01*Vp ) 
1404             break;
1405           start = Standard_False;
1406           PrevY = CurY;
1407           PrevOnBound = OnBound;
1408         }
1409         // if found, adjust seam part
1410         if ( ind <= nbrPnt ) {
1411           PrevY = ( myAdjustOverDegen ? vf : vl );
1412           Standard_Real dV = Vp/2 + Precision::PConfusion();
1413           if ( PrevOnBound ) {
1414             pnt2d(ind-1).SetY ( PrevY );
1415             for ( Standard_Integer j=ind-2; j >0; j-- ) {
1416               Standard_Real CurY = pnt2d(j).Y();
1417               while ( CurY < PrevY - dV ) pnt2d(j).SetY ( CurY += Vp );
1418               while ( CurY > PrevY + dV ) pnt2d(j).SetY ( CurY -= Vp );
1419             }
1420           }
1421           else if ( OnBound ) {
1422             pnt2d(ind).SetY ( PrevY );
1423             for ( Standard_Integer j=ind+1; j <= nbrPnt; j++ ) {
1424               Standard_Real CurY = pnt2d(j).Y();
1425               while ( CurY < PrevY - dV ) pnt2d(j).SetY ( CurY += Vp );
1426               while ( CurY > PrevY + dV ) pnt2d(j).SetY ( CurY -= Vp );
1427             }
1428           }
1429           myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE4);
1430         }
1431       }
1432     }
1433   }
1434
1435   //:q9: fill cashe
1436   myNbCashe = 2;
1437   if(ChangeCycle) {  // msv 10.08.04: avoid using of uninitialised field
1438   //if(myCashe3d[0].Distance(points(1))>Precision::Confusion() &&
1439   //   myCashe3d[1].Distance(points(1))>Precision::Confusion()) {
1440     myCashe3d[0] = points(1);
1441     myCashe3d[1] = points(nbrPnt);
1442     myCashe2d[0] = pnt2d(1);
1443     myCashe2d[1] = pnt2d(nbrPnt);
1444   }
1445   else {
1446     myCashe3d[1] = points(1);
1447     myCashe3d[0] = points(nbrPnt);
1448     myCashe2d[1] = pnt2d(1);
1449     myCashe2d[0] = pnt2d(nbrPnt);
1450   }
1451   return isDone;
1452 }
1453
1454 //=======================================================================
1455 //function : ApproximatePCurve
1456 //purpose  : 
1457 //=======================================================================
1458
1459 Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(const Standard_Integer /*nbrPnt*/,
1460                                                                              Handle(TColgp_HArray1OfPnt2d)& points2d, 
1461                                                                              Handle(TColStd_HArray1OfReal)& params,
1462                                                                              const Handle(Geom_Curve)& /*orig*/) const
1463 {
1464 //  Standard_Real resol = Min(mySurf->Adaptor3d()->VResolution(myPreci), mySurf->Adaptor3d()->UResolution(myPreci));
1465   Standard_Real theTolerance2d = myPreci; // (100*nbrPnt);//resol;
1466   Handle(Geom2d_Curve) C2d;
1467   try {
1468     OCC_CATCH_SIGNALS
1469     CheckPoints2d (points2d, params, theTolerance2d);
1470     Standard_Integer numberPnt = points2d->Length();
1471     
1472     TColgp_Array1OfPnt points3d(1,numberPnt);
1473     gp_Pnt2d pnt2d;
1474     gp_Pnt pnt;
1475     Standard_Integer i; // svv #1 
1476     for( i = 1; i <= numberPnt; i++) {
1477       pnt2d = points2d->Value(i);
1478       pnt.SetCoord(pnt2d.X(),pnt2d.Y(),0);
1479       points3d(i) = pnt;
1480     }
1481     
1482     GeomAPI_PointsToBSpline appr(points3d, params->Array1(), 1, 10, GeomAbs_C1, theTolerance2d);
1483     Handle(Geom_BSplineCurve) crv3d = appr.Curve();
1484     Standard_Integer NbPoles = crv3d->NbPoles();
1485     TColgp_Array1OfPnt poles3d (1, NbPoles);
1486     TColgp_Array1OfPnt2d poles2d (1, NbPoles);
1487     crv3d->Poles(poles3d);
1488     for( i = 1; i <= NbPoles; i++) {
1489       pnt2d.SetCoord(poles3d(i).X(),poles3d(i).Y());
1490       poles2d(i) = pnt2d;
1491     }
1492     TColStd_Array1OfReal weights (1,NbPoles);
1493     TColStd_Array1OfInteger multiplicities (1,crv3d->NbKnots());
1494     TColStd_Array1OfReal knots(1,crv3d->NbKnots());
1495     crv3d->Knots(knots);
1496     crv3d->Weights(weights);
1497     crv3d->Multiplicities(multiplicities);
1498     C2d = new Geom2d_BSplineCurve  ( poles2d, weights, knots, multiplicities, crv3d->Degree(), crv3d->IsPeriodic());
1499     return C2d;
1500   }
1501   catch(Standard_Failure const& anException) {
1502 #ifdef OCCT_DEBUG
1503  //:s5
1504     //    debug ...
1505     Standard_Integer nbp = params->Length();
1506     Standard_Integer nb2 = points2d->Length();
1507     cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(): Exception: ";
1508     anException.Print(cout); 
1509     cout<<"Pb Geom2dAPI_Approximate, tol2d="<<theTolerance2d<<" NbParams="<<nbp<<" NbPnts="<<nb2<<endl;
1510 //     if (nb2 > nbp) nb2 = nbp;
1511 //     Standard_Real rbp,rb2; rbp = nbp; rb2 = nb2;
1512 //     //    dbl.AddString ("NbP2d/NbParams puis  X Y Param -> mini");
1513 //     dbl.AddReals (rb2,rbp);
1514 //     for (Standard_Integer i = 1; i <= nb2; i ++) {
1515 //       gp_XYZ quoi (points2d->Value(i).X(),points2d->Value(i).Y(),params->Value(i) );
1516 //       dbl.AddXYZ (quoi);
1517 //     }
1518 #endif
1519      (void)anException;
1520      C2d.Nullify();
1521   }
1522   return C2d;
1523 }  
1524
1525 //=======================================================================
1526 //function : InterpolatePCurve
1527 //purpose  : 
1528 //=======================================================================
1529
1530  Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolatePCurve(const Standard_Integer nbrPnt,
1531                                                                               Handle(TColgp_HArray1OfPnt2d)& points2d, 
1532                                                                               Handle(TColStd_HArray1OfReal)& params,
1533                                                                               const Handle(Geom_Curve)& /*orig*/) const
1534 {
1535   Handle(Geom2d_Curve) C2d;    // NULL si echec
1536   Standard_Real theTolerance2d = myPreci / (100 * nbrPnt);
1537   try {
1538     OCC_CATCH_SIGNALS
1539     // on verifie d abord s il n y a pas de points confondus
1540     // si besoin on retouche les valeurs ...
1541     CheckPoints2d (points2d, params, theTolerance2d);
1542     Geom2dAPI_Interpolate myInterPol2d (points2d, params, 
1543                                         Standard_False, theTolerance2d);
1544     myInterPol2d.Perform();
1545     if (myInterPol2d.IsDone()) C2d = myInterPol2d.Curve();
1546   }
1547   catch(Standard_Failure const& anException) {
1548 #ifdef OCCT_DEBUG
1549 //:s5
1550 // //    debug ...
1551     Standard_Integer nbp = params->Length();
1552     Standard_Integer nb2 = points2d->Length();
1553     cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::InterpolatePCurve(): Exception: ";
1554     anException.Print(cout); 
1555     cout<<"Pb Geom2dAPI_Interpolate, tol2d="<<theTolerance2d<<" NbParams="<<nbp<<" NbPnts="<<nb2<<endl;
1556 //     if (nb2 > nbp) nb2 = nbp;
1557 //     Standard_Real rbp,rb2; rbp = nbp; rb2 = nb2;
1558 // //    dbl.AddString ("NbP2d/NbParams puis  X Y Param -> mini");
1559 //     dbl.AddReals (rb2,rbp);
1560 //     for (Standard_Integer i = 1; i <= nb2; i ++) {
1561 //       gp_XYZ quoi (points2d->Value(i).X(),points2d->Value(i).Y(),params->Value(i) );
1562 //       dbl.AddXYZ (quoi);
1563 //     }
1564 #endif
1565     (void)anException;
1566     C2d.Nullify();
1567   }
1568   return C2d;
1569 }
1570
1571 //=======================================================================
1572 //function : InterpolateCurve3d
1573 //purpose  : 
1574 //=======================================================================
1575
1576 Handle(Geom_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(const Standard_Integer,
1577                                                                             Handle(TColgp_HArray1OfPnt)& points, 
1578                                                                             Handle(TColStd_HArray1OfReal)& params,
1579                                                                             const Handle(Geom_Curve)& /*orig*/) const
1580 {
1581   Handle(Geom_Curve) C3d;    // NULL si echec
1582   try {
1583     OCC_CATCH_SIGNALS
1584     Standard_Real Tol = myPreci;
1585     CheckPoints(points, params, Tol);
1586     GeomAPI_Interpolate myInterPol(points, params, Standard_False, Tol);
1587     myInterPol.Perform();
1588     if (myInterPol.IsDone()) C3d = myInterPol.Curve();
1589   }
1590   catch(Standard_Failure const& anException) {
1591 #ifdef OCCT_DEBUG
1592  //:s5
1593     cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(): Exception: ";
1594     anException.Print(cout); cout << endl;
1595 #endif
1596     (void)anException;
1597     C3d.Nullify();
1598   }
1599   return C3d;
1600 }
1601
1602 //=======================================================================
1603 //function : CheckPoints
1604 //purpose  : 
1605 //=======================================================================
1606
1607  void ShapeConstruct_ProjectCurveOnSurface::CheckPoints(Handle(TColgp_HArray1OfPnt)& points,Handle(TColStd_HArray1OfReal)& params,Standard_Real& preci) const
1608 {
1609   Standard_Integer firstElem = points->Lower();
1610   Standard_Integer lastElem  = points->Upper();
1611   Standard_Integer i;
1612   Standard_Integer nbPntDropped = 0;
1613   Standard_Integer lastValid = firstElem; // indice of last undropped point
1614
1615   // will store 0 when the point is to be removed, 1 otherwise
1616   TColStd_Array1OfInteger tmpParam(firstElem, lastElem);
1617   for (i = firstElem; i<=lastElem ; i++) tmpParam.SetValue(i,1);
1618   Standard_Real DistMin2 = RealLast();
1619   gp_Pnt Prev = points->Value (lastValid);
1620   gp_Pnt Curr;
1621   for (i = firstElem + 1; i <= lastElem ; i ++) {
1622     Curr = points->Value(i);
1623     Standard_Real CurDist2 = Prev.SquareDistance(Curr);
1624     if (CurDist2 < gp::Resolution()) {  // test 0
1625       nbPntDropped ++;
1626       if ( i == lastElem ) tmpParam.SetValue(lastValid, 0); // last point kept
1627       else tmpParam.SetValue(i, 0);    // current dropped, lastValid unchanged
1628     } else {
1629       if (CurDist2 < DistMin2) 
1630         DistMin2 = CurDist2;
1631       // lastValid becomes the current (i.e. i)
1632       lastValid = i;
1633       Prev = Curr;
1634     }
1635   }
1636   if (DistMin2 < RealLast())
1637     preci = 0.9 * Sqrt (DistMin2); // preci est la distance min entre les points on la reduit un peu
1638   if (nbPntDropped == 0)
1639     return;
1640
1641 #ifdef OCCT_DEBUG
1642   cout << "Warning : removing 3d points for interpolation" << endl;
1643 #endif
1644   // Build new HArrays
1645   Standard_Integer newLast = lastElem - nbPntDropped;
1646   if ((newLast - firstElem + 1)  < 2) {
1647 #ifdef OCCT_DEBUG
1648     cout << "Too many degenerated points for 3D interpolation" << endl;
1649 #endif
1650     return;
1651   }
1652   Handle(TColgp_HArray1OfPnt) newPnts = 
1653     new TColgp_HArray1OfPnt(firstElem, newLast);
1654   Handle(TColStd_HArray1OfReal) newParams =
1655     new TColStd_HArray1OfReal(firstElem, newLast);
1656   Standard_Integer newCurr = 1;
1657   for (i = firstElem; i<= lastElem ; i++) {
1658     if (tmpParam.Value(i) == 1) { 
1659       newPnts->SetValue(newCurr, points->Value(i));
1660       newParams->SetValue(newCurr, params->Value(i));
1661       newCurr ++;
1662     }
1663   }
1664   points = newPnts;
1665   params = newParams;
1666   // on la reduit un peu
1667 }
1668
1669 //=======================================================================
1670 //function : CheckPoints2d
1671 //purpose  : 
1672 //=======================================================================
1673
1674  void ShapeConstruct_ProjectCurveOnSurface::CheckPoints2d(Handle(TColgp_HArray1OfPnt2d)& points,
1675                                                           Handle(TColStd_HArray1OfReal)& params,
1676                                                           Standard_Real& preci) const
1677 {
1678   Standard_Integer firstElem = points->Lower();
1679   Standard_Integer lastElem  = points->Upper();
1680   Standard_Integer i;
1681   Standard_Integer nbPntDropped = 0;
1682   Standard_Integer lastValid = firstElem; // indice of last undropped point
1683
1684   // will store 0 when the point is to be removed, 1 otherwise
1685   TColStd_Array1OfInteger tmpParam(firstElem, lastElem);
1686   for (i = firstElem; i<=lastElem ; i++) {
1687     tmpParam.SetValue(i,1);
1688   }
1689   Standard_Real DistMin2 = RealLast();
1690   gp_Pnt2d Prev = points->Value(lastValid);
1691   gp_Pnt2d Curr;
1692   for (i = firstElem + 1; i<=lastElem ; i++) {
1693     Curr = points->Value(i);
1694     Standard_Real CurDist2 = Prev.SquareDistance(Curr);
1695     if (CurDist2 < gp::Resolution()) {  // test 0
1696       nbPntDropped ++;
1697       if ( i == lastElem ) tmpParam.SetValue(lastValid, 0); // last point kept
1698       else tmpParam.SetValue(i, 0);    // current dropped, lastValid unchanged
1699     } else {
1700       if (CurDist2 < DistMin2) 
1701         DistMin2 = CurDist2;
1702       // lastValid becomes the current (i.e. i)
1703       lastValid = i;
1704       Prev = Curr;
1705     }
1706   }
1707   if (DistMin2 < RealLast())
1708     preci = 0.9 * Sqrt (DistMin2);
1709   if (nbPntDropped == 0)
1710     return;
1711
1712 #ifdef OCCT_DEBUG
1713   cout << "Warning : removing 2d points for interpolation" << endl;
1714 #endif
1715   // Build new HArrays
1716   Standard_Integer newLast = lastElem - nbPntDropped;
1717   if ((newLast - firstElem + 1)  < 2) {
1718 #ifdef OCCT_DEBUG
1719     cout << "Too many degenerated points for 2D interpolation" << endl;
1720 #endif
1721     //pdn 12.02.99 S4135 Creating pcurve with minimal length.
1722     tmpParam.SetValue(firstElem,1);
1723     tmpParam.SetValue(lastElem,1);
1724     gp_XY  lastPnt = points->Value(lastElem).XY();
1725     lastPnt.Add(gp_XY(preci,preci));
1726     points->SetValue(lastElem,lastPnt);
1727     newLast = firstElem+1;
1728     //return;
1729   }
1730   Handle(TColgp_HArray1OfPnt2d) newPnts = 
1731     new TColgp_HArray1OfPnt2d(firstElem, newLast);
1732   Handle(TColStd_HArray1OfReal) newParams =
1733     new TColStd_HArray1OfReal(firstElem, newLast);
1734   Standard_Integer newCurr = 1;
1735   for (i = firstElem; i <= lastElem ; i++) {
1736     if (tmpParam.Value(i) == 1) { 
1737 #ifdef OCCT_DEBUG
1738       cout << "Point " << i << " : " << points->Value(i).X() << " " << points->Value(i).Y() << " at param " <<  params->Value(i) << endl;
1739 #endif
1740       newPnts->SetValue(newCurr, points->Value(i));
1741       newParams->SetValue(newCurr, params->Value(i));
1742       newCurr ++;
1743     }
1744     else {
1745 #ifdef OCCT_DEBUG
1746       cout << "Removed " << i << " : " << points->Value(i).X() << " " << points->Value(i).Y() << " at param " <<  params->Value(i) << endl;
1747 #endif
1748     }
1749   }
1750   points = newPnts;
1751   params = newParams;
1752 }
1753
1754 //=======================================================================
1755 //function : IsAnIsoparametric
1756 //purpose  : 
1757 //=======================================================================
1758 //:S4030: modified for optimization
1759 //:p9 abv 11 Mar 99: PRO7226 #489490: find nearest boundary instead of first one
1760
1761  Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::IsAnIsoparametric(const Standard_Integer nbrPnt,
1762                                                                           const TColgp_Array1OfPnt& points,
1763                                                                           const TColStd_Array1OfReal& params,
1764                                                                           Standard_Boolean& isoTypeU,
1765                                                                           Standard_Boolean& p1OnIso,
1766                                                                           gp_Pnt2d& valueP1,
1767                                                                           Standard_Boolean& p2OnIso,
1768                                                                           gp_Pnt2d& valueP2,
1769                                                                           Standard_Boolean& isoPar2d3d,
1770                                                                           Handle(Geom_Curve)& cIso,
1771                                                                           Standard_Real& t1,
1772                                                                           Standard_Real& t2,
1773                                                                           TColStd_Array1OfReal& pout) const
1774 {
1775   try {    // RAJOUT
1776     OCC_CATCH_SIGNALS
1777     
1778   Standard_Real prec = Precision::Confusion();//myPreci;
1779     
1780   Standard_Boolean isoParam = Standard_False;
1781   isoPar2d3d = Standard_False;
1782   
1783   Standard_Real U1, U2, V1, V2;
1784   mySurf->Bounds(U1, U2, V1, V2);
1785   
1786   if ( mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
1787     Handle(Geom_RectangularTrimmedSurface) sTrim =
1788       Handle(Geom_RectangularTrimmedSurface)::DownCast(mySurf->Surface());
1789     sTrim->Bounds(U1, U2, V1, V2);
1790   }
1791   
1792   gp_Pnt pt;
1793   Standard_Integer mpt[2]; mpt[0] = mpt[1] = 0;
1794   Standard_Real t, tpar[2] = { 0.0, 0.0 }, isoValue=0.;
1795   Standard_Real mindist2;
1796   Standard_Real mind2[2];
1797   mindist2 = mind2[0] = mind2[1] = 4*prec*prec;
1798   
1799   p1OnIso = Standard_False;
1800   p2OnIso = Standard_False;
1801   const Bnd_Box* aBox = 0;
1802   
1803   for (Standard_Integer j=1; (j<=4) /*&& !isoParam*/; j++) {
1804     Standard_Real isoVal=0.;
1805     Standard_Boolean isoU=Standard_False; //szv#4:S4163:12Mar99 `isoU` must be Standard_Boolean
1806     Handle(Geom_Curve) cI;
1807     Standard_Real tt1, tt2;
1808     
1809     if      (j == 1 ) {
1810       if (Precision::IsInfinite(U1)) continue;
1811       cI = mySurf->UIso(U1);
1812       isoU = Standard_True;
1813       isoVal = U1;
1814       aBox = & mySurf->GetBoxUF();
1815     }
1816     else if (j == 2) {
1817       if (Precision::IsInfinite(U2)) continue;
1818       cI = mySurf->UIso(U2);
1819       isoU = Standard_True;
1820       isoVal = U2;
1821       aBox = & mySurf->GetBoxUL();
1822     }
1823     else if (j == 3) {
1824       if (Precision::IsInfinite(V1)) continue;
1825       cI = mySurf->VIso(V1);
1826       isoU = Standard_False;
1827       isoVal = V1;
1828       aBox = & mySurf->GetBoxVF();
1829     }
1830     else if (j == 4) {
1831       if (Precision::IsInfinite(V2)) continue;
1832       cI = mySurf->VIso(V2);
1833       isoU = Standard_False;
1834       isoVal = V2;
1835       aBox = & mySurf->GetBoxVL();
1836     }
1837     if(cI.IsNull())
1838       continue;
1839     
1840     if (isoU)  {  tt1 = V1;  tt2 = V2;  }
1841     else       {  tt1 = U1;  tt2 = U2;  }
1842
1843     gp_Pnt ext1, ext2;
1844     cI->D0(tt1, ext1);
1845     cI->D0(tt2, ext2);
1846
1847 // PATCH CKY 9-JUL-1998 : protection contre singularite
1848     gp_Pnt extmi;
1849     cI->D0( (tt1+tt2)/2,extmi);
1850     if (ext1.IsEqual(ext2,prec) && ext1.IsEqual(extmi,prec)) continue;
1851
1852     Standard_Boolean PtEQext1 = Standard_False;
1853     Standard_Boolean PtEQext2 = Standard_False;
1854
1855     Standard_Real currd2[2], tp[2] = {0, 0};
1856     Standard_Integer mp[2];
1857     
1858     for (Standard_Integer i=0; i<2; i++) {
1859       mp[i] = 0;
1860       Standard_Integer k = (i == 0 ? 1 : nbrPnt);
1861
1862       // si ext1 == ext2 => valueP1 == valueP2 => vect null plus tard
1863       currd2[i] = points(k).SquareDistance ( ext1 );
1864       if ( currd2[i] <= prec*prec && !PtEQext1) {
1865         mp[i] = 1;
1866         tp[i] = tt1;
1867         PtEQext1 = Standard_True;
1868         continue;
1869       } 
1870       
1871       currd2[i] = points(k).SquareDistance ( ext2 );
1872       if ( currd2[i] <= prec*prec && !PtEQext2) {
1873         mp[i] = 2;
1874         tp[i] = tt2;
1875         PtEQext2 = Standard_True;
1876         continue;
1877       }  
1878       
1879       // On evite de projecter sur un iso degenere
1880       // on doit egalement le faire pour l apex du cone
1881       if (mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) && !isoU) {
1882         continue;
1883       }
1884       
1885       if(aBox->IsOut(points(k))) continue;
1886       
1887       Standard_Real Cf = cI->FirstParameter();
1888       Standard_Real Cl = cI->LastParameter();
1889       if (Precision::IsInfinite(Cf))  Cf = -1000;
1890       if (Precision::IsInfinite(Cl))  Cl = +1000;
1891
1892       ShapeAnalysis_Curve sac;
1893       Standard_Real dist = sac.Project (cI,points(k),prec,pt,t,Cf,Cl);
1894       currd2[i] = dist * dist;
1895       if ((dist <= prec) && (t>= Cf) && (t<=Cl)) {
1896         mp[i]  = 3;
1897         tp[i] = t;
1898       }
1899     }
1900
1901     //:e7 abv 21 Apr 98: ProSTEP TR8, r0501_pe #56679:
1902     // avoid possible null-length curves
1903     if ( mp[0] >0 && mp[1] >0 &&
1904          Abs ( tp[0] - tp[1] ) < Precision::PConfusion() ) continue;
1905     
1906     
1907     if (mp[0] > 0 && 
1908         ( ! p1OnIso || currd2[0] < mind2[0] ) ) {
1909       p1OnIso = Standard_True;
1910       mind2[0] = currd2[0]; // LP2.stp #105899: FLT_INVALID_OPERATION on Windows 7 VC 9 Release mode on the whole file 
1911       if (isoU) valueP1.SetCoord(isoVal, tp[0]);
1912       else      valueP1.SetCoord(tp[0], isoVal);
1913     }
1914
1915     if (mp[1] > 0 &&
1916         ( ! p2OnIso || currd2[1] < mind2[1] ) ) {
1917       p2OnIso = Standard_True;
1918       mind2[1] = currd2[1];
1919       if (isoU) valueP2.SetCoord(isoVal, tp[1]);
1920       else      valueP2.SetCoord(tp[1], isoVal);
1921     }
1922
1923     if ( mp[0] <=0 || mp[1] <=0 ) continue;
1924
1925     Standard_Real md2 = currd2[0] + currd2[1];
1926     if ( mindist2 <= md2 ) continue;
1927     
1928     mindist2 = md2;
1929     mpt[0] = mp[0];
1930     mpt[1] = mp[1];
1931     tpar[0] = tp[0];
1932     tpar[1] = tp[1];
1933     isoTypeU = isoU;
1934     isoValue = isoVal;
1935     cIso = cI;
1936     t1 = tt1;
1937     t2 = tt2;
1938   }
1939     
1940   // probablely it concerns an isoparametrics
1941   if ( mpt[0] >0 && mpt[1] >0 ) {
1942
1943     p1OnIso = p2OnIso = Standard_True;
1944     if (isoTypeU) {
1945       valueP1.SetCoord(isoValue, tpar[0]);
1946       valueP2.SetCoord(isoValue, tpar[1]);
1947     }
1948     else {
1949       valueP1.SetCoord(tpar[0], isoValue);
1950       valueP2.SetCoord(tpar[1], isoValue);
1951     }
1952
1953     if ( mpt[0] != 3 && mpt[1] != 3 ) {
1954       isoPar2d3d = Standard_True;
1955       for (Standard_Integer i=2; i < nbrPnt && isoPar2d3d; i++){
1956         if (tpar[1] > tpar[0])  t = params(i);
1957         else                    t = t1+t2-params(i);
1958         cIso->D0(t, pt);
1959         if (!points(i).IsEqual(pt, prec)) isoPar2d3d = Standard_False;
1960       }
1961     }
1962
1963     if (isoPar2d3d) isoParam = Standard_True;
1964     else {
1965       Standard_Real prevParam = tpar[0];
1966       Standard_Real Cf, Cl;
1967       Standard_Boolean isoByDistance = Standard_True;
1968       Cf = cIso->FirstParameter();
1969       Cl = cIso->LastParameter();
1970       if (Precision::IsInfinite(Cf))  Cf = -1000;
1971       if (Precision::IsInfinite(Cl))  Cl = +1000;
1972         
1973       ShapeAnalysis_Curve sac;
1974       for (Standard_Integer i=2; i < nbrPnt && isoByDistance; i++) {
1975         Standard_Real dist = sac.NextProject (prevParam,cIso,points(i),
1976                                               prec,pt,t,Cf,Cl,
1977                                               Standard_False); //:j8 abv 10.12.98: TR10 r0501_db.stp #9423: avoid adjusting to ends
1978         prevParam = t;
1979         pout(i)=t;
1980         if( (dist > prec) || (t < Cf) || (t > Cl) ) 
1981           isoByDistance = Standard_False;
1982       }
1983       if (isoByDistance) isoParam = Standard_True;
1984     }
1985   }
1986 /*  if (!isoParam) {    CKY 29-mai-1997 : garder tout ce qu on peut ?
1987     p1OnIso = Standard_False;
1988     p2OnIso = Standard_False;
1989   }  */
1990   return isoParam;
1991   }  // RAJOUT
1992   catch(Standard_Failure const& anException) {
1993 #ifdef OCCT_DEBUG
1994 //  pb : on affiche ce qu on peut
1995     for (Standard_Integer numpnt = 1; numpnt <= nbrPnt; numpnt ++) {
1996       cout<<"["<<numpnt<<"]param="<<params(numpnt)<<" point=("<<
1997         points(numpnt).X()<<"  "<<points(numpnt).Y()<<"  "<<points(numpnt).Z()<<")"<<endl;
1998     }
1999     cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::IsAnIsoparametric(): Exception: ";
2000     anException.Print(cout); cout << endl;
2001 #endif
2002     (void)anException;
2003     return Standard_False;
2004   }
2005 }
2006
2007 /* S4135 : BestExtremum is commented after IsAnIsoparametric works with Precision::Confusion()
2008 //=======================================================================
2009 //function : BestExtremum
2010 //purpose  : auxiliaire prenant le meilleur extremum si ISO car doute possible
2011 //=======================================================================
2012
2013  gp_Pnt2d ShapeConstruct_ProjectCurveOnSurface::BestExtremum(const gp_Pnt2d& P2iso,const gp_Pnt& P3ext,const gp_Pnt& P3next) const
2014 {
2015 //  P2iso a ete calcule depuis P3ext sur une iso externe de la surface
2016 //  En principe bon mais circularite possible ... et IsU/VClosed faillible
2017 //    (si baillement 1e-4 ou 1e-5, on est dedans !). DONC
2018 //  1/ on privilegie l iso mais a tout hasard on verifie si Surf meilleur
2019 //  2/ si iso, attention a la circularite (cas limite)
2020
2021 //  NB : si isoParam, on suppose que P2iso est bon (car il y en a 2). A voir...
2022
2023 //  D abord, calcul p2ext depuis la surface. choix surface/iso
2024   return P2iso;
2025   Standard_Real prec = Precision::Confusion();//myPreci;
2026   gp_Pnt2d P2cal = mySurf->ValueOfUV(P3ext, prec);
2027   gp_Pnt   P3cal = mySurf->Value (P2cal);
2028   Standard_Real dcal = P3ext.Distance (P3cal);
2029   Standard_Real dnxt = P3ext.Distance (P3next);
2030   if (dcal > dnxt) return P2iso;    // en fait protection sur BUG (PRO8468)
2031
2032 //  On choisit entre P2iso et P2cal, le plus proche de P2next ... !!!
2033   gp_Pnt2d P2next = mySurf->ValueOfUV(P3next, prec);
2034   if (P2next.Distance(P2cal) < P2next.Distance(P2iso)) return P2cal;
2035   return P2iso;
2036 }
2037 */