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