3f9f2ca388575982888a2f432e02b612a5a54974
[occt.git] / src / ProjLib / ProjLib_CompProjectedCurve.cxx
1 // Created on: 1997-09-23
2 // Created by: Roman BORISOV
3 // Copyright (c) 1997-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17
18 #include <Adaptor2d_HCurve2d.hxx>
19 #include <Adaptor3d_HCurve.hxx>
20 #include <Adaptor3d_HSurface.hxx>
21 #include <Extrema_ExtCS.hxx>
22 #include <Extrema_ExtPS.hxx>
23 #include <Extrema_GenLocateExtPS.hxx>
24 #include <Extrema_POnCurv.hxx>
25 #include <Extrema_POnSurf.hxx>
26 #include <GeomAbs_CurveType.hxx>
27 #include <GeomLib.hxx>
28 #include <gp_Mat2d.hxx>
29 #include <gp_Pnt2d.hxx>
30 #include <gp_Vec2d.hxx>
31 #include <gp_XY.hxx>
32 #include <Precision.hxx>
33 #include <ProjLib_CompProjectedCurve.hxx>
34 #include <ProjLib_HCompProjectedCurve.hxx>
35 #include <ProjLib_PrjResolve.hxx>
36 #include <Standard_DomainError.hxx>
37 #include <Standard_NoSuchObject.hxx>
38 #include <Standard_NotImplemented.hxx>
39 #include <Standard_OutOfRange.hxx>
40 #include <TColgp_HSequenceOfPnt.hxx>
41
42 #define FuncTol 1.e-10
43
44 #ifdef OCCT_DEBUG_CHRONO
45 #include <OSD_Timer.hxx>
46
47 static OSD_Chronometer chr_init_point, chr_dicho_bound;
48
49 Standard_EXPORT Standard_Real t_init_point, t_dicho_bound;
50 Standard_EXPORT Standard_Integer init_point_count, dicho_bound_count;
51
52 static void InitChron(OSD_Chronometer& ch)
53
54   ch.Reset();
55   ch.Start();
56 }
57
58 static void ResultChron( OSD_Chronometer & ch, Standard_Real & time) 
59 {
60   Standard_Real tch ;
61   ch.Stop();
62   ch.Show(tch);
63   time=time +tch;
64 }
65 #endif
66
67
68 //=======================================================================
69 //function : d1
70 //purpose  : computes first derivative of the projected curve
71 //=======================================================================
72
73 static void d1(const Standard_Real t,
74   const Standard_Real u,
75   const Standard_Real v,
76   gp_Vec2d& V, 
77   const Handle(Adaptor3d_HCurve)& Curve, 
78   const Handle(Adaptor3d_HSurface)& Surface)
79 {
80   gp_Pnt S, C;
81   gp_Vec DS1_u, DS1_v, DS2_u, DS2_uv, DS2_v, DC1_t;
82   Surface->D2(u, v, S, DS1_u, DS1_v, DS2_u, DS2_v, DS2_uv);
83   Curve->D1(t, C, DC1_t);
84   gp_Vec Ort(C, S);// Ort = S - C
85
86   gp_Vec2d dE_dt(-DC1_t*DS1_u, -DC1_t*DS1_v);
87   gp_XY dE_du(DS1_u*DS1_u + Ort*DS2_u, 
88     DS1_u*DS1_v + Ort*DS2_uv);
89   gp_XY dE_dv(DS1_v*DS1_u + Ort*DS2_uv, 
90     DS1_v*DS1_v + Ort*DS2_v);
91
92   Standard_Real det = dE_du.X()*dE_dv.Y() - dE_du.Y()*dE_dv.X();
93   if (fabs(det) < gp::Resolution()) Standard_ConstructionError::Raise();
94
95   gp_Mat2d M(gp_XY(dE_dv.Y()/det, -dE_du.Y()/det), 
96     gp_XY(-dE_dv.X()/det, dE_du.X()/det));
97
98   V = - gp_Vec2d(gp_Vec2d(M.Row(1))*dE_dt, gp_Vec2d(M.Row(2))*dE_dt);
99 }
100
101 //=======================================================================
102 //function : d2
103 //purpose  : computes second derivative of the projected curve
104 //=======================================================================
105
106 static void d2(const Standard_Real t,
107   const Standard_Real u,
108   const Standard_Real v,
109   gp_Vec2d& V1, gp_Vec2d& V2,
110   const Handle(Adaptor3d_HCurve)& Curve, 
111   const Handle(Adaptor3d_HSurface)& Surface)
112 {
113   gp_Pnt S, C;
114   gp_Vec DS1_u, DS1_v, DS2_u, DS2_uv, DS2_v, 
115     DS3_u, DS3_v, DS3_uuv, DS3_uvv, 
116     DC1_t, DC2_t;
117   Surface->D3(u, v, S, DS1_u, DS1_v, DS2_u, DS2_v, DS2_uv, 
118     DS3_u, DS3_v, DS3_uuv, DS3_uvv);
119   Curve->D2(t, C, DC1_t, DC2_t);
120   gp_Vec Ort(C, S);
121
122   gp_Vec2d dE_dt(-DC1_t*DS1_u, -DC1_t*DS1_v);
123   gp_XY dE_du(DS1_u*DS1_u + Ort*DS2_u, 
124     DS1_u*DS1_v + Ort*DS2_uv);
125   gp_XY dE_dv(DS1_v*DS1_u + Ort*DS2_uv, 
126     DS1_v*DS1_v + Ort*DS2_v);
127
128   Standard_Real det = dE_du.X()*dE_dv.Y() - dE_du.Y()*dE_dv.X();
129   if (fabs(det) < gp::Resolution()) Standard_ConstructionError::Raise();
130
131   gp_Mat2d M(gp_XY(dE_dv.Y()/det, -dE_du.Y()/det), 
132     gp_XY(-dE_dv.X()/det, dE_du.X()/det));
133
134   // First derivative
135   V1 = - gp_Vec2d(gp_Vec2d(M.Row(1))*dE_dt, gp_Vec2d(M.Row(2))*dE_dt);
136
137   /* Second derivative */
138
139   // Computation of d2E_dt2 = S1
140   gp_Vec2d d2E_dt(-DC2_t*DS1_u, -DC2_t*DS1_v);
141
142   // Computation of 2*(d2E/dtdX)(dX/dt) = S2
143   gp_Vec2d d2E1_dtdX(-DC1_t*DS2_u,
144     -DC1_t*DS2_uv);
145   gp_Vec2d d2E2_dtdX(-DC1_t*DS2_uv,
146     -DC1_t*DS2_v);
147   gp_Vec2d S2 = 2*gp_Vec2d(d2E1_dtdX*V1, d2E2_dtdX*V1);
148
149   // Computation of (d2E/dX2)*(dX/dt)2 = S3
150
151   // Row11 = (d2E1/du2, d2E1/dudv)
152   Standard_Real tmp;
153   gp_Vec2d Row11(3*DS1_u*DS2_u + Ort*DS3_u,
154     tmp = 2*DS1_u*DS2_uv + 
155     DS1_v*DS2_u + Ort*DS3_uuv);  
156
157   // Row12 = (d2E1/dudv, d2E1/dv2)
158   gp_Vec2d Row12(tmp, DS2_v*DS1_u + 2*DS1_v*DS2_uv + 
159     Ort*DS3_uvv);
160
161   // Row21 = (d2E2/du2, d2E2/dudv)
162   gp_Vec2d Row21(DS2_u*DS1_v + 2*DS1_u*DS2_uv + Ort*DS3_uuv, 
163     tmp = 2*DS2_uv*DS1_v + DS1_u*DS2_v + Ort*DS3_uvv);
164
165   // Row22 = (d2E2/duv, d2E2/dvdv)
166   gp_Vec2d Row22(tmp, 3*DS1_v*DS2_v + Ort*DS3_v);
167
168   gp_Vec2d S3(V1*gp_Vec2d(Row11*V1, Row12*V1),
169     V1*gp_Vec2d(Row21*V1, Row22*V1));
170
171   gp_Vec2d Sum = d2E_dt + S2 + S3;
172
173   V2 = - gp_Vec2d(gp_Vec2d(M.Row(1))*Sum, gp_Vec2d(M.Row(2))*Sum);
174 }
175 //=======================================================================
176 //function : d1CurveOnSurf
177 //purpose  : computes first derivative of the 3d projected curve
178 //=======================================================================
179
180 #if 0
181 static void d1CurvOnSurf(const Standard_Real t,
182   const Standard_Real u,
183   const Standard_Real v,
184   gp_Vec& V, 
185   const Handle(Adaptor3d_HCurve)& Curve, 
186   const Handle(Adaptor3d_HSurface)& Surface)
187 {
188   gp_Pnt S, C;
189   gp_Vec2d V2d;
190   gp_Vec DS1_u, DS1_v, DS2_u, DS2_uv, DS2_v, DC1_t;
191   Surface->D2(u, v, S, DS1_u, DS1_v, DS2_u, DS2_v, DS2_uv);
192   Curve->D1(t, C, DC1_t);
193   gp_Vec Ort(C, S);// Ort = S - C
194
195   gp_Vec2d dE_dt(-DC1_t*DS1_u, -DC1_t*DS1_v);
196   gp_XY dE_du(DS1_u*DS1_u + Ort*DS2_u, 
197     DS1_u*DS1_v + Ort*DS2_uv);
198   gp_XY dE_dv(DS1_v*DS1_u + Ort*DS2_uv, 
199     DS1_v*DS1_v + Ort*DS2_v);
200
201   Standard_Real det = dE_du.X()*dE_dv.Y() - dE_du.Y()*dE_dv.X();
202   if (fabs(det) < gp::Resolution()) Standard_ConstructionError::Raise();
203
204   gp_Mat2d M(gp_XY(dE_dv.Y()/det, -dE_du.Y()/det), 
205     gp_XY(-dE_dv.X()/det, dE_du.X()/det));
206
207   V2d = - gp_Vec2d(gp_Vec2d(M.Row(1))*dE_dt, gp_Vec2d(M.Row(2))*dE_dt);
208
209   V = DS1_u * V2d.X() + DS1_v * V2d.Y();
210
211 }
212 #endif
213
214 //=======================================================================
215 //function : d2CurveOnSurf
216 //purpose  : computes second derivative of the 3D projected curve
217 //=======================================================================
218
219 static void d2CurvOnSurf(const Standard_Real t,
220   const Standard_Real u,
221   const Standard_Real v,
222   gp_Vec& V1 , gp_Vec& V2 ,
223   const Handle(Adaptor3d_HCurve)& Curve, 
224   const Handle(Adaptor3d_HSurface)& Surface)
225 {
226   gp_Pnt S, C;
227   gp_Vec2d V12d,V22d;
228   gp_Vec DS1_u, DS1_v, DS2_u, DS2_uv, DS2_v, 
229     DS3_u, DS3_v, DS3_uuv, DS3_uvv, 
230     DC1_t, DC2_t;
231   Surface->D3(u, v, S, DS1_u, DS1_v, DS2_u, DS2_v, DS2_uv, 
232     DS3_u, DS3_v, DS3_uuv, DS3_uvv);
233   Curve->D2(t, C, DC1_t, DC2_t);
234   gp_Vec Ort(C, S);
235
236   gp_Vec2d dE_dt(-DC1_t*DS1_u, -DC1_t*DS1_v);
237   gp_XY dE_du(DS1_u*DS1_u + Ort*DS2_u, 
238     DS1_u*DS1_v + Ort*DS2_uv);
239   gp_XY dE_dv(DS1_v*DS1_u + Ort*DS2_uv, 
240     DS1_v*DS1_v + Ort*DS2_v);
241
242   Standard_Real det = dE_du.X()*dE_dv.Y() - dE_du.Y()*dE_dv.X();
243   if (fabs(det) < gp::Resolution()) Standard_ConstructionError::Raise();
244
245   gp_Mat2d M(gp_XY(dE_dv.Y()/det, -dE_du.Y()/det), 
246     gp_XY(-dE_dv.X()/det, dE_du.X()/det));
247
248   // First derivative
249   V12d = - gp_Vec2d(gp_Vec2d(M.Row(1))*dE_dt, gp_Vec2d(M.Row(2))*dE_dt);
250
251   /* Second derivative */
252
253   // Computation of d2E_dt2 = S1
254   gp_Vec2d d2E_dt(-DC2_t*DS1_u, -DC2_t*DS1_v);
255
256   // Computation of 2*(d2E/dtdX)(dX/dt) = S2
257   gp_Vec2d d2E1_dtdX(-DC1_t*DS2_u,
258     -DC1_t*DS2_uv);
259   gp_Vec2d d2E2_dtdX(-DC1_t*DS2_uv,
260     -DC1_t*DS2_v);
261   gp_Vec2d S2 = 2*gp_Vec2d(d2E1_dtdX*V12d, d2E2_dtdX*V12d);
262
263   // Computation of (d2E/dX2)*(dX/dt)2 = S3
264
265   // Row11 = (d2E1/du2, d2E1/dudv)
266   Standard_Real tmp;
267   gp_Vec2d Row11(3*DS1_u*DS2_u + Ort*DS3_u,
268     tmp = 2*DS1_u*DS2_uv + 
269     DS1_v*DS2_u + Ort*DS3_uuv);  
270
271   // Row12 = (d2E1/dudv, d2E1/dv2)
272   gp_Vec2d Row12(tmp, DS2_v*DS1_u + 2*DS1_v*DS2_uv + 
273     Ort*DS3_uvv);
274
275   // Row21 = (d2E2/du2, d2E2/dudv)
276   gp_Vec2d Row21(DS2_u*DS1_v + 2*DS1_u*DS2_uv + Ort*DS3_uuv, 
277     tmp = 2*DS2_uv*DS1_v + DS1_u*DS2_v + Ort*DS3_uvv);
278
279   // Row22 = (d2E2/duv, d2E2/dvdv)
280   gp_Vec2d Row22(tmp, 3*DS1_v*DS2_v + Ort*DS3_v);
281
282   gp_Vec2d S3(V12d*gp_Vec2d(Row11*V12d, Row12*V12d),
283     V12d*gp_Vec2d(Row21*V12d, Row22*V12d));
284
285   gp_Vec2d Sum = d2E_dt + S2 + S3;
286
287   V22d = - gp_Vec2d(gp_Vec2d(M.Row(1))*Sum, gp_Vec2d(M.Row(2))*Sum);
288
289   V1 = DS1_u * V12d.X() + DS1_v * V12d.Y();
290   V2 =     DS2_u * V12d.X() *V12d.X()  
291     +     DS1_u * V22d.X() 
292     + 2 * DS2_uv * V12d.X() *V12d.Y()
293     +     DS2_v * V12d.Y() * V12d.Y()
294     +     DS1_v * V22d.Y();
295 }
296
297 //=======================================================================
298 //function : ExactBound
299 //purpose  : computes exact boundary point
300 //=======================================================================
301
302 static Standard_Boolean ExactBound(gp_Pnt& Sol, 
303   const Standard_Real NotSol, 
304   const Standard_Real Tol, 
305   const Standard_Real TolU, 
306   const Standard_Real TolV,  
307   const Handle(Adaptor3d_HCurve)& Curve, 
308   const Handle(Adaptor3d_HSurface)& Surface)
309 {
310   Standard_Real U0, V0, t, t1, t2, FirstU, LastU, FirstV, LastV;
311   gp_Pnt2d POnS;
312   U0 = Sol.Y();
313   V0 = Sol.Z();
314   FirstU = Surface->FirstUParameter();
315   LastU = Surface->LastUParameter();
316   FirstV = Surface->FirstVParameter();
317   LastV = Surface->LastVParameter();
318   // Here we have to compute the boundary that projection is going to intersect
319   gp_Vec2d D2d;
320   //these variables are to estimate which boundary has more apportunity 
321   //to be intersected
322   Standard_Real RU1, RU2, RV1, RV2; 
323   d1(Sol.X(), U0, V0, D2d, Curve, Surface);
324   // Here we assume that D2d != (0, 0)
325   if(Abs(D2d.X()) < gp::Resolution()) 
326   {
327     RU1 = Precision::Infinite();
328     RU2 = Precision::Infinite();
329     RV1 = V0 - FirstV;
330     RV2 = LastV - V0;
331   }
332   else if(Abs(D2d.Y()) < gp::Resolution()) 
333   {
334     RU1 = U0 - FirstU;
335     RU2 = LastU - U0;
336     RV1 = Precision::Infinite();
337     RV2 = Precision::Infinite();    
338   }
339   else 
340   {
341     RU1 = gp_Pnt2d(U0, V0).
342       Distance(gp_Pnt2d(FirstU, V0 + (FirstU - U0)*D2d.Y()/D2d.X()));
343     RU2 = gp_Pnt2d(U0, V0).
344       Distance(gp_Pnt2d(LastU, V0 + (LastU - U0)*D2d.Y()/D2d.X()));
345     RV1 = gp_Pnt2d(U0, V0).
346       Distance(gp_Pnt2d(U0 + (FirstV - V0)*D2d.X()/D2d.Y(), FirstV));
347     RV2 = gp_Pnt2d(U0, V0).
348       Distance(gp_Pnt2d(U0 + (LastV - V0)*D2d.X()/D2d.Y(), LastV));
349   }
350   TColgp_SequenceOfPnt Seq;
351   Seq.Append(gp_Pnt(FirstU, RU1, 2));
352   Seq.Append(gp_Pnt(LastU, RU2, 2));
353   Seq.Append(gp_Pnt(FirstV, RV1, 3));
354   Seq.Append(gp_Pnt(LastV, RV2, 3));
355   Standard_Integer i, j;
356   for(i = 1; i <= 3; i++)
357     for(j = 1; j <= 4-i; j++)
358       if(Seq(j).Y() < Seq(j+1).Y()) 
359       {
360         gp_Pnt swp;
361         swp = Seq.Value(j+1);
362         Seq.ChangeValue(j+1) = Seq.Value(j);
363         Seq.ChangeValue(j) = swp;
364       }
365
366       t = Sol.X();
367       t1 = Min(Sol.X(), NotSol);
368       t2 = Max(Sol.X(), NotSol);
369
370       Standard_Boolean isDone = Standard_False;
371       while (!Seq.IsEmpty()) 
372       {
373         gp_Pnt P;
374         P = Seq.Last();
375         Seq.Remove(Seq.Length());
376         ProjLib_PrjResolve aPrjPS(Curve->Curve(), 
377           Surface->Surface(), 
378           Standard_Integer(P.Z()));
379         if(Standard_Integer(P.Z()) == 2) 
380         {
381           aPrjPS.Perform(t, P.X(), V0, gp_Pnt2d(Tol, TolV), 
382             gp_Pnt2d(t1, Surface->FirstVParameter()), 
383             gp_Pnt2d(t2, Surface->LastVParameter()), FuncTol);
384           if(!aPrjPS.IsDone()) continue;
385           POnS = aPrjPS.Solution();
386           Sol = gp_Pnt(POnS.X(), P.X(), POnS.Y());
387           isDone = Standard_True;
388           break;
389         }
390         else 
391         {
392           aPrjPS.Perform(t, U0, P.X(), gp_Pnt2d(Tol, TolU), 
393             gp_Pnt2d(t1, Surface->FirstUParameter()), 
394             gp_Pnt2d(t2, Surface->LastUParameter()), FuncTol);
395           if(!aPrjPS.IsDone()) continue;
396           POnS = aPrjPS.Solution();
397           Sol = gp_Pnt(POnS.X(), POnS.Y(), P.X());
398           isDone = Standard_True;
399           break;
400         }
401       }
402
403       return isDone;
404 }
405
406 //=======================================================================
407 //function : DichExactBound
408 //purpose  : computes exact boundary point
409 //=======================================================================
410
411 static void DichExactBound(gp_Pnt& Sol, 
412   const Standard_Real NotSol, 
413   const Standard_Real Tol, 
414   const Standard_Real TolU, 
415   const Standard_Real TolV,  
416   const Handle(Adaptor3d_HCurve)& Curve, 
417   const Handle(Adaptor3d_HSurface)& Surface)
418 {
419 #ifdef OCCT_DEBUG_CHRONO
420   InitChron(chr_dicho_bound);
421 #endif
422
423   Standard_Real U0, V0, t;
424   gp_Pnt2d POnS;
425   U0 = Sol.Y();
426   V0 = Sol.Z();
427   ProjLib_PrjResolve aPrjPS(Curve->Curve(), Surface->Surface(), 1);
428
429   Standard_Real aNotSol = NotSol;
430   while (fabs(Sol.X() - aNotSol) > Tol) 
431   {
432     t = (Sol.X() + aNotSol)/2;
433     aPrjPS.Perform(t, U0, V0, gp_Pnt2d(TolU, TolV), 
434       gp_Pnt2d(Surface->FirstUParameter(),Surface->FirstVParameter()), 
435       gp_Pnt2d(Surface->LastUParameter(),Surface->LastVParameter()), 
436       FuncTol, Standard_True);
437
438     if (aPrjPS.IsDone()) 
439     {
440       POnS = aPrjPS.Solution();
441       Sol = gp_Pnt(t, POnS.X(), POnS.Y());
442       U0=Sol.Y();
443       V0=Sol.Z();
444     }
445     else aNotSol = t; 
446   }
447 #ifdef OCCT_DEBUG_CHRONO
448   ResultChron(chr_dicho_bound,t_dicho_bound);
449   dicho_bound_count++;
450 #endif
451 }
452
453 //=======================================================================
454 //function : InitialPoint
455 //purpose  : 
456 //=======================================================================
457
458 static Standard_Boolean InitialPoint(const gp_Pnt& Point, 
459   const Standard_Real t,
460   const Handle(Adaptor3d_HCurve)& C,
461   const Handle(Adaptor3d_HSurface)& S, 
462   const Standard_Real TolU, 
463   const Standard_Real TolV, 
464   Standard_Real& U, 
465   Standard_Real& V)
466 {
467
468   ProjLib_PrjResolve aPrjPS(C->Curve(), S->Surface(), 1);
469   Standard_Real ParU,ParV;
470   Extrema_ExtPS aExtPS;
471   aExtPS.Initialize(S->Surface(), S->FirstUParameter(), 
472     S->LastUParameter(), S->FirstVParameter(), 
473     S->LastVParameter(), TolU, TolV);
474
475   aExtPS.Perform(Point);
476   Standard_Integer argmin = 0;
477   if (aExtPS.IsDone() && aExtPS.NbExt()) 
478   {
479     Standard_Integer i, Nend;
480     // Search for the nearest solution which is also a normal projection
481     Nend = aExtPS.NbExt();
482     for(i = 1; i <= Nend; i++)
483     {
484       Extrema_POnSurf POnS = aExtPS.Point(i);
485       POnS.Parameter(ParU, ParV);
486       aPrjPS.Perform(t, ParU, ParV, gp_Pnt2d(TolU, TolV), 
487         gp_Pnt2d(S->FirstUParameter(), S->FirstVParameter()), 
488         gp_Pnt2d(S->LastUParameter(), S->LastVParameter()), 
489         FuncTol, Standard_True);
490       if(aPrjPS.IsDone() )
491         if  (argmin == 0 || aExtPS.SquareDistance(i) < aExtPS.SquareDistance(argmin)) argmin = i;  
492     }
493   }
494   if( argmin == 0 ) return Standard_False;
495   else
496   {  
497     Extrema_POnSurf POnS = aExtPS.Point(argmin);
498     POnS.Parameter(U, V);
499     return Standard_True;
500   }
501 }
502
503 //=======================================================================
504 //function : ProjLib_CompProjectedCurve
505 //purpose  : 
506 //=======================================================================
507
508 ProjLib_CompProjectedCurve::ProjLib_CompProjectedCurve()
509 : myNbCurves(0),
510   myTolU    (0.0),
511   myTolV    (0.0),
512   myMaxDist (0.0)
513 {
514 }
515
516 //=======================================================================
517 //function : ProjLib_CompProjectedCurve
518 //purpose  : 
519 //=======================================================================
520
521 ProjLib_CompProjectedCurve::ProjLib_CompProjectedCurve
522                            (const Handle(Adaptor3d_HSurface)& theSurface,
523                             const Handle(Adaptor3d_HCurve)&   theCurve,
524                             const Standard_Real               theTolU,
525                             const Standard_Real               theTolV)
526 : mySurface (theSurface),
527   myCurve   (theCurve),
528   myNbCurves(0),
529   mySequence(new ProjLib_HSequenceOfHSequenceOfPnt()),
530   myTolU    (theTolU),
531   myTolV    (theTolV),
532   myMaxDist (-1.0)
533 {
534   Init();
535 }
536
537 //=======================================================================
538 //function : ProjLib_CompProjectedCurve
539 //purpose  : 
540 //=======================================================================
541
542 ProjLib_CompProjectedCurve::ProjLib_CompProjectedCurve
543                            (const Handle(Adaptor3d_HSurface)& theSurface,
544                             const Handle(Adaptor3d_HCurve)&   theCurve,
545                             const Standard_Real               theTolU,
546                             const Standard_Real               theTolV,
547                             const Standard_Real               theMaxDist)
548 : mySurface (theSurface),
549   myCurve   (theCurve),
550   myNbCurves(0),
551   mySequence(new ProjLib_HSequenceOfHSequenceOfPnt()),
552   myTolU    (theTolU),
553   myTolV    (theTolV),
554   myMaxDist (theMaxDist)
555 {
556   Init();
557 }
558
559 //=======================================================================
560 //function : Init
561 //purpose  : 
562 //=======================================================================
563
564 void ProjLib_CompProjectedCurve::Init() 
565 {
566   myTabInt.Nullify();
567
568   Standard_Real Tol;// Tolerance for ExactBound
569   Standard_Integer i, Nend = 0;
570   Standard_Boolean FromLastU=Standard_False;
571
572   //new part (to discard far solutions)
573   Standard_Real TolC = Precision::Confusion(), TolS = Precision::Confusion();
574   Extrema_ExtCS CExt(myCurve->Curve(),
575     mySurface->Surface(),
576     TolC,
577     TolS);
578   if (CExt.IsDone() && CExt.NbExt()) 
579   {
580     // Search for the minimum solution
581     Nend = CExt.NbExt();
582     if(myMaxDist > 0 &&
583        // Avoid usage of extrema result that can be wrong for extrusion
584        mySurface->GetType() != GeomAbs_SurfaceOfExtrusion)
585     {
586       Standard_Real min_val2;
587       min_val2 = CExt.SquareDistance(1);
588       for(i = 2; i <= Nend; i++)
589         if (CExt.SquareDistance(i) < min_val2) min_val2 = CExt.SquareDistance(i);
590       if (min_val2 > myMaxDist * myMaxDist)
591         return;
592     }
593   }
594   // end of new part
595
596   Standard_Real FirstU, LastU, Step, SearchStep, WalkStep, t;
597
598   FirstU = myCurve->FirstParameter();
599   LastU  = myCurve->LastParameter();
600   const Standard_Real GlobalMinStep = 1.e-4;
601   //<GlobalMinStep> is sufficiently small to provide solving from initial point
602   //and, on the other hand, it is sufficiently large to avoid too close solutions.
603   const Standard_Real MinStep = 0.01*(LastU - FirstU), 
604     MaxStep = 0.1*(LastU - FirstU);
605   SearchStep = 10*MinStep;
606   Step = SearchStep;
607
608   //Initialization of aPrjPS
609   Standard_Real Uinf = mySurface->FirstUParameter();
610   Standard_Real Usup = mySurface->LastUParameter();
611   Standard_Real Vinf = mySurface->FirstVParameter();
612   Standard_Real Vsup = mySurface->LastVParameter();
613
614   ProjLib_PrjResolve aPrjPS(myCurve->Curve(), mySurface->Surface(), 1);
615
616   t = FirstU;
617   Standard_Boolean new_part; 
618   Standard_Real prevDeb=0.;
619   Standard_Boolean SameDeb=Standard_False;
620
621
622   gp_Pnt Triple, prevTriple;
623
624   //Basic loop  
625   while(t <= LastU) 
626   {
627     //Search for the begining a new continuous part
628     //To avoid infinite computation in some difficult cases
629     new_part = Standard_False;
630     if(t > FirstU && Abs(t-prevDeb) <= Precision::PConfusion()) SameDeb=Standard_True;
631     while(t <= LastU && !new_part && !FromLastU && !SameDeb)
632     {
633       prevDeb=t;
634       if (t == LastU) FromLastU=Standard_True;
635       Standard_Boolean initpoint=Standard_False;
636       Standard_Real U = 0., V = 0.;
637       gp_Pnt CPoint;
638       Standard_Real ParT,ParU,ParV; 
639
640       // Search an initpoint in the list of Extrema Curve-Surface
641       if(Nend != 0 && !CExt.IsParallel()) 
642       {
643         for (i=1;i<=Nend;i++)
644         {
645           Extrema_POnCurv P1;
646           Extrema_POnSurf P2;
647           CExt.Points(i,P1,P2);
648           ParT=P1.Parameter();
649           P2.Parameter(ParU, ParV);
650
651           aPrjPS.Perform(ParT, ParU, ParV, gp_Pnt2d(myTolU, myTolV), 
652             gp_Pnt2d(mySurface->FirstUParameter(),mySurface->FirstVParameter()), 
653             gp_Pnt2d(mySurface->LastUParameter(), mySurface->LastVParameter()), 
654             FuncTol, Standard_True);           
655           if ( aPrjPS.IsDone() && P1.Parameter() > Max(FirstU,t-Step+Precision::PConfusion()) 
656             && P1.Parameter() <= t) 
657           {
658             t=ParT;
659             U=ParU;
660             V=ParV;
661             CPoint=P1.Value();
662             initpoint = Standard_True;
663             break;
664           }
665         }
666       }
667       if (!initpoint) 
668       {        
669         myCurve->D0(t,CPoint);
670 #ifdef OCCT_DEBUG_CHRONO
671         InitChron(chr_init_point);
672 #endif
673         initpoint=InitialPoint(CPoint, t,myCurve,mySurface, myTolU, myTolV, U, V);
674 #ifdef OCCT_DEBUG_CHRONO
675         ResultChron(chr_init_point,t_init_point);
676         init_point_count++;
677 #endif
678       }
679       if(initpoint) 
680       {
681         // When U or V lie on surface joint in some cases we cannot use them 
682         // as initial point for aPrjPS, so we switch them
683         gp_Vec2d D;
684
685         if ((mySurface->IsUPeriodic() &&
686             Abs(Usup - Uinf - mySurface->UPeriod()) < Precision::Confusion()) ||
687             (mySurface->IsVPeriodic() && 
688             Abs(Vsup - Vinf - mySurface->VPeriod()) < Precision::Confusion()))
689         {
690           if((Abs(U - Uinf) < mySurface->UResolution(Precision::PConfusion())) &&
691             mySurface->IsUPeriodic())
692           { 
693             d1(t, U, V, D, myCurve, mySurface);
694             if (D.X() < 0 ) U = Usup;
695           }
696           else if((Abs(U - Usup) < mySurface->UResolution(Precision::PConfusion())) &&
697             mySurface->IsUPeriodic())
698           {
699             d1(t, U, V, D, myCurve, mySurface);
700             if (D.X() > 0) U = Uinf;
701           }
702
703           if((Abs(V - Vinf) < mySurface->VResolution(Precision::PConfusion())) && 
704             mySurface->IsVPeriodic()) 
705           {
706             d1(t, U, V, D, myCurve, mySurface);
707             if (D.Y() < 0) V = Vsup;
708           }
709           else if((Abs(V - Vsup) <= mySurface->VResolution(Precision::PConfusion())) &&
710             mySurface->IsVPeriodic())
711           {
712             d1(t, U, V, D, myCurve, mySurface);
713             if (D.Y() > 0) V = Vinf;
714           }
715         }
716
717
718         if (myMaxDist > 0) 
719         {
720           // Here we are going to stop if the distance between projection and 
721           // corresponding curve point is greater than myMaxDist
722           gp_Pnt POnS;
723           Standard_Real d;
724           mySurface->D0(U, V, POnS);
725           d = CPoint.Distance(POnS);
726           if (d > myMaxDist) 
727           {
728             mySequence->Clear();
729             myNbCurves = 0;
730             return;
731           }
732         }
733         Triple = gp_Pnt(t, U, V);
734         if (t != FirstU) 
735         {
736           //Search for exact boundary point
737           Tol = Min(myTolU, myTolV);
738           gp_Vec2d aD;
739           d1(Triple.X(), Triple.Y(), Triple.Z(), aD, myCurve, mySurface);
740           Tol /= Max(Abs(aD.X()), Abs(aD.Y()));
741
742           if(!ExactBound(Triple, t - Step, Tol, 
743             myTolU, myTolV, myCurve, mySurface)) 
744           {
745 #ifdef OCCT_DEBUG
746             cout<<"There is a problem with ExactBound computation"<<endl;
747 #endif
748             DichExactBound(Triple, t - Step, Tol, myTolU, myTolV, 
749               myCurve, mySurface);
750           }
751         }
752         new_part = Standard_True;
753       }
754       else 
755       {
756         if(t == LastU) break;
757         t += Step;
758         if(t>LastU) 
759         { 
760           Step =Step+LastU-t;
761           t=LastU;
762         }  
763       }
764     }
765     if (!new_part) break;
766
767
768     //We have found a new continuous part
769     Handle(TColgp_HSequenceOfPnt) hSeq = new TColgp_HSequenceOfPnt();    
770     mySequence->Append(hSeq);
771     myNbCurves++;
772     mySequence->Value(myNbCurves)->Append(Triple);
773     prevTriple = Triple;
774
775     if (Triple.X() == LastU) break;//return;
776
777     //Computation of WalkStep
778     gp_Vec D1, D2;
779     Standard_Real MagnD1, MagnD2;
780     d2CurvOnSurf(Triple.X(), Triple.Y(), Triple.Z(), D1, D2, myCurve, mySurface);
781     MagnD1 = D1.Magnitude();
782     MagnD2 = D2.Magnitude();
783     if(MagnD2 < Precision::Confusion()) WalkStep = MaxStep;
784     else WalkStep = Min(MaxStep, Max(MinStep, 0.1*MagnD1/MagnD2));
785
786     Step = WalkStep;
787
788     t = Triple.X() + Step;
789     if (t > LastU) t = LastU;
790     Standard_Real prevStep = Step;
791     Standard_Real U0, V0;
792     gp_Pnt2d aLowBorder(mySurface->FirstUParameter(),mySurface->FirstVParameter());
793     gp_Pnt2d aUppBorder(mySurface->LastUParameter(), mySurface->LastVParameter());
794     gp_Pnt2d aTol(myTolU, myTolV);
795     //Here we are trying to prolong continuous part
796     while (t <= LastU && new_part) 
797     {
798
799       U0 = Triple.Y() + (Step / prevStep) * (Triple.Y() - prevTriple.Y());
800       V0 = Triple.Z() + (Step / prevStep) * (Triple.Z() - prevTriple.Z());
801       // adjust U0 to be in [mySurface->FirstUParameter(),mySurface->LastUParameter()]
802       U0 = Min(Max(U0, aLowBorder.X()), aUppBorder.X()); 
803       // adjust V0 to be in [mySurface->FirstVParameter(),mySurface->LastVParameter()]
804       V0 = Min(Max(V0, aLowBorder.Y()), aUppBorder.Y()); 
805
806
807       aPrjPS.Perform(t, U0, V0, aTol,
808                      aLowBorder, aUppBorder, FuncTol, Standard_True);
809       if(!aPrjPS.IsDone()) 
810       {
811         if (Step <= GlobalMinStep)
812         {
813           //Search for exact boundary point
814           Tol = Min(myTolU, myTolV);
815           gp_Vec2d D;
816           d1(Triple.X(), Triple.Y(), Triple.Z(), D, myCurve, mySurface);
817           Tol /= Max(Abs(D.X()), Abs(D.Y()));
818
819           if(!ExactBound(Triple, t, Tol, myTolU, myTolV, 
820             myCurve, mySurface)) 
821           {
822 #ifdef OCCT_DEBUG
823             cout<<"There is a problem with ExactBound computation"<<endl;
824 #endif
825             DichExactBound(Triple, t, Tol, myTolU, myTolV, 
826               myCurve, mySurface);
827           }
828
829           if((Triple.X() - mySequence->Value(myNbCurves)->Value(mySequence->Value(myNbCurves)->Length()).X()) > 1.e-10)
830             mySequence->Value(myNbCurves)->Append(Triple);
831           if((LastU - Triple.X()) < Tol) {t = LastU + 1; break;}//return;
832
833           Step = SearchStep;
834           t = Triple.X() + Step;
835           if (t > (LastU-MinStep/2) ) 
836           { 
837             Step =Step+LastU-t;
838             t = LastU;
839           }
840           new_part = Standard_False;
841         }
842         else 
843         {
844           // decrease step
845           Standard_Real SaveStep = Step;
846           Step /= 2.;
847           t = Triple .X() + Step;
848           if (t > (LastU-MinStep/4) ) 
849           { 
850             Step =Step+LastU-t;
851             if (Abs(Step - SaveStep) <= Precision::PConfusion())
852               Step = GlobalMinStep; //to avoid looping
853             t = LastU;
854           }
855         }
856       }
857       // Go further
858       else 
859       {
860         prevTriple = Triple;
861         prevStep = Step;
862         Triple = gp_Pnt(t, aPrjPS.Solution().X(), aPrjPS.Solution().Y());
863
864         if (mySurface->GetType() == GeomAbs_SurfaceOfRevolution &&
865            (Abs (Triple.Z() - mySurface->FirstVParameter()) < Precision::Confusion() ||
866             Abs (Triple.Z() - mySurface->LastVParameter() ) < Precision::Confusion() ))
867         {
868           // Go out from possible attraktor.
869
870           Standard_Real U,V;
871           InitialPoint(myCurve->Value(t), t, myCurve, mySurface, myTolU, myTolV, U, V);
872           if (Abs (Abs(U - Triple.Y()) - mySurface->UPeriod()) < Precision::Confusion())
873           {
874             // Handle period jump.
875             U = Triple.Y();
876           }
877           Triple.SetY(U);
878           Triple.SetZ(V);
879         }
880
881         if((Triple.X() - mySequence->Value(myNbCurves)->Value(mySequence->Value(myNbCurves)->Length()).X()) > 1.e-10)
882           mySequence->Value(myNbCurves)->Append(Triple);
883         if (t == LastU) {t = LastU + 1; break;}//return;
884
885         //Computation of WalkStep
886         d2CurvOnSurf(Triple.X(), Triple.Y(), Triple.Z(), D1, D2, myCurve, mySurface);
887         MagnD1 = D1.Magnitude();
888         MagnD2 = D2.Magnitude();
889         if(MagnD2 < Precision::Confusion() ) WalkStep = MaxStep;
890         else WalkStep = Min(MaxStep, Max(MinStep, 0.1*MagnD1/MagnD2));
891
892         Step = WalkStep;
893         t += Step;
894         if (t > (LastU-MinStep/2) ) 
895         {
896           Step =Step+LastU-t;
897           t = LastU;
898         }       
899       }
900     }
901   }
902   // Sequence postproceeding
903   Standard_Integer j;
904
905   // 1. Removing poor parts
906   Standard_Integer NbPart=myNbCurves;
907   Standard_Integer ipart=1;
908   for(i = 1; i <= NbPart; i++) {
909     //    Standard_Integer NbPoints = mySequence->Value(i)->Length();
910     if(mySequence->Value(ipart)->Length() < 2) {
911       mySequence->Remove(ipart);
912       myNbCurves--;
913     }
914     else ipart++;
915   }
916
917   if(myNbCurves == 0) return;
918
919   // 2. Removing common parts of bounds
920   for(i = 1; i < myNbCurves; i++) 
921   {
922     if(mySequence->Value(i)->Value(mySequence->Value(i)->Length()).X() >= 
923       mySequence->Value(i+1)->Value(1).X())
924       mySequence->ChangeValue(i+1)->ChangeValue(1).SetX(mySequence->Value(i)->Value(mySequence->Value(i)->Length()).X() + 1.e-12);
925   }
926
927   // 3. Computation of the maximum distance from each part of curve to surface
928
929   myMaxDistance = new TColStd_HArray1OfReal(1, myNbCurves);
930   myMaxDistance->Init(0);
931   for(i = 1; i <= myNbCurves; i++)
932     for(j = 1; j <= mySequence->Value(i)->Length(); j++) 
933     {
934       gp_Pnt POnC, POnS, aTriple;
935       Standard_Real Distance;
936       aTriple = mySequence->Value(i)->Value(j);
937       myCurve->D0(aTriple.X(), POnC);
938       mySurface->D0(aTriple.Y(), aTriple.Z(), POnS);
939       Distance = POnC.Distance(POnS);
940       if (myMaxDistance->Value(i) < Distance)
941         myMaxDistance->ChangeValue(i) = Distance;
942     } 
943
944
945     // 4. Check the projection to be a single point
946
947     gp_Pnt2d Pmoy, Pcurr, P;
948     Standard_Real AveU, AveV;
949     mySnglPnts = new TColStd_HArray1OfBoolean(1, myNbCurves);
950     for(i = 1; i <= myNbCurves; i++) mySnglPnts->SetValue(i, Standard_True);
951
952     for(i = 1; i <= myNbCurves; i++)
953     {    
954       //compute an average U and V
955
956       for(j = 1, AveU = 0., AveV = 0.; j <= mySequence->Value(i)->Length(); j++)
957       {
958         AveU += mySequence->Value(i)->Value(j).Y();
959         AveV += mySequence->Value(i)->Value(j).Z();
960       }
961       AveU /= mySequence->Value(i)->Length();
962       AveV /= mySequence->Value(i)->Length();
963
964       Pmoy.SetCoord(AveU,AveV);
965       for(j = 1; j <= mySequence->Value(i)->Length(); j++)
966       {
967         Pcurr = 
968           gp_Pnt2d(mySequence->Value(i)->Value(j).Y(), mySequence->Value(i)->Value(j).Z());
969         if (Pcurr.Distance(Pmoy) > ((myTolU < myTolV) ? myTolV : myTolU))
970         {
971           mySnglPnts->SetValue(i, Standard_False);
972           break;
973         }
974       }
975     }
976
977     // 5. Check the projection to be an isoparametric curve of the surface
978
979     myUIso = new TColStd_HArray1OfBoolean(1, myNbCurves);
980     for(i = 1; i <= myNbCurves; i++) myUIso->SetValue(i, Standard_True);
981
982     myVIso = new TColStd_HArray1OfBoolean(1, myNbCurves);
983     for(i = 1; i <= myNbCurves; i++) myVIso->SetValue(i, Standard_True);
984
985     for(i = 1; i <= myNbCurves; i++) {
986       if (IsSinglePnt(i, P)|| mySequence->Value(i)->Length() <=2) {
987         myUIso->SetValue(i, Standard_False);
988         myVIso->SetValue(i, Standard_False);
989         continue;
990       }
991
992       // new test for isoparametrics
993
994       if ( mySequence->Value(i)->Length() > 2) {
995         //compute an average U and V
996
997         for(j = 1, AveU = 0., AveV = 0.; j <= mySequence->Value(i)->Length(); j++) {
998           AveU += mySequence->Value(i)->Value(j).Y();
999           AveV += mySequence->Value(i)->Value(j).Z();
1000         }
1001         AveU /= mySequence->Value(i)->Length();
1002         AveV /= mySequence->Value(i)->Length();
1003
1004         // is i-part U-isoparametric ?
1005         for(j = 1; j <= mySequence->Value(i)->Length(); j++) 
1006         {
1007           if(Abs(mySequence->Value(i)->Value(j).Y() - AveU) > myTolU) 
1008           {
1009             myUIso->SetValue(i, Standard_False);
1010             break;
1011           }
1012         }
1013
1014         // is i-part V-isoparametric ?
1015         for(j = 1; j <= mySequence->Value(i)->Length(); j++) 
1016         {
1017           if(Abs(mySequence->Value(i)->Value(j).Z() - AveV) > myTolV) 
1018           {
1019             myVIso->SetValue(i, Standard_False);
1020             break;
1021           }
1022         }
1023         //
1024       }
1025     }
1026 }
1027 //=======================================================================
1028 //function : Load
1029 //purpose  : 
1030 //=======================================================================
1031
1032 void ProjLib_CompProjectedCurve::Load(const Handle(Adaptor3d_HSurface)& S) 
1033 {
1034   mySurface = S;
1035 }
1036
1037 //=======================================================================
1038 //function : Load
1039 //purpose  : 
1040 //=======================================================================
1041
1042 void ProjLib_CompProjectedCurve::Load(const Handle(Adaptor3d_HCurve)& C) 
1043 {
1044   myCurve = C;
1045 }
1046
1047 //=======================================================================
1048 //function : GetSurface
1049 //purpose  : 
1050 //=======================================================================
1051
1052 const Handle(Adaptor3d_HSurface)& ProjLib_CompProjectedCurve::GetSurface() const
1053 {
1054   return mySurface;
1055 }
1056
1057
1058 //=======================================================================
1059 //function : GetCurve
1060 //purpose  : 
1061 //=======================================================================
1062
1063 const Handle(Adaptor3d_HCurve)& ProjLib_CompProjectedCurve::GetCurve() const
1064 {
1065   return myCurve;
1066 }
1067
1068 //=======================================================================
1069 //function : GetTolerance
1070 //purpose  : 
1071 //=======================================================================
1072
1073 void ProjLib_CompProjectedCurve::GetTolerance(Standard_Real& TolU,
1074   Standard_Real& TolV) const
1075 {
1076   TolU = myTolU;
1077   TolV = myTolV;
1078 }
1079
1080 //=======================================================================
1081 //function : NbCurves
1082 //purpose  : 
1083 //=======================================================================
1084
1085 Standard_Integer ProjLib_CompProjectedCurve::NbCurves() const
1086 {
1087   return myNbCurves;
1088 }
1089 //=======================================================================
1090 //function : Bounds
1091 //purpose  : 
1092 //=======================================================================
1093
1094 void ProjLib_CompProjectedCurve::Bounds(const Standard_Integer Index,
1095   Standard_Real& Udeb,
1096   Standard_Real& Ufin) const
1097 {
1098   if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1099   Udeb = mySequence->Value(Index)->Value(1).X();
1100   Ufin = mySequence->Value(Index)->Value(mySequence->Value(Index)->Length()).X();
1101 }
1102 //=======================================================================
1103 //function : IsSinglePnt
1104 //purpose  : 
1105 //=======================================================================
1106
1107 Standard_Boolean ProjLib_CompProjectedCurve::IsSinglePnt(const Standard_Integer Index, gp_Pnt2d& P) const
1108 {
1109   if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1110   P = gp_Pnt2d(mySequence->Value(Index)->Value(1).Y(), mySequence->Value(Index)->Value(1).Z());
1111   return mySnglPnts->Value(Index);
1112 }
1113
1114 //=======================================================================
1115 //function : IsUIso
1116 //purpose  : 
1117 //=======================================================================
1118
1119 Standard_Boolean ProjLib_CompProjectedCurve::IsUIso(const Standard_Integer Index, Standard_Real& U) const
1120 {
1121   if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1122   U = mySequence->Value(Index)->Value(1).Y();
1123   return myUIso->Value(Index);
1124 }
1125 //=======================================================================
1126 //function : IsVIso
1127 //purpose  : 
1128 //=======================================================================
1129
1130 Standard_Boolean ProjLib_CompProjectedCurve::IsVIso(const Standard_Integer Index, Standard_Real& V) const
1131 {
1132   if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1133   V = mySequence->Value(Index)->Value(1).Z();
1134   return myVIso->Value(Index);
1135 }
1136 //=======================================================================
1137 //function : Value
1138 //purpose  : 
1139 //=======================================================================
1140
1141 gp_Pnt2d ProjLib_CompProjectedCurve::Value(const Standard_Real t) const
1142 {
1143   gp_Pnt2d P;
1144   D0(t, P);
1145   return P;
1146 }
1147 //=======================================================================
1148 //function : D0
1149 //purpose  : 
1150 //=======================================================================
1151
1152 void ProjLib_CompProjectedCurve::D0(const Standard_Real U,gp_Pnt2d& P) const
1153 {
1154   Standard_Integer i, j;
1155   Standard_Real Udeb, Ufin;
1156   Standard_Boolean found = Standard_False;
1157
1158   for(i = 1; i <= myNbCurves; i++) 
1159   {
1160     Bounds(i, Udeb, Ufin);
1161     if (U >= Udeb && U <= Ufin) 
1162     {
1163       found = Standard_True;
1164       break;
1165     }
1166   }
1167   if (!found) Standard_DomainError::Raise("ProjLib_CompProjectedCurve::D0");
1168
1169   Standard_Real U0, V0;
1170
1171   Standard_Integer End = mySequence->Value(i)->Length();
1172   for(j = 1; j < End; j++)
1173     if ((U >= mySequence->Value(i)->Value(j).X()) && (U <= mySequence->Value(i)->Value(j + 1).X())) break;
1174
1175   //  U0 = mySequence->Value(i)->Value(j).Y();
1176   //  V0 = mySequence->Value(i)->Value(j).Z();
1177
1178   //  Cubic Interpolation
1179   if(mySequence->Value(i)->Length() < 4 || 
1180     (Abs(U-mySequence->Value(i)->Value(j).X()) <= Precision::PConfusion()) ) 
1181   {
1182     U0 = mySequence->Value(i)->Value(j).Y();
1183     V0 = mySequence->Value(i)->Value(j).Z();
1184   }
1185   else if (Abs(U-mySequence->Value(i)->Value(j+1).X()) 
1186     <= Precision::PConfusion())
1187   {
1188     U0 = mySequence->Value(i)->Value(j+1).Y();
1189     V0 = mySequence->Value(i)->Value(j+1).Z();
1190   }
1191   else 
1192   {
1193     if (j == 1) j = 2;
1194     if (j > mySequence->Value(i)->Length() - 2) 
1195       j = mySequence->Value(i)->Length() - 2;
1196
1197     gp_Vec2d I1, I2, I3, I21, I22, I31, Y1, Y2, Y3, Y4, Res;
1198     Standard_Real X1, X2, X3, X4;
1199
1200     X1 = mySequence->Value(i)->Value(j - 1).X();
1201     X2 = mySequence->Value(i)->Value(j).X();
1202     X3 = mySequence->Value(i)->Value(j + 1).X();
1203     X4 = mySequence->Value(i)->Value(j + 2).X();
1204
1205     Y1 = gp_Vec2d(mySequence->Value(i)->Value(j - 1).Y(), 
1206       mySequence->Value(i)->Value(j - 1).Z());
1207     Y2 = gp_Vec2d(mySequence->Value(i)->Value(j).Y(), 
1208       mySequence->Value(i)->Value(j).Z());
1209     Y3 = gp_Vec2d(mySequence->Value(i)->Value(j + 1).Y(), 
1210       mySequence->Value(i)->Value(j + 1).Z());
1211     Y4 = gp_Vec2d(mySequence->Value(i)->Value(j + 2).Y(), 
1212       mySequence->Value(i)->Value(j + 2).Z());
1213
1214     I1 = (Y1 - Y2)/(X1 - X2);
1215     I2 = (Y2 - Y3)/(X2 - X3);
1216     I3 = (Y3 - Y4)/(X3 - X4);
1217
1218     I21 = (I1 - I2)/(X1 - X3);
1219     I22 = (I2 - I3)/(X2 - X4);
1220
1221     I31 = (I21 - I22)/(X1 - X4);
1222
1223     Res = Y1 + (U - X1)*(I1 + (U - X2)*(I21 + (U - X3)*I31));
1224
1225     U0 = Res.X();
1226     V0 = Res.Y();
1227
1228     if(U0 < mySurface->FirstUParameter()) U0 = mySurface->FirstUParameter();
1229     else if(U0 > mySurface->LastUParameter()) U0 = mySurface->LastUParameter();
1230
1231     if(V0 < mySurface->FirstVParameter()) V0 = mySurface->FirstVParameter();
1232     else if(V0 > mySurface->LastVParameter()) V0 = mySurface->LastVParameter();
1233   }
1234   //End of cubic interpolation
1235
1236   ProjLib_PrjResolve aPrjPS(myCurve->Curve(), mySurface->Surface(), 1);
1237   aPrjPS.Perform(U, U0, V0, gp_Pnt2d(myTolU, myTolV), 
1238     gp_Pnt2d(mySurface->FirstUParameter(), mySurface->FirstVParameter()), 
1239     gp_Pnt2d(mySurface->LastUParameter(), mySurface->LastVParameter()));
1240   if (aPrjPS.IsDone())
1241     P = aPrjPS.Solution();
1242   else
1243   {
1244     gp_Pnt thePoint = myCurve->Value(U);
1245     Extrema_ExtPS aExtPS(thePoint, mySurface->Surface(), myTolU, myTolV);
1246     if (aExtPS.IsDone() && aExtPS.NbExt()) 
1247     {
1248       Standard_Integer k, Nend, imin = 1;
1249       // Search for the nearest solution which is also a normal projection
1250       Nend = aExtPS.NbExt();
1251       for(k = 2; k <= Nend; k++)
1252         if (aExtPS.SquareDistance(k) < aExtPS.SquareDistance(imin))
1253           imin = k;
1254       const Extrema_POnSurf& POnS = aExtPS.Point(imin);
1255       Standard_Real ParU,ParV;
1256       POnS.Parameter(ParU, ParV);
1257       P.SetCoord(ParU, ParV);
1258     }
1259     else
1260       P.SetCoord(U0,V0);
1261   }
1262 }
1263 //=======================================================================
1264 //function : D1
1265 //purpose  : 
1266 //=======================================================================
1267
1268 void ProjLib_CompProjectedCurve::D1(const Standard_Real t,
1269   gp_Pnt2d& P,
1270   gp_Vec2d& V) const
1271 {
1272   Standard_Real u, v;
1273   D0(t, P);
1274   u = P.X();
1275   v = P.Y();
1276   d1(t, u, v, V, myCurve, mySurface);
1277 }
1278 //=======================================================================
1279 //function : D2
1280 //purpose  : 
1281 //=======================================================================
1282
1283 void ProjLib_CompProjectedCurve::D2(const Standard_Real t,
1284   gp_Pnt2d& P,
1285   gp_Vec2d& V1,
1286   gp_Vec2d& V2) const
1287 {
1288   Standard_Real u, v;
1289   D0(t, P);
1290   u = P.X();
1291   v = P.Y();
1292   d2(t, u, v, V1, V2, myCurve, mySurface);
1293 }
1294 //=======================================================================
1295 //function : DN
1296 //purpose  : 
1297 //=======================================================================
1298
1299 gp_Vec2d ProjLib_CompProjectedCurve::DN(const Standard_Real t, 
1300   const Standard_Integer N) const 
1301 {
1302   if (N < 1 ) Standard_OutOfRange::Raise("ProjLib_CompProjectedCurve : N must be greater than 0");
1303   else if (N ==1) 
1304   {
1305     gp_Pnt2d P;
1306     gp_Vec2d V;
1307     D1(t,P,V);
1308     return V;
1309   }
1310   else if ( N==2)
1311   {
1312     gp_Pnt2d P;
1313     gp_Vec2d V1,V2;
1314     D2(t,P,V1,V2);
1315     return V2;
1316   }
1317   else if (N > 2 ) 
1318     Standard_NotImplemented::Raise("ProjLib_CompProjectedCurve::DN");
1319   return gp_Vec2d();
1320 }
1321
1322 //=======================================================================
1323 //function : GetSequence
1324 //purpose  : 
1325 //=======================================================================
1326
1327 const Handle(ProjLib_HSequenceOfHSequenceOfPnt)& ProjLib_CompProjectedCurve::GetSequence() const
1328 {
1329   return mySequence;
1330 }
1331 //=======================================================================
1332 //function : FirstParameter
1333 //purpose  : 
1334 //=======================================================================
1335
1336 Standard_Real ProjLib_CompProjectedCurve::FirstParameter() const
1337 {
1338   return myCurve->FirstParameter();
1339 }
1340
1341 //=======================================================================
1342 //function : LastParameter
1343 //purpose  : 
1344 //=======================================================================
1345
1346 Standard_Real ProjLib_CompProjectedCurve::LastParameter() const
1347 {
1348   return myCurve->LastParameter();
1349 }
1350
1351 //=======================================================================
1352 //function : MaxDistance
1353 //purpose  : 
1354 //=======================================================================
1355
1356 Standard_Real ProjLib_CompProjectedCurve::MaxDistance(const Standard_Integer Index) const
1357 {
1358   if(Index < 1 || Index > myNbCurves) Standard_NoSuchObject::Raise();
1359   return myMaxDistance->Value(Index);
1360 }
1361
1362 //=======================================================================
1363 //function : NbIntervals
1364 //purpose  : 
1365 //=======================================================================
1366
1367 Standard_Integer ProjLib_CompProjectedCurve::NbIntervals(const GeomAbs_Shape S) const
1368 {
1369   const_cast<ProjLib_CompProjectedCurve*>(this)->myTabInt.Nullify();
1370   BuildIntervals(S);
1371   return myTabInt->Length() - 1;
1372 }
1373
1374 //=======================================================================
1375 //function : Intervals
1376 //purpose  : 
1377 //=======================================================================
1378
1379 void ProjLib_CompProjectedCurve::Intervals(TColStd_Array1OfReal& T,const GeomAbs_Shape S) const
1380 {
1381   if (myTabInt.IsNull()) BuildIntervals (S);
1382   T = myTabInt->Array1();
1383 }
1384
1385 //=======================================================================
1386 //function : BuildIntervals
1387 //purpose  : 
1388 //=======================================================================
1389
1390 void ProjLib_CompProjectedCurve::BuildIntervals(const GeomAbs_Shape S) const
1391 {
1392   GeomAbs_Shape SforS = GeomAbs_CN;
1393   switch(S) {
1394   case GeomAbs_C0: 
1395     SforS = GeomAbs_C1; 
1396     break;    
1397   case GeomAbs_C1: 
1398     SforS = GeomAbs_C2; 
1399     break;
1400   case GeomAbs_C2: 
1401     SforS = GeomAbs_C3; 
1402     break;
1403   case GeomAbs_C3:
1404     SforS = GeomAbs_CN; 
1405     break;
1406   case GeomAbs_CN: 
1407     SforS = GeomAbs_CN; 
1408     break;
1409   default: 
1410     Standard_OutOfRange::Raise();
1411   }
1412   Standard_Integer i, j, k;
1413   Standard_Integer NbIntCur = myCurve->NbIntervals(S);
1414   Standard_Integer NbIntSurU = mySurface->NbUIntervals(SforS);
1415   Standard_Integer NbIntSurV = mySurface->NbVIntervals(SforS);
1416
1417   TColStd_Array1OfReal CutPntsT(1, NbIntCur+1);
1418   TColStd_Array1OfReal CutPntsU(1, NbIntSurU+1);
1419   TColStd_Array1OfReal CutPntsV(1, NbIntSurV+1);
1420
1421   myCurve->Intervals(CutPntsT, S);
1422   mySurface->UIntervals(CutPntsU, SforS);
1423   mySurface->VIntervals(CutPntsV, SforS);
1424
1425   Standard_Real Tl, Tr, Ul, Ur, Vl, Vr, Tol;
1426
1427   Handle(TColStd_HArray1OfReal) BArr = NULL, 
1428     CArr = NULL, 
1429     UArr = NULL, 
1430     VArr = NULL;
1431
1432   // proccessing projection bounds
1433   BArr = new TColStd_HArray1OfReal(1, 2*myNbCurves);
1434   for(i = 1; i <= myNbCurves; i++)
1435     Bounds(i, BArr->ChangeValue(2*i - 1), BArr->ChangeValue(2*i));
1436
1437   // proccessing curve discontinuities
1438   if(NbIntCur > 1) {
1439     CArr = new TColStd_HArray1OfReal(1, NbIntCur - 1);
1440     for(i = 1; i <= CArr->Length(); i++)
1441       CArr->ChangeValue(i) = CutPntsT(i + 1);
1442   }
1443
1444   // proccessing U-surface discontinuities  
1445   TColStd_SequenceOfReal TUdisc;
1446
1447   for(k = 2; k <= NbIntSurU; k++) {
1448     //    cout<<"CutPntsU("<<k<<") = "<<CutPntsU(k)<<endl;
1449     for(i = 1; i <= myNbCurves; i++)
1450       for(j = 1; j < mySequence->Value(i)->Length(); j++) {
1451         Ul = mySequence->Value(i)->Value(j).Y();
1452         Ur = mySequence->Value(i)->Value(j + 1).Y();
1453
1454         if(Abs(Ul - CutPntsU(k)) <= myTolU) 
1455           TUdisc.Append(mySequence->Value(i)->Value(j).X());
1456         else if(Abs(Ur - CutPntsU(k)) <= myTolU) 
1457           TUdisc.Append(mySequence->Value(i)->Value(j + 1).X());
1458         else if((Ul < CutPntsU(k) && CutPntsU(k) < Ur) ||
1459           (Ur < CutPntsU(k) && CutPntsU(k) < Ul)) 
1460         {
1461           Standard_Real V;
1462           V = (mySequence->Value(i)->Value(j).Z() 
1463             + mySequence->Value(i)->Value(j +1).Z())/2;
1464           ProjLib_PrjResolve Solver(myCurve->Curve(), mySurface->Surface(), 2);
1465
1466           gp_Vec2d D;
1467           gp_Pnt Triple;
1468           Triple = mySequence->Value(i)->Value(j);
1469           d1(Triple.X(), Triple.Y(), Triple.Z(), D, myCurve, mySurface);
1470           if (Abs(D.X()) < Precision::Confusion()) 
1471             Tol = myTolU;
1472           else 
1473             Tol = Min(myTolU, myTolU / Abs(D.X()));
1474
1475           Tl = mySequence->Value(i)->Value(j).X();
1476           Tr = mySequence->Value(i)->Value(j + 1).X();
1477
1478           Solver.Perform((Tl + Tr)/2, CutPntsU(k), V, 
1479             gp_Pnt2d(Tol, myTolV), 
1480             gp_Pnt2d(Tl, mySurface->FirstVParameter()), 
1481             gp_Pnt2d(Tr, mySurface->LastVParameter()));
1482           //
1483           if(Solver.IsDone()) 
1484           {
1485             TUdisc.Append(Solver.Solution().X());
1486           }
1487         }
1488       }
1489   }
1490   for(i = 2; i <= TUdisc.Length(); i++)
1491     if(TUdisc(i) - TUdisc(i-1) < Precision::PConfusion())
1492       TUdisc.Remove(i--);
1493
1494   if(TUdisc.Length()) 
1495   {
1496     UArr = new TColStd_HArray1OfReal(1, TUdisc.Length());
1497     for(i = 1; i <= UArr->Length(); i++)
1498       UArr->ChangeValue(i) = TUdisc(i);
1499   }
1500   // proccessing V-surface discontinuities
1501   TColStd_SequenceOfReal TVdisc;
1502
1503   for(k = 2; k <= NbIntSurV; k++)
1504     for(i = 1; i <= myNbCurves; i++) 
1505     {
1506       //      cout<<"CutPntsV("<<k<<") = "<<CutPntsV(k)<<endl;
1507       for(j = 1; j < mySequence->Value(i)->Length(); j++) {
1508
1509         Vl = mySequence->Value(i)->Value(j).Z();
1510         Vr = mySequence->Value(i)->Value(j + 1).Z();
1511
1512         if(Abs(Vl - CutPntsV(k)) <= myTolV) 
1513           TVdisc.Append(mySequence->Value(i)->Value(j).X());
1514         else if (Abs(Vr - CutPntsV(k)) <= myTolV) 
1515           TVdisc.Append(mySequence->Value(i)->Value(j + 1).X());
1516         else if((Vl < CutPntsV(k) && CutPntsV(k) < Vr) ||
1517           (Vr < CutPntsV(k) && CutPntsV(k) < Vl)) 
1518         {
1519           Standard_Real U;
1520           U = (mySequence->Value(i)->Value(j).Y() 
1521             + mySequence->Value(i)->Value(j +1).Y())/2;
1522           ProjLib_PrjResolve Solver(myCurve->Curve(), mySurface->Surface(), 3);
1523
1524           gp_Vec2d D;
1525           gp_Pnt Triple;
1526           Triple = mySequence->Value(i)->Value(j);
1527           d1(Triple.X(), Triple.Y(), Triple.Z(), D, myCurve, mySurface);
1528           if (Abs(D.Y()) < Precision::Confusion()) 
1529             Tol = myTolV;
1530           else 
1531             Tol = Min(myTolV, myTolV / Abs(D.Y()));
1532
1533           Tl = mySequence->Value(i)->Value(j).X();
1534           Tr = mySequence->Value(i)->Value(j + 1).X();
1535
1536           Solver.Perform((Tl + Tr)/2, U, CutPntsV(k), 
1537             gp_Pnt2d(Tol, myTolV), 
1538             gp_Pnt2d(Tl, mySurface->FirstUParameter()), 
1539             gp_Pnt2d(Tr, mySurface->LastUParameter()));
1540           //
1541           if(Solver.IsDone()) 
1542           {
1543             TVdisc.Append(Solver.Solution().X());
1544           }
1545         }
1546       }
1547     }
1548     for(i = 2; i <= TVdisc.Length(); i++)
1549       if(TVdisc(i) - TVdisc(i-1) < Precision::PConfusion())
1550         TVdisc.Remove(i--);
1551
1552     if(TVdisc.Length()) 
1553     {
1554       VArr = new TColStd_HArray1OfReal(1, TVdisc.Length());
1555       for(i = 1; i <= VArr->Length(); i++)
1556         VArr->ChangeValue(i) = TVdisc(i);
1557     }
1558
1559     // fusion
1560     TColStd_SequenceOfReal Fusion;
1561     if(!CArr.IsNull()) 
1562     {
1563       GeomLib::FuseIntervals(BArr->ChangeArray1(), 
1564         CArr->ChangeArray1(), 
1565         Fusion, Precision::PConfusion());
1566       BArr = new TColStd_HArray1OfReal(1, Fusion.Length());
1567       for(i = 1; i <= BArr->Length(); i++)
1568         BArr->ChangeValue(i) = Fusion(i);
1569       Fusion.Clear();
1570     }
1571
1572     if(!UArr.IsNull()) 
1573     {
1574       GeomLib::FuseIntervals(BArr->ChangeArray1(), 
1575         UArr->ChangeArray1(), 
1576         Fusion, Precision::PConfusion());
1577       BArr = new TColStd_HArray1OfReal(1, Fusion.Length());
1578       for(i = 1; i <= BArr->Length(); i++)
1579         BArr->ChangeValue(i) = Fusion(i);
1580       Fusion.Clear();
1581     }
1582
1583     if(!VArr.IsNull()) 
1584     {
1585       GeomLib::FuseIntervals(BArr->ChangeArray1(), 
1586         VArr->ChangeArray1(), 
1587         Fusion, Precision::PConfusion());
1588       BArr = new TColStd_HArray1OfReal(1, Fusion.Length());
1589       for(i = 1; i <= BArr->Length(); i++)
1590         BArr->ChangeValue(i) = Fusion(i);
1591     }
1592
1593     const_cast<ProjLib_CompProjectedCurve*>(this)->myTabInt = new TColStd_HArray1OfReal(1, BArr->Length());
1594     for(i = 1; i <= BArr->Length(); i++)
1595       myTabInt->ChangeValue(i) = BArr->Value(i);
1596
1597 }
1598
1599 //=======================================================================
1600 //function : Trim
1601 //purpose  : 
1602 //=======================================================================
1603
1604 Handle(Adaptor2d_HCurve2d) ProjLib_CompProjectedCurve::Trim
1605   (const Standard_Real First,
1606   const Standard_Real Last,
1607   const Standard_Real Tol) const 
1608 {
1609   Handle(ProjLib_HCompProjectedCurve) HCS = 
1610     new ProjLib_HCompProjectedCurve(*this);
1611   HCS->ChangeCurve2d().Load(mySurface);
1612   HCS->ChangeCurve2d().Load(myCurve->Trim(First,Last,Tol));
1613   return HCS;
1614 }
1615
1616 //=======================================================================
1617 //function : GetType
1618 //purpose  : 
1619 //=======================================================================
1620
1621 GeomAbs_CurveType ProjLib_CompProjectedCurve::GetType() const 
1622 {
1623   return GeomAbs_OtherCurve;
1624 }