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