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