e5b36a4bd9e9b66df4d7ae27881d9c30b61b3b00
[occt.git] / src / IntWalk / IntWalk_PWalking.cxx
1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
15
16 #include <Adaptor3d_HSurface.hxx>
17 #include <Adaptor3d_HSurfaceTool.hxx>
18 #include <Extrema_GenLocateExtPS.hxx>
19 #include <Geom_Surface.hxx>
20 #include <gp_Dir.hxx>
21 #include <gp_Pnt.hxx>
22 #include <gp_Pnt2d.hxx>
23 #include <IntImp_ComputeTangence.hxx>
24 #include <IntSurf_LineOn2S.hxx>
25 #include <IntSurf_PntOn2S.hxx>
26 #include <IntWalk_PWalking.hxx>
27 #include <IntWalk_StatusDeflection.hxx>
28 #include <math_FunctionSetRoot.hxx>
29 #include <Precision.hxx>
30 #include <Standard_Failure.hxx>
31 #include <Standard_OutOfRange.hxx>
32 #include <StdFail_NotDone.hxx>
33 #include <TColgp_Array1OfPnt.hxx>
34 #include <TColStd_Array1OfReal.hxx>
35
36 //==================================================================================
37 // function : ComputePasInit
38 // purpose  : estimate of max step : To avoid abrupt changes during change of isos 
39 //==================================================================================
40 void IntWalk_PWalking::ComputePasInit(const Standard_Real theDeltaU1,
41                                       const Standard_Real theDeltaV1,
42                                       const Standard_Real theDeltaU2,
43                                       const Standard_Real theDeltaV2)
44 {
45   const Standard_Real aRangePart = 0.01;
46   const Standard_Real Increment = 2.0*pasMax;
47   const Handle(Adaptor3d_HSurface)& 
48           Caro1 = myIntersectionOn2S.Function().AuxillarSurface1();
49   const Handle(Adaptor3d_HSurface)& 
50           Caro2 = myIntersectionOn2S.Function().AuxillarSurface2();
51
52   const Standard_Real aDeltaU1=Abs(UM1-Um1);
53   const Standard_Real aDeltaV1=Abs(VM1-Vm1);
54   const Standard_Real aDeltaU2=Abs(UM2-Um2);
55   const Standard_Real aDeltaV2=Abs(VM2-Vm2);
56
57   //-- limit the reduction of uv box estimate to 0.01 natural box
58   //--  theDeltaU1 : On box of Inter
59   //-- aDeltaU1 : On parametric space
60   if(!Precision::IsInfinite(aDeltaU1))
61     pasuv[0]=Max(Increment*Max(theDeltaU1, aRangePart*aDeltaU1), pasuv[0]);
62   else
63     pasuv[0]=Max(Increment*theDeltaU1, pasuv[0]);
64
65   if(!Precision::IsInfinite(aDeltaV1))
66     pasuv[1]=Max(Increment*Max(theDeltaV1, aRangePart*aDeltaV1), pasuv[1]);
67   else
68     pasuv[1]=Max(Increment*theDeltaV1, pasuv[1]);
69
70   if(!Precision::IsInfinite(aDeltaU2))
71     pasuv[2]=Max(Increment*Max(theDeltaU2, aRangePart*aDeltaU2), pasuv[2]);
72   else
73     pasuv[2]=Max(Increment*theDeltaU2, pasuv[2]);
74
75   if(!Precision::IsInfinite(aDeltaV2))
76     pasuv[3]=Max(Increment*Max(theDeltaV2, aRangePart*aDeltaV2), pasuv[3]);
77   else
78     pasuv[3]=Max(Increment*theDeltaV2, pasuv[3]);
79
80   const Standard_Real ResoU1tol = Adaptor3d_HSurfaceTool::UResolution(Caro1, tolconf);
81   const Standard_Real ResoV1tol = Adaptor3d_HSurfaceTool::VResolution(Caro1, tolconf);
82   const Standard_Real ResoU2tol = Adaptor3d_HSurfaceTool::UResolution(Caro2, tolconf);
83   const Standard_Real ResoV2tol = Adaptor3d_HSurfaceTool::VResolution(Caro2, tolconf);
84
85   myStepMin[0] = Max(myStepMin[0], 2.0*ResoU1tol);
86   myStepMin[1] = Max(myStepMin[1], 2.0*ResoV1tol);
87   myStepMin[2] = Max(myStepMin[2], 2.0*ResoU2tol);
88   myStepMin[3] = Max(myStepMin[3], 2.0*ResoV2tol);
89
90   for(Standard_Integer i  = 0; i < 4; i++)
91   {
92     pasuv[i]=Max(myStepMin[i], pasuv[i]);
93   }
94 }
95
96 //=======================================================================
97 //function : IsParallel
98 //purpose  : Checks if theLine is parallel of some boundary of given
99 //            surface (it is determined by theCheckSurf1 flag). 
100 //            Parallelism assumes small oscillations (swing is less or 
101 //            equal than theToler).
102 //            Small lines (if first and last parameters in the Surface 
103 //            are almost equal) are classified as parallel (as same as 
104 //            any point can be considered as parallel of any line).
105 //=======================================================================
106 static void IsParallel(const Handle(IntSurf_LineOn2S)& theLine,
107                        const Standard_Boolean theCheckSurf1,
108                        const Standard_Real theToler,
109                        Standard_Boolean& theIsUparallel,
110                        Standard_Boolean& theIsVparallel)
111 {
112   const Standard_Integer aNbPointsMAX = 23;
113
114   theIsUparallel = theIsVparallel = Standard_True;
115
116   const Standard_Integer aNbLinePnts = theLine->NbPoints();
117   Standard_Integer aNbPoints = aNbLinePnts;
118   if(aNbPoints > aNbPointsMAX)
119   {
120     aNbPoints = aNbPointsMAX;
121   }
122   else if(aNbPoints < 3)
123   {
124     //Here we cannot estimate parallelism.
125     //Do all same as for small lines 
126     return;
127   }
128
129   Standard_Real aStep = IntToReal(theLine->NbPoints()) / aNbPoints;
130   Standard_Real aNPoint = 1.0;
131
132   Standard_Real aUmin = RealLast(), aUmax = RealFirst(), aVmin = RealLast(), aVmax = RealFirst();
133   for(Standard_Integer aNum = 1; aNum <= aNbPoints; aNum++, aNPoint += aStep)
134   {
135     // Fix possible "out of parameter" case.
136     if (aNPoint > aNbLinePnts)
137       aNPoint = aNbLinePnts;
138
139     Standard_Real u, v;
140     if(theCheckSurf1)
141       theLine->Value(RealToInt(aNPoint)).ParametersOnS1(u, v);
142     else
143       theLine->Value(RealToInt(aNPoint)).ParametersOnS2(u, v);
144
145     if(u < aUmin)
146       aUmin = u;
147
148     if(u > aUmax)
149       aUmax = u;
150
151     if(v < aVmin)
152       aVmin = v;
153
154     if(v > aVmax)
155       aVmax = v;
156   }
157
158   theIsVparallel = ((aUmax - aUmin) < theToler);
159   theIsUparallel = ((aVmax - aVmin) < theToler);
160 }
161
162 //=======================================================================
163 //function : AdjustToDomain
164 //purpose  : Returns TRUE if theP has been changed (i.e. initial value
165 //            was out of the domain)
166 //=======================================================================
167 static Standard_Boolean AdjustToDomain(const Standard_Integer theNbElem,
168                                        Standard_Real* theParam,
169                                        const Standard_Real* const theLowBorder,
170                                        const Standard_Real* const theUppBorder)
171 {
172   Standard_Boolean aRetVal = Standard_False;
173   for (Standard_Integer i = 0; i < theNbElem; i++)
174   {
175     if ((theParam[i] - theLowBorder[i]) < -Precision::PConfusion())
176     {
177       theParam[i] = theLowBorder[i];
178       aRetVal = Standard_True;
179     }
180
181     if ((theParam[i] - theUppBorder[i]) > Precision::PConfusion())
182     {
183       theParam[i] = theUppBorder[i];
184       aRetVal = Standard_True;
185     }
186   }
187
188   return aRetVal;
189 }
190
191 //==================================================================================
192 // function : IntWalk_PWalking::IntWalk_PWalking
193 // purpose  : 
194 //==================================================================================
195 IntWalk_PWalking::IntWalk_PWalking(const Handle(Adaptor3d_HSurface)& Caro1,
196                                    const Handle(Adaptor3d_HSurface)& Caro2,
197                                    const Standard_Real TolTangency,
198                                    const Standard_Real Epsilon,
199                                    const Standard_Real Deflection,
200                                    const Standard_Real Increment ) 
201                                    :
202
203 done(Standard_True),
204 close(Standard_False),
205 fleche(Deflection),
206 tolconf(Epsilon),
207 myTolTang(TolTangency),
208 sensCheminement(1),
209 myIntersectionOn2S(Caro1,Caro2,TolTangency),
210 STATIC_BLOCAGE_SUR_PAS_TROP_GRAND(0),
211 STATIC_PRECEDENT_INFLEXION(0)
212 {
213   Standard_Real KELARG=20.;
214   //
215   pasMax=Increment*0.2; //-- June 25 99 after problems with precision 
216   Um1 = Adaptor3d_HSurfaceTool::FirstUParameter(Caro1);
217   Vm1 = Adaptor3d_HSurfaceTool::FirstVParameter(Caro1);
218   UM1 = Adaptor3d_HSurfaceTool::LastUParameter(Caro1);
219   VM1 = Adaptor3d_HSurfaceTool::LastVParameter(Caro1);
220
221   Um2 = Adaptor3d_HSurfaceTool::FirstUParameter(Caro2);
222   Vm2 = Adaptor3d_HSurfaceTool::FirstVParameter(Caro2);
223   UM2 = Adaptor3d_HSurfaceTool::LastUParameter(Caro2);
224   VM2 = Adaptor3d_HSurfaceTool::LastVParameter(Caro2);
225
226   ResoU1 = Adaptor3d_HSurfaceTool::UResolution(Caro1,Precision::Confusion());
227   ResoV1 = Adaptor3d_HSurfaceTool::VResolution(Caro1,Precision::Confusion());
228
229   ResoU2 = Adaptor3d_HSurfaceTool::UResolution(Caro2,Precision::Confusion());
230   ResoV2 = Adaptor3d_HSurfaceTool::VResolution(Caro2,Precision::Confusion());
231
232   Standard_Real NEWRESO;
233   Standard_Real MAXVAL;
234   Standard_Real MAXVAL2;
235   //
236   MAXVAL  = Abs(Um1);  MAXVAL2 = Abs(UM1);
237   if(MAXVAL2 > MAXVAL) MAXVAL = MAXVAL2;
238   NEWRESO = ResoU1 * MAXVAL ;
239   if(NEWRESO > ResoU1 &&NEWRESO<10) {    ResoU1 = NEWRESO;  }
240
241
242   MAXVAL  = Abs(Um2);   MAXVAL2 = Abs(UM2);
243   if(MAXVAL2 > MAXVAL) MAXVAL = MAXVAL2;
244   NEWRESO = ResoU2 * MAXVAL ;
245   if(NEWRESO > ResoU2 && NEWRESO<10) {     ResoU2 = NEWRESO;  }
246
247
248   MAXVAL  = Abs(Vm1);  MAXVAL2 = Abs(VM1);
249   if(MAXVAL2 > MAXVAL) MAXVAL = MAXVAL2;
250   NEWRESO = ResoV1 * MAXVAL ;
251   if(NEWRESO > ResoV1 && NEWRESO<10) {     ResoV1 = NEWRESO;  }
252
253
254   MAXVAL  = Abs(Vm2);  MAXVAL2 = Abs(VM2);
255   if(MAXVAL2 > MAXVAL) MAXVAL = MAXVAL2;
256   NEWRESO = ResoV2 * MAXVAL ;
257   if(NEWRESO > ResoV2 && NEWRESO<10) {     ResoV2 = NEWRESO;  }
258
259   pasuv[0]=pasMax*Abs(UM1-Um1);
260   pasuv[1]=pasMax*Abs(VM1-Vm1);
261   pasuv[2]=pasMax*Abs(UM2-Um2);
262   pasuv[3]=pasMax*Abs(VM2-Vm2);
263
264   if(ResoU1>0.0001*pasuv[0]) ResoU1=0.00001*pasuv[0];
265   if(ResoV1>0.0001*pasuv[1]) ResoV1=0.00001*pasuv[1];
266   if(ResoU2>0.0001*pasuv[2]) ResoU2=0.00001*pasuv[2];
267   if(ResoV2>0.0001*pasuv[3]) ResoV2=0.00001*pasuv[3];
268
269
270   if(Adaptor3d_HSurfaceTool::IsUPeriodic(Caro1)==Standard_False) { 
271     //UM1+=KELARG*pasuv[0];  Um1-=KELARG*pasuv[0];
272   }
273   else { 
274     Standard_Real t = UM1-Um1; 
275     if(t<Adaptor3d_HSurfaceTool::UPeriod(Caro1)) { 
276       t=0.5*(Adaptor3d_HSurfaceTool::UPeriod(Caro1)-t);
277       t=(t>KELARG*pasuv[0])? KELARG*pasuv[0] : t;
278       UM1+=t;  Um1-=t;
279     }
280   }
281
282   if(Adaptor3d_HSurfaceTool::IsVPeriodic(Caro1)==Standard_False) { 
283     //VM1+=KELARG*pasuv[1];  Vm1-=KELARG*pasuv[1];
284   }
285   else { 
286     Standard_Real t = VM1-Vm1; 
287     if(t<Adaptor3d_HSurfaceTool::VPeriod(Caro1)) { 
288       t=0.5*(Adaptor3d_HSurfaceTool::VPeriod(Caro1)-t);
289       t=(t>KELARG*pasuv[1])? KELARG*pasuv[1] : t;
290       VM1+=t;  Vm1-=t;
291     }
292   }
293
294   if(Adaptor3d_HSurfaceTool::IsUPeriodic(Caro2)==Standard_False) { 
295     //UM2+=KELARG*pasuv[2];  Um2-=KELARG*pasuv[2];
296   }
297   else { 
298     Standard_Real t = UM2-Um2; 
299     if(t<Adaptor3d_HSurfaceTool::UPeriod(Caro2)) { 
300       t=0.5*(Adaptor3d_HSurfaceTool::UPeriod(Caro2)-t);
301       t=(t>KELARG*pasuv[2])? KELARG*pasuv[2] : t;
302       UM2+=t;  Um2-=t;
303     }
304   }
305
306   if(Adaptor3d_HSurfaceTool::IsVPeriodic(Caro2)==Standard_False) {   
307     //VM2+=KELARG*pasuv[3];  Vm2-=KELARG*pasuv[3];
308   }
309   else { 
310     Standard_Real t = VM2-Vm2; 
311     if(t<Adaptor3d_HSurfaceTool::VPeriod(Caro2)) { 
312       t=0.5*(Adaptor3d_HSurfaceTool::VPeriod(Caro2)-t);
313       t=(t>KELARG*pasuv[3])? KELARG*pasuv[3] : t;
314       VM2+=t;  Vm2-=t;
315     }
316   }
317
318   myStepMin[0] = 100.0*ResoU1;
319   myStepMin[1] = 100.0*ResoV1;
320   myStepMin[2] = 100.0*ResoU2;
321   myStepMin[3] = 100.0*ResoV2;
322
323   //-- ComputePasInit(pasuv,Um1,UM1,Vm1,VM1,Um2,UM2,Vm2,VM2,Caro1,Caro2);
324
325   for (Standard_Integer i = 0; i<=3;i++) {
326     if(pasuv[i]>10) 
327       pasuv[i] = 10; 
328     pasInit[i] = pasSav[i] = pasuv[i]; 
329   }
330
331
332 }
333 //==================================================================================
334 // function : IntWalk_PWalking
335 // purpose  : 
336 //==================================================================================
337 IntWalk_PWalking::IntWalk_PWalking(const Handle(Adaptor3d_HSurface)& Caro1,
338                                    const Handle(Adaptor3d_HSurface)& Caro2,
339                                    const Standard_Real TolTangency,
340                                    const Standard_Real Epsilon,
341                                    const Standard_Real Deflection,
342                                    const Standard_Real Increment, 
343                                    const Standard_Real U1,
344                                    const Standard_Real V1,
345                                    const Standard_Real U2, 
346                                    const Standard_Real V2)
347                                    :
348
349 done(Standard_True),
350 close(Standard_False),
351 fleche(Deflection),
352 tolconf(Epsilon),
353 myTolTang(TolTangency),
354 sensCheminement(1),
355 myIntersectionOn2S(Caro1,Caro2,TolTangency),
356 STATIC_BLOCAGE_SUR_PAS_TROP_GRAND(0),
357 STATIC_PRECEDENT_INFLEXION(0)
358 {
359   Standard_Real KELARG=20.;
360   //
361   pasMax=Increment*0.2; //-- June 25 99 after problems with precision 
362   //
363   Um1 = Adaptor3d_HSurfaceTool::FirstUParameter(Caro1);
364   Vm1 = Adaptor3d_HSurfaceTool::FirstVParameter(Caro1);
365   UM1 = Adaptor3d_HSurfaceTool::LastUParameter(Caro1);
366   VM1 = Adaptor3d_HSurfaceTool::LastVParameter(Caro1);
367
368   Um2 = Adaptor3d_HSurfaceTool::FirstUParameter(Caro2);
369   Vm2 = Adaptor3d_HSurfaceTool::FirstVParameter(Caro2);
370   UM2 = Adaptor3d_HSurfaceTool::LastUParameter(Caro2);
371   VM2 = Adaptor3d_HSurfaceTool::LastVParameter(Caro2);
372
373   ResoU1 = Adaptor3d_HSurfaceTool::UResolution(Caro1,Precision::Confusion());
374   ResoV1 = Adaptor3d_HSurfaceTool::VResolution(Caro1,Precision::Confusion());
375
376   ResoU2 = Adaptor3d_HSurfaceTool::UResolution(Caro2,Precision::Confusion());
377   ResoV2 = Adaptor3d_HSurfaceTool::VResolution(Caro2,Precision::Confusion());
378   //
379   Standard_Real NEWRESO, MAXVAL, MAXVAL2;
380   //
381   MAXVAL  = Abs(Um1);  
382   MAXVAL2 = Abs(UM1);
383   if(MAXVAL2 > MAXVAL) {
384     MAXVAL = MAXVAL2;
385   }
386   NEWRESO = ResoU1 * MAXVAL ;
387   if(NEWRESO > ResoU1) {
388     ResoU1 = NEWRESO;  
389   }
390   //
391   MAXVAL  = Abs(Um2);   
392   MAXVAL2 = Abs(UM2);
393   if(MAXVAL2 > MAXVAL){
394     MAXVAL = MAXVAL2;
395   }  
396   NEWRESO = ResoU2 * MAXVAL ;
397   if(NEWRESO > ResoU2) {
398     ResoU2 = NEWRESO;  
399   }
400   //
401   MAXVAL  = Abs(Vm1);  
402   MAXVAL2 = Abs(VM1);
403   if(MAXVAL2 > MAXVAL) {
404     MAXVAL = MAXVAL2;
405   }
406   NEWRESO = ResoV1 * MAXVAL ;
407   if(NEWRESO > ResoV1) {    
408     ResoV1 = NEWRESO; 
409   }
410   //
411   MAXVAL  = Abs(Vm2);  
412   MAXVAL2 = Abs(VM2);
413   if(MAXVAL2 > MAXVAL){
414     MAXVAL = MAXVAL2;
415   }  
416   NEWRESO = ResoV2 * MAXVAL ;
417   if(NEWRESO > ResoV2) {  
418     ResoV2 = NEWRESO;
419   }
420   //
421   pasuv[0]=pasMax*Abs(UM1-Um1);
422   pasuv[1]=pasMax*Abs(VM1-Vm1);
423   pasuv[2]=pasMax*Abs(UM2-Um2);
424   pasuv[3]=pasMax*Abs(VM2-Vm2);
425   //
426   if(Adaptor3d_HSurfaceTool::IsUPeriodic(Caro1)==Standard_False) { 
427     UM1+=KELARG*pasuv[0];  
428     Um1-=KELARG*pasuv[0];
429   }
430   else { 
431     Standard_Real t = UM1-Um1; 
432     if(t<Adaptor3d_HSurfaceTool::UPeriod(Caro1)) { 
433       t=0.5*(Adaptor3d_HSurfaceTool::UPeriod(Caro1)-t);
434       t=(t>KELARG*pasuv[0])? KELARG*pasuv[0] : t;
435       UM1+=t;  
436       Um1-=t;
437     }
438   }
439   //
440   if(Adaptor3d_HSurfaceTool::IsVPeriodic(Caro1)==Standard_False) { 
441     VM1+=KELARG*pasuv[1];
442     Vm1-=KELARG*pasuv[1];
443   }
444   else { 
445     Standard_Real t = VM1-Vm1; 
446     if(t<Adaptor3d_HSurfaceTool::VPeriod(Caro1)) { 
447       t=0.5*(Adaptor3d_HSurfaceTool::VPeriod(Caro1)-t);
448       t=(t>KELARG*pasuv[1])? KELARG*pasuv[1] : t;
449       VM1+=t;  Vm1-=t;
450     }
451   }
452   //
453   if(Adaptor3d_HSurfaceTool::IsUPeriodic(Caro2)==Standard_False) { 
454     UM2+=KELARG*pasuv[2];  
455     Um2-=KELARG*pasuv[2];
456   }
457   else { 
458     Standard_Real t = UM2-Um2; 
459     if(t<Adaptor3d_HSurfaceTool::UPeriod(Caro2)) { 
460       t=0.5*(Adaptor3d_HSurfaceTool::UPeriod(Caro2)-t);
461       t=(t>KELARG*pasuv[2])? KELARG*pasuv[2] : t;
462       UM2+=t;  
463       Um2-=t;
464     }
465   }
466
467   if(Adaptor3d_HSurfaceTool::IsVPeriodic(Caro2)==Standard_False) {   
468     VM2+=KELARG*pasuv[3];  
469     Vm2-=KELARG*pasuv[3];
470   }
471   else { 
472     Standard_Real t = VM2-Vm2; 
473     if(t<Adaptor3d_HSurfaceTool::VPeriod(Caro2)) { 
474       t=0.5*(Adaptor3d_HSurfaceTool::VPeriod(Caro2)-t);
475       t=(t>KELARG*pasuv[3])? KELARG*pasuv[3] : t;
476       VM2+=t;  
477       Vm2-=t;
478     }
479   }
480   //-- ComputePasInit(pasuv,Um1,UM1,Vm1,VM1,Um2,UM2,Vm2,VM2,Caro1,Caro2);
481
482   for (Standard_Integer i = 0; i<=3;i++) {
483     pasInit[i] = pasSav[i] = pasuv[i]; 
484   }  
485
486   if(ResoU1>0.0001*pasuv[0]) ResoU1=0.00001*pasuv[0];
487   if(ResoV1>0.0001*pasuv[1]) ResoV1=0.00001*pasuv[1];
488   if(ResoU2>0.0001*pasuv[2]) ResoU2=0.00001*pasuv[2];
489   if(ResoV2>0.0001*pasuv[3]) ResoV2=0.00001*pasuv[3];
490   
491   myStepMin[0] = 100.0*ResoU1;
492   myStepMin[1] = 100.0*ResoV1;
493   myStepMin[2] = 100.0*ResoU2;
494   myStepMin[3] = 100.0*ResoV2;
495
496   //
497   TColStd_Array1OfReal Par(1,4);
498   Par(1) = U1;
499   Par(2) = V1;
500   Par(3) = U2;
501   Par(4) = V2;
502   Perform(Par);
503 }
504
505 //==================================================================================
506 // function : PerformFirstPoint
507 // purpose  : 
508 //==================================================================================
509 Standard_Boolean IntWalk_PWalking::PerformFirstPoint  (const TColStd_Array1OfReal& ParDep,
510                                                        IntSurf_PntOn2S& FirstPoint)   
511 {
512   sensCheminement = 1;
513   close = Standard_False;
514   //
515   Standard_Integer i;
516   TColStd_Array1OfReal Param(1,4);
517   //
518   for (i=1; i<=4; ++i) {
519     Param(i) = ParDep(i);
520   }
521   //-- calculate the first solution point
522   math_FunctionSetRoot  Rsnld(myIntersectionOn2S.Function());
523   //
524   myIntersectionOn2S.Perform(Param,Rsnld);
525   if (!myIntersectionOn2S.IsDone())  { 
526     return Standard_False;
527   }
528
529   if (myIntersectionOn2S.IsEmpty()) {
530     return Standard_False;
531   }
532
533   FirstPoint = myIntersectionOn2S.Point();
534   return Standard_True;
535 }
536 //==================================================================================
537 // function : Perform
538 // purpose  : 
539 //==================================================================================
540 void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep)    
541 {
542   Perform(ParDep,Um1,Vm1,Um2,Vm2,UM1,VM1,UM2,VM2);
543 }
544
545 //=======================================================================
546 //function : SQDistPointSurface
547 //purpose  : Returns square distance between thePnt and theSurf.
548 //            (theU0, theV0) is initial point for extrema
549 //=======================================================================
550 static Standard_Real SQDistPointSurface(const gp_Pnt &thePnt,
551                                         const Adaptor3d_Surface& theSurf,
552                                         const Standard_Real theU0,
553                                         const Standard_Real theV0)
554 {
555   Extrema_GenLocateExtPS aExtPS(theSurf);
556   aExtPS.Perform(thePnt, theU0, theV0);
557
558   if(!aExtPS.IsDone())
559     return RealLast();
560   
561   return aExtPS.SquareDistance();
562 }
563
564 //==================================================================================
565 // function : IsTangentExtCheck
566 // purpose  : Additional check if the surfaces are tangent.
567 //            Checks if any point in one surface lie in another surface
568 //            (with given tolerance)
569 //==================================================================================
570 static Standard_Boolean IsTangentExtCheck(const Handle(Adaptor3d_HSurface)& theSurf1,
571                                           const Handle(Adaptor3d_HSurface)& theSurf2,
572                                           const Standard_Real theU10,
573                                           const Standard_Real theV10,
574                                           const Standard_Real theU20,
575                                           const Standard_Real theV20,
576                                           const Standard_Real theToler,
577                                           const Standard_Real theArrStep[])
578 {
579   {
580     gp_Pnt aPt;
581     gp_Vec aDu1, aDv1, aDu2, aDv2;
582     theSurf1->D1(theU10, theV10, aPt, aDu1, aDv1);
583     theSurf2->D1(theU20, theV20, aPt, aDu2, aDv2);
584
585     const gp_Vec  aN1(aDu1.Crossed(aDv1)),
586                   aN2(aDu2.Crossed(aDv2));
587     const Standard_Real aDP = aN1.Dot(aN2),
588                         aSQ1 = aN1.SquareMagnitude(),
589                         aSQ2 = aN2.SquareMagnitude();
590
591     if((aSQ1 < RealSmall()) || (aSQ2 < RealSmall()))
592       return Standard_True; //Tangent
593
594     if(aDP*aDP < 0.9998*aSQ1*aSQ2)
595     {//cos(ang N1<->N2) < 0.9999
596       return Standard_False; //Not tangent
597     }
598   }
599
600   //For two faces (2^2 = 4)
601   const Standard_Real aSQToler = 4.0*theToler*theToler;
602   const Standard_Integer aNbItems = 4;
603   const Standard_Real aParUS1[aNbItems] = { theU10 + theArrStep[0],
604                                             theU10 - theArrStep[0],
605                                             theU10, theU10};
606   const Standard_Real aParVS1[aNbItems] = { theV10, theV10,
607                                             theV10 + theArrStep[1],
608                                             theV10 - theArrStep[1]};
609   const Standard_Real aParUS2[aNbItems] = { theU20 + theArrStep[2],
610                                             theU20 - theArrStep[2],
611                                             theU20, theU20};
612   const Standard_Real aParVS2[aNbItems] = { theV20, theV20,
613                                             theV20 + theArrStep[3],
614                                             theV20 - theArrStep[3]};
615
616   for(Standard_Integer i = 0; i < aNbItems; i++)
617   {
618     gp_Pnt aP(theSurf1->Value(aParUS1[i], aParVS1[i]));
619     const Standard_Real aSqDist = SQDistPointSurface(aP, theSurf2->Surface(), theU20, theV20);
620     if(aSqDist > aSQToler)
621       return Standard_False;
622   }
623
624   for(Standard_Integer i = 0; i < aNbItems; i++)
625   {
626     gp_Pnt aP(theSurf2->Value(aParUS2[i], aParVS2[i]));
627     const Standard_Real aSqDist = SQDistPointSurface(aP, theSurf1->Surface(), theU10, theV10);
628     if(aSqDist > aSQToler)
629       return Standard_False;
630   }
631
632   return Standard_True;
633 }
634
635 //==================================================================================
636 // function : Perform
637 // purpose  : 
638 //==================================================================================
639 void IntWalk_PWalking::Perform(const TColStd_Array1OfReal& ParDep,
640                                const Standard_Real u1min,
641                                const Standard_Real v1min,
642                                const Standard_Real u2min,
643                                const Standard_Real v2min,
644                                const Standard_Real u1max,
645                                const Standard_Real v1max,
646                                const Standard_Real u2max,
647                                const Standard_Real v2max)
648 {
649   const Standard_Real aSQDistMax = 1.0e-14;
650   //xf
651
652   Standard_Integer NbPasOKConseq=0;
653   TColStd_Array1OfReal Param(1,4);
654   IntImp_ConstIsoparametric ChoixIso;
655   //xt
656   //
657   done = Standard_False;
658   //
659   // Caro1 and Caro2
660   const Handle(Adaptor3d_HSurface)& Caro1 =myIntersectionOn2S.Function().AuxillarSurface1();
661   const Handle(Adaptor3d_HSurface)& Caro2 =myIntersectionOn2S.Function().AuxillarSurface2();
662   //
663   const Standard_Real UFirst1 = Adaptor3d_HSurfaceTool::FirstUParameter(Caro1);
664   const Standard_Real VFirst1 = Adaptor3d_HSurfaceTool::FirstVParameter(Caro1);
665   const Standard_Real ULast1  = Adaptor3d_HSurfaceTool::LastUParameter (Caro1);
666   const Standard_Real VLast1  = Adaptor3d_HSurfaceTool::LastVParameter (Caro1);
667
668   const Standard_Real UFirst2 = Adaptor3d_HSurfaceTool::FirstUParameter(Caro2);
669   const Standard_Real VFirst2 = Adaptor3d_HSurfaceTool::FirstVParameter(Caro2);
670   const Standard_Real ULast2  = Adaptor3d_HSurfaceTool::LastUParameter (Caro2);
671   const Standard_Real VLast2  = Adaptor3d_HSurfaceTool::LastVParameter (Caro2);
672   //
673   ComputePasInit(u1max - u1min,v1max - v1min,u2max - u2min,v2max - v2min);
674
675   for (Standard_Integer i=0; i<4; ++i)
676   {
677     if(pasuv[i]>10)
678     {
679       pasuv[i] = 10;
680     }
681
682     pasInit[i] = pasSav[i] = pasuv[i]; 
683   }
684   //
685   line = new IntSurf_LineOn2S ();
686   //
687   for (Standard_Integer i=1; i<=4; ++i)
688   {
689     Param(i)=ParDep(i);
690   }
691   //-- reproduce steps uv connected to surfaces Caro1 and Caro2
692   //-- pasuv[] and pasSav[] are modified during the marching
693   for(Standard_Integer i = 0; i < 4; ++i)
694   {
695     pasSav[i] = pasuv[i] = pasInit[i];
696   }
697
698   //-- calculate the first solution point
699   math_FunctionSetRoot  Rsnld(myIntersectionOn2S.Function());
700   //
701   ChoixIso = myIntersectionOn2S.Perform(Param,Rsnld);
702   if (!myIntersectionOn2S.IsDone())
703   {
704     return;
705   }
706
707   //
708   if (myIntersectionOn2S.IsEmpty())
709   {
710     return;
711   }
712   //
713   if(myIntersectionOn2S.IsTangent())
714   {
715     return;
716   }
717   //
718   Standard_Boolean Arrive, DejaReparti;
719   const Standard_Integer RejectIndexMAX = 250000;
720   Standard_Integer IncKey, RejectIndex;
721   gp_Pnt pf,pl;
722   //
723   DejaReparti = Standard_False;
724   IncKey = 0;
725   RejectIndex = 0;
726   //
727   previousPoint = myIntersectionOn2S.Point();
728   previoustg = Standard_False;
729   previousd  = myIntersectionOn2S.Direction();
730   previousd1 = myIntersectionOn2S.DirectionOnS1();
731   previousd2 = myIntersectionOn2S.DirectionOnS2();
732   myTangentIdx = 1;
733   tgdir   = previousd;
734   firstd1 = previousd1;
735   firstd2 = previousd2;
736   tgfirst = tglast = Standard_False;
737   choixIsoSav  =  ChoixIso;
738   //------------------------------------------------------------
739   //-- Test if the first point of marching corresponds 
740   //-- to a point on borders. 
741   //-- In this case, DejaReparti is initialized as True
742   //-- 
743   pf = previousPoint.Value();
744   Standard_Boolean bTestFirstPoint = Standard_True;
745
746   previousPoint.Parameters(Param(1),Param(2),Param(3),Param(4));
747
748   if(IsTangentExtCheck(Caro1, Caro2, Param(1), Param(2), Param(3), Param(4), myTolTang, pasuv))
749     return;
750
751   AddAPoint(previousPoint);
752   //
753   IntWalk_StatusDeflection aStatus = IntWalk_OK, aPrevStatus = IntWalk_OK;
754   Standard_Boolean NoTestDeflection = Standard_False;
755   Standard_Real SvParam[4], f;
756   Standard_Integer LevelOfEmptyInmyIntersectionOn2S=0;
757   Standard_Integer LevelOfPointConfondu = 0; 
758   Standard_Integer LevelOfIterWithoutAppend = -1;
759   //
760
761   const Standard_Real aTol[4] = { Epsilon(u1max - u1min),
762                                   Epsilon(v1max - v1min),
763                                   Epsilon(u2max - u2min),
764                                   Epsilon(v2max - v2min)};
765   Arrive = Standard_False;
766   while(!Arrive) //010
767   {
768     aPrevStatus = aStatus;
769
770     LevelOfIterWithoutAppend++;
771     if(LevelOfIterWithoutAppend>20)
772     {
773       Arrive = Standard_True; 
774       if(DejaReparti) {
775         break;
776       }
777       RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
778       LevelOfIterWithoutAppend = 0;
779     }
780     //
781     // compute f
782     f = 0.;
783     switch (ChoixIso) { 
784       case IntImp_UIsoparametricOnCaro1: f = Abs(previousd1.X()); break;
785       case IntImp_VIsoparametricOnCaro1: f = Abs(previousd1.Y()); break;
786       case IntImp_UIsoparametricOnCaro2: f = Abs(previousd2.X()); break;
787       case IntImp_VIsoparametricOnCaro2: f = Abs(previousd2.Y()); break;
788       default:break;
789     }
790     //
791     if(f<0.1) {
792       f=0.1;
793     }
794     //
795     previousPoint.Parameters(Param(1),Param(2),Param(3),Param(4));
796     //
797     //--ofv.begin
798     Standard_Real aIncKey, aEps, dP1, dP2, dP3, dP4;
799     //
800     dP1 = sensCheminement * pasuv[0] * previousd1.X() /f;
801     dP2 = sensCheminement * pasuv[1] * previousd1.Y() /f;
802     dP3 = sensCheminement * pasuv[2] * previousd2.X() /f; 
803     dP4 = sensCheminement * pasuv[3] * previousd2.Y() /f;
804     //
805     aIncKey=5.*(Standard_Real)IncKey;
806     aEps=1.e-7;
807     if(ChoixIso == IntImp_UIsoparametricOnCaro1 && Abs(dP1) < aEps)
808     {
809       dP1 *= aIncKey;
810     }
811
812     if(ChoixIso == IntImp_VIsoparametricOnCaro1 && Abs(dP2) < aEps)
813     {
814       dP2 *= aIncKey;
815     }
816
817     if(ChoixIso == IntImp_UIsoparametricOnCaro2 && Abs(dP3) < aEps)
818     {
819       dP3 *= aIncKey;
820     }
821
822     if(ChoixIso == IntImp_VIsoparametricOnCaro2 && Abs(dP4) < aEps)
823     {
824       dP4 *= aIncKey;
825     }
826     //--ofv.end
827     //
828     Param(1) += dP1;
829     Param(2) += dP2;
830     Param(3) += dP3; 
831     Param(4) += dP4;
832     //==========================
833     SvParam[0]=Param(1); 
834     SvParam[1]=Param(2);
835     SvParam[2]=Param(3);
836     SvParam[3]=Param(4);
837     //
838     Standard_Integer aTryNumber = 0;
839     Standard_Boolean isBadPoint = Standard_False;
840     IntImp_ConstIsoparametric aBestIso = ChoixIso;
841     do
842     {
843       isBadPoint = Standard_False;
844
845       ChoixIso= myIntersectionOn2S.Perform(Param, Rsnld, aBestIso);
846
847       if (myIntersectionOn2S.IsDone() && !myIntersectionOn2S.IsEmpty())
848       {
849         //If we go along any surface boundary then it is possible 
850         //to find "outboundaried" point.
851         //Nevertheless, if this deflection is quite small, we will be
852         //able to adjust this point to the boundary.
853
854         Standard_Real aNewPnt[4], anAbsParamDist[4];
855         myIntersectionOn2S.Point().Parameters(aNewPnt[0], aNewPnt[1], aNewPnt[2], aNewPnt[3]);
856         const Standard_Real aParMin[4] = {u1min, v1min, u2min, v2min};
857         const Standard_Real aParMax[4] = {u1max, v1max, u2max, v2max};
858
859         for(Standard_Integer i = 0; i < 4; i++)
860         {
861           if(Abs(aNewPnt[i] - aParMin[i]) < aTol[i])
862             aNewPnt[i] = aParMin[i];
863           else if(Abs(aNewPnt[i] - aParMax[i]) < aTol[i])
864             aNewPnt[i] = aParMax[i];
865         }
866
867         if (aNewPnt[0] < u1min || aNewPnt[0] > u1max ||
868             aNewPnt[1] < v1min || aNewPnt[1] > v1max ||
869             aNewPnt[2] < u2min || aNewPnt[2] > u2max ||
870             aNewPnt[3] < v2min || aNewPnt[3] > v2max)
871         {
872           break; // Out of borders, handle this later.
873         }
874
875         myIntersectionOn2S.ChangePoint().SetValue(aNewPnt[0],
876                                                   aNewPnt[1],
877                                                   aNewPnt[2],
878                                                   aNewPnt[3]);
879
880         anAbsParamDist[0] = Abs(Param(1) - dP1 - aNewPnt[0]);
881         anAbsParamDist[1] = Abs(Param(2) - dP2 - aNewPnt[1]);
882         anAbsParamDist[2] = Abs(Param(3) - dP3 - aNewPnt[2]);
883         anAbsParamDist[3] = Abs(Param(4) - dP4 - aNewPnt[3]);
884         if (anAbsParamDist[0] < ResoU1 &&
885             anAbsParamDist[1] < ResoV1 &&
886             anAbsParamDist[2] < ResoU2 &&
887             anAbsParamDist[3] < ResoV2 &&
888             aStatus != IntWalk_PasTropGrand)
889         {
890           isBadPoint = Standard_True;
891           aBestIso = IntImp_ConstIsoparametric((aBestIso + 1) % 4);
892         }
893       }
894     } while (isBadPoint && ++aTryNumber <= 4);
895     //
896     if (!myIntersectionOn2S.IsDone())
897     {
898       //end of line, division
899       Arrive = Standard_False;
900       Param(1)=SvParam[0]; 
901       Param(2)=SvParam[1]; 
902       Param(3)=SvParam[2];
903       Param(4)=SvParam[3];
904       RepartirOuDiviser(DejaReparti, ChoixIso, Arrive);
905     }
906     else  //009 
907     {
908       //== Calculation of exact point from Param(.) is possible
909       if (myIntersectionOn2S.IsEmpty())
910       {
911         Standard_Real u1,v1,u2,v2;
912         previousPoint.Parameters(u1,v1,u2,v2);
913         //
914         Arrive = Standard_False;
915         if(u1<UFirst1 || u1>ULast1)
916         {
917           Arrive=Standard_True;
918         }       
919
920         if(u2<UFirst2 || u2>ULast2)
921         {
922           Arrive=Standard_True;
923         }
924
925         if(v1<VFirst1 || v1>VLast1)
926         {
927           Arrive=Standard_True;
928         }
929
930         if(v2<VFirst2 || v2>VLast2)
931         {
932           Arrive=Standard_True;
933         }
934
935         RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
936         LevelOfEmptyInmyIntersectionOn2S++;
937         //
938         if(LevelOfEmptyInmyIntersectionOn2S>10)
939         {
940           pasuv[0]=pasSav[0]; 
941           pasuv[1]=pasSav[1]; 
942           pasuv[2]=pasSav[2]; 
943           pasuv[3]=pasSav[3];
944         }
945       }
946       else //008
947       {
948         //============================================================
949         //== A point has been found :  T E S T   D E F L E C T I O N 
950         //============================================================
951         if(NoTestDeflection)
952         {
953           NoTestDeflection = Standard_False;
954         }                 
955         else
956         {
957           if(--LevelOfEmptyInmyIntersectionOn2S<=0)
958           {
959             LevelOfEmptyInmyIntersectionOn2S=0;
960             if(LevelOfIterWithoutAppend < 10)
961             {
962               aStatus = TestDeflection(ChoixIso);
963             }                   
964             else
965             {
966               pasuv[0]*=0.5; 
967               pasuv[1]*=0.5; 
968               pasuv[2]*=0.5; 
969               pasuv[3]*=0.5;
970             }
971           }
972         }
973
974         //============================================================
975         //==       T r a i t e m e n t   s u r   S t a t u s        ==
976         //============================================================
977         if(LevelOfPointConfondu > 5)
978         { 
979           aStatus = IntWalk_ArretSurPoint;
980           LevelOfPointConfondu = 0;  
981         }
982         //
983         if(aStatus==IntWalk_OK)
984         { 
985           NbPasOKConseq++;
986           if(NbPasOKConseq >= 5)
987           {
988             NbPasOKConseq=0;
989             Standard_Boolean pastroppetit;
990             Standard_Real t;
991             //
992             do
993             {
994               pastroppetit=Standard_True;
995               //
996               if(pasuv[0]<pasInit[0])
997               {
998                 t = (pasInit[0]-pasuv[0])*0.25;
999                 if(t>0.1*pasInit[0])
1000                 {
1001                   t=0.1*pasuv[0];
1002                 }
1003
1004                 pasuv[0]+=t; 
1005                 pastroppetit=Standard_False;
1006               }
1007
1008               if(pasuv[1]<pasInit[1])
1009               {
1010                 t = (pasInit[1]-pasuv[1])*0.25;
1011                 if(t>0.1*pasInit[1]) {
1012                   t=0.1*pasuv[1];
1013                 }               
1014
1015                 pasuv[1]+=t; 
1016                 pastroppetit=Standard_False;
1017               }
1018
1019               if(pasuv[2]<pasInit[2])
1020               {
1021                 t = (pasInit[2]-pasuv[2])*0.25;
1022                 if(t>0.1*pasInit[2])
1023                 {
1024                   t=0.1*pasuv[2];
1025                 }
1026
1027                 pasuv[2]+=t; 
1028                 pastroppetit=Standard_False;
1029               }
1030
1031               if(pasuv[3]<pasInit[3])
1032               {
1033                 t = (pasInit[3]-pasuv[3])*0.25;
1034                 if(t>0.1*pasInit[3]) {
1035                   t=0.1*pasuv[3];
1036                 }
1037                 pasuv[3]+=t; 
1038                 pastroppetit=Standard_False;
1039               }
1040               if(pastroppetit)
1041               {
1042                 if(pasMax<0.1)
1043                 {
1044                   pasMax*=1.1;
1045                   pasInit[0]*=1.1; 
1046                   pasInit[1]*=1.1; 
1047                   pasInit[2]*=1.1; 
1048                   pasInit[3]*=1.1; 
1049                 }
1050                 else
1051                 {
1052                   pastroppetit=Standard_False;
1053                 }
1054               }
1055             }
1056             while(pastroppetit);
1057           }
1058         }//aStatus==IntWalk_OK
1059         else
1060           NbPasOKConseq=0;
1061
1062         //
1063         switch(aStatus)//007
1064         {
1065         case IntWalk_ArretSurPointPrecedent:
1066           {
1067             Arrive = Standard_False;
1068             RepartirOuDiviser(DejaReparti, ChoixIso, Arrive);
1069             break;
1070           }
1071         case IntWalk_PasTropGrand:
1072           {
1073             Param(1)=SvParam[0];
1074             Param(2)=SvParam[1]; 
1075             Param(3)=SvParam[2]; 
1076             Param(4)=SvParam[3];
1077
1078             if(LevelOfIterWithoutAppend > 5)
1079             {
1080               for (Standard_Integer i = 0; i < 4; i++)
1081               {
1082                 if (pasSav[i] > pasInit[i])
1083                   continue;
1084
1085                 const Standard_Real aDelta = (pasInit[i]-pasSav[i])*0.25;
1086
1087                 if(aDelta > Epsilon(pasInit[i]))
1088                 {
1089                   pasInit[i] -= aDelta;
1090                   LevelOfIterWithoutAppend=0;
1091                 }
1092               }
1093             }
1094
1095             break;
1096           }
1097         case IntWalk_PointConfondu:
1098           {
1099             LevelOfPointConfondu++;
1100
1101             if(LevelOfPointConfondu>5)
1102             {
1103               Standard_Boolean pastroppetit;
1104               //
1105               do
1106               {
1107                 pastroppetit=Standard_True;
1108
1109                 for(Standard_Integer i = 0; i < 4; i++)
1110                 {
1111                   if(pasuv[i]<pasInit[i])
1112                   {
1113                     pasuv[i]+=(pasInit[i]-pasuv[i])*0.25;
1114                     pastroppetit=Standard_False;
1115                   }
1116                 }
1117
1118                 if(pastroppetit)
1119                 {
1120                   if(pasMax<0.1)
1121                   {
1122                     pasMax*=1.1;
1123                     pasInit[0]*=1.1;
1124                     pasInit[1]*=1.1;
1125                     pasInit[2]*=1.1;
1126                     pasInit[3]*=1.1; 
1127                   }
1128                   else
1129                   {
1130                     pastroppetit=Standard_False;
1131                   }
1132                 }
1133               }
1134               while(pastroppetit);
1135             }
1136
1137             break;
1138           }
1139         case IntWalk_StepTooSmall:
1140           {
1141             Standard_Boolean hasStepBeenIncreased = Standard_False;
1142
1143             for(Standard_Integer i = 0; i < 4; i++)
1144             {
1145               const Standard_Real aNewStep = Min(1.5*pasuv[i], pasInit[i]);
1146               if(aNewStep > pasuv[i])
1147               {
1148                 pasuv[i] = aNewStep;
1149                 hasStepBeenIncreased = Standard_True;
1150               }
1151             }
1152
1153             if(hasStepBeenIncreased)
1154             {
1155               Param(1)=SvParam[0];
1156               Param(2)=SvParam[1];
1157               Param(3)=SvParam[2];
1158               Param(4)=SvParam[3];
1159
1160               // In order to avoid cyclic changes
1161               // (PasTropGrand --> Decrease step --> 
1162               // StepTooSmall --> Increase step --> PasTropGrand...)
1163               // nullify LevelOfIterWithoutAppend only if the condition
1164               // is satisfied:
1165               if (aPrevStatus != IntWalk_PasTropGrand)
1166                 LevelOfIterWithoutAppend = 0;
1167
1168               break;
1169             }
1170           }
1171           Standard_FALLTHROUGH
1172         case IntWalk_OK:
1173         case IntWalk_ArretSurPoint://006
1174           {
1175             //=======================================================
1176             //== Stop Test t   :  Frame on Param(.)     ==
1177             //=======================================================
1178             //xft arrive here
1179             Arrive = TestArret(DejaReparti,Param,ChoixIso); 
1180             // JMB 30th December 1999. 
1181             // Some statement below should not be put in comment because they are useful.
1182             // See grid CTO 909 A1 which infinitely loops 
1183             if(Arrive==Standard_False && aStatus==IntWalk_ArretSurPoint)
1184             {
1185               Arrive=Standard_True;
1186 #ifdef OCCT_DEBUG
1187               cout << "IntWalk_PWalking_1.gxx: Problems with intersection"<<endl;
1188 #endif
1189             }
1190
1191             if(Arrive)
1192             {
1193               NbPasOKConseq = -10;
1194             }
1195
1196             if(!Arrive)//005
1197             {
1198               //=====================================================
1199               //== Param(.) is in the limits                       ==
1200               //==  and does not end a closed  line                ==
1201               //=====================================================
1202               //== Check on the current point of myInters
1203               Standard_Boolean pointisvalid = Standard_False;
1204               {
1205                 Standard_Real u1,v1,u2,v2; 
1206                 myIntersectionOn2S.Point().Parameters(u1,v1,u2,v2);
1207
1208                 //
1209                 if(u1 <= UM1  && u2 <= UM2 && v1 <= VM1 && 
1210                   v2 <= VM2  && u1 >= Um1 && u2 >= Um2 &&
1211                   v1 >= Vm1  && v2 >= Vm2)
1212                 {
1213                   pointisvalid=Standard_True;
1214                 }
1215               }
1216
1217               //
1218               if(pointisvalid)
1219               {
1220                 previousPoint = myIntersectionOn2S.Point();
1221                 previoustg = myIntersectionOn2S.IsTangent();
1222
1223                 if(!previoustg)
1224                 {
1225                   previousd  = myIntersectionOn2S.Direction();
1226                   previousd1 = myIntersectionOn2S.DirectionOnS1();
1227                   previousd2 = myIntersectionOn2S.DirectionOnS2();
1228                 }
1229                 //=====================================================
1230                 //== Check on the previous Point
1231                 {
1232                   Standard_Real u1,v1,u2,v2;
1233                   previousPoint.Parameters(u1,v1,u2,v2); 
1234                   if( u1 <= UM1  && u2 <= UM2 && v1 <= VM1 &&
1235                     v2 <= VM2  && u1 >= Um1 && u2 >= Um2 &&
1236                     v1 >= Vm1  && v2 >= Vm2)
1237                   {
1238                     pl = previousPoint.Value();
1239                     if(bTestFirstPoint)
1240                     {
1241                       if(pf.SquareDistance(pl) < aSQDistMax)
1242                       {
1243                         IncKey++;
1244                         if(IncKey == 5000)
1245                           return;
1246                         else
1247                           continue;
1248                       }
1249                       else
1250                       {
1251                         bTestFirstPoint = Standard_False;
1252                       }
1253                     }
1254                     //
1255                     AddAPoint(previousPoint);
1256                     RejectIndex++;
1257
1258                     if(RejectIndex >= RejectIndexMAX)
1259                     {
1260                       Arrive = Standard_True;
1261                       break;
1262                     }
1263
1264                     //
1265                     LevelOfIterWithoutAppend = 0;
1266                   }
1267                 }
1268               }//pointisvalid
1269               //====================================================
1270
1271               if (aStatus == IntWalk_ArretSurPoint)
1272               {
1273                 RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
1274               }
1275               else
1276               {
1277                 if (line->NbPoints() == 2)
1278                 {
1279                   pasSav[0] = pasuv[0];
1280                   pasSav[1] = pasuv[1];
1281                   pasSav[2] = pasuv[2];
1282                   pasSav[3] = pasuv[3];
1283                   
1284                   if ((aPrevStatus == IntWalk_PasTropGrand) &&
1285                       (LevelOfIterWithoutAppend > 0))
1286                   {
1287                     pasInit[0] = pasuv[0];
1288                     pasInit[1] = pasuv[1];
1289                     pasInit[2] = pasuv[2];
1290                     pasInit[3] = pasuv[3];
1291                   }
1292                 }
1293               }
1294             }//005 if(!Arrive)
1295             else  //004
1296             {
1297               if(close)
1298               {
1299                 //================= la ligne est fermee ===============
1300                 AddAPoint(line->Value(1)); //ligne fermee
1301                 LevelOfIterWithoutAppend=0;
1302               }
1303               else    //$$$
1304               {
1305                 //====================================================
1306                 //== Param was not in the limits (was reframed)
1307                 //====================================================
1308                 Standard_Boolean bPrevNotTangent = !previoustg || !myIntersectionOn2S.IsTangent();
1309
1310                 IntImp_ConstIsoparametric SauvChoixIso = ChoixIso;
1311                 ChoixIso = myIntersectionOn2S.Perform(Param,Rsnld,ChoixIso);
1312                 //
1313                 if(!myIntersectionOn2S.IsEmpty()) //002
1314                 {
1315                   // mutially outpasses in the square or intersection in corner
1316
1317                   if(TestArret(Standard_True,Param,ChoixIso))
1318                   {
1319                     NbPasOKConseq = -10;
1320                     ChoixIso = myIntersectionOn2S.Perform(Param,Rsnld,ChoixIso);
1321
1322                     if(!myIntersectionOn2S.IsEmpty())
1323                     {
1324                       previousPoint = myIntersectionOn2S.Point();
1325                       previoustg = myIntersectionOn2S.IsTangent();
1326
1327                       if (!previoustg)
1328                       {
1329                         previousd  = myIntersectionOn2S.Direction();
1330                         previousd1 = myIntersectionOn2S.DirectionOnS1();
1331                         previousd2 = myIntersectionOn2S.DirectionOnS2();
1332                       }
1333
1334                       pl = previousPoint.Value();
1335
1336                       if(bTestFirstPoint)
1337                       {
1338                         if(pf.SquareDistance(pl) < aSQDistMax)
1339                         {
1340                           IncKey++;
1341                           if(IncKey == 5000)
1342                             return;
1343                           else
1344                             continue;
1345                         }
1346                         else
1347                         {
1348                           bTestFirstPoint = Standard_False;
1349                         }
1350                       }
1351                       //
1352                       AddAPoint(previousPoint);
1353                       RejectIndex++;
1354
1355                       if(RejectIndex >= RejectIndexMAX)
1356                       {
1357                         Arrive = Standard_True;
1358                         break;
1359                       }
1360
1361                       //
1362                       LevelOfIterWithoutAppend=0;
1363                       RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
1364                     }
1365                     else
1366                     {
1367                       //fail framing divides the step
1368                       Arrive = Standard_False;
1369                       RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
1370                       NoTestDeflection = Standard_True;
1371                       ChoixIso = SauvChoixIso;
1372                     }
1373                   }//if(TestArret())
1374                   else
1375                   {
1376                     // save the last point
1377                     // to revert to it if the current point is out of bounds
1378
1379                     IntSurf_PntOn2S previousPointSave = previousPoint;
1380                     Standard_Boolean previoustgSave   = previoustg;
1381                     gp_Dir previousdSave              = previousd;
1382                     gp_Dir2d previousd1Save           = previousd1;
1383                     gp_Dir2d previousd2Save           = previousd2;
1384
1385                     previousPoint = myIntersectionOn2S.Point();
1386                     previoustg = myIntersectionOn2S.IsTangent();
1387                     Arrive = Standard_False;
1388
1389                     if(!previoustg)
1390                     {
1391                       previousd  = myIntersectionOn2S.Direction();
1392                       previousd1 = myIntersectionOn2S.DirectionOnS1();
1393                       previousd2 = myIntersectionOn2S.DirectionOnS2();
1394                     }
1395
1396                     //========================================
1397                     //== Check on PreviousPoint @@
1398
1399                     {
1400                       Standard_Real u1,v1,u2,v2;
1401                       previousPoint.Parameters(u1,v1,u2,v2);
1402
1403                       //To save initial 2d points
1404                       gp_Pnt2d ParamPntOnS1(Param(1), Param(2));
1405                       gp_Pnt2d ParamPntOnS2(Param(3), Param(4));
1406
1407                       ///////////////////////////
1408                       Param(1) = u1;
1409                       Param(2) = v1;
1410                       Param(3) = u2;
1411                       Param(4) = v2;
1412                       //
1413
1414                       //xf
1415                       Standard_Boolean bFlag1, bFlag2;
1416                       Standard_Real aTol2D=1.e-11;
1417                       //
1418                       bFlag1=u1 >= Um1-aTol2D && v1 >= Vm1-aTol2D && u1 <= UM1+aTol2D && v1 <= VM1+aTol2D;
1419                       bFlag2=u2 >= Um2-aTol2D && v2 >= Vm2-aTol2D && u2 <= UM2+aTol2D && v2 <= VM2+aTol2D;
1420                       if (bFlag1 && bFlag2)
1421                       {
1422                         if (line->NbPoints() > 1)
1423                         {
1424                           IntSurf_PntOn2S prevprevPoint = line->Value(line->NbPoints()-1);
1425                           Standard_Real ppU1, ppV1, ppU2, ppV2;
1426                           prevprevPoint.Parameters(ppU1, ppV1, ppU2, ppV2);
1427                           Standard_Real pU1, pV1, pU2, pV2;
1428                           previousPointSave.Parameters(pU1, pV1, pU2, pV2);
1429                           gp_Vec2d V1onS1(gp_Pnt2d(ppU1, ppV1), gp_Pnt2d(pU1, pV1));
1430                           gp_Vec2d V2onS1(gp_Pnt2d(pU1, pV1), gp_Pnt2d(u1, v1));
1431                           gp_Vec2d V1onS2(gp_Pnt2d(ppU2, ppV2), gp_Pnt2d(pU2, pV2));
1432                           gp_Vec2d V2onS2(gp_Pnt2d(pU2, pV2), gp_Pnt2d(u2, v2));
1433
1434                           const Standard_Real aDot1 = V1onS1 * V2onS1;
1435                           const Standard_Real aDot2 = V1onS2 * V2onS2;
1436
1437                           if ((aDot1 < 0.0) || (aDot2 < 0.0))
1438                           {
1439                             Arrive = Standard_True;
1440                             break;
1441                           }
1442                         }
1443                         /*
1444                         if(u1 <= UM1  && u2 <= UM2 && v1 <= VM1 &&
1445                         v2 <= VM2  && u1 >= Um1 && u2 >= Um2 &&
1446                         v1 >= Vm1  && v2 >= Vm2)  {
1447                         */                      
1448                         //xt
1449                         pl = previousPoint.Value();
1450
1451                         if(bTestFirstPoint)
1452                         {
1453                           if(pf.SquareDistance(pl) < aSQDistMax)
1454                           {
1455                             IncKey++;
1456
1457                             if(IncKey == 5000)
1458                               return;
1459                             else
1460                               continue;
1461                           }
1462                           else
1463                           {
1464                             bTestFirstPoint = Standard_False;
1465                           }
1466                         }
1467
1468                         //To avoid walking around the same point
1469                         //in the tangent zone near a border
1470
1471                         if (previoustg)
1472                         {
1473                           //There are three consecutive points:
1474                           //previousPointSave -> ParamPnt -> curPnt.
1475
1476                           Standard_Real prevU1, prevV1, prevU2, prevV2;
1477                           previousPointSave.Parameters(prevU1, prevV1, prevU2, prevV2);
1478                           gp_Pnt2d prevPntOnS1(prevU1, prevV1), prevPntOnS2(prevU2, prevV2);
1479                           gp_Pnt2d curPntOnS1(u1, v1), curPntOnS2(u2, v2);
1480                           gp_Vec2d PrevToParamOnS1(prevPntOnS1, ParamPntOnS1);
1481                           gp_Vec2d PrevToCurOnS1(prevPntOnS1, curPntOnS1);
1482                           gp_Vec2d PrevToParamOnS2(prevPntOnS2, ParamPntOnS2);
1483                           gp_Vec2d PrevToCurOnS2(prevPntOnS2, curPntOnS2);
1484                           Standard_Real MaxAngle = 3*M_PI/4;
1485                           Standard_Real anAngleS1 = 0.0, anAngleS2 = 0.0;
1486                           const Standard_Real aSQMParS1 = PrevToParamOnS1.SquareMagnitude();
1487                           const Standard_Real aSQMParS2 = PrevToParamOnS2.SquareMagnitude();
1488                           const Standard_Real aSQMCurS1 = PrevToCurOnS1.SquareMagnitude();
1489                           const Standard_Real aSQMCurS2 = PrevToCurOnS2.SquareMagnitude();
1490
1491                           if(aSQMCurS1 < gp::Resolution())
1492                           {
1493                             //We came back to the one of previos point.
1494                             //Therefore, we must break;
1495
1496                             anAngleS1 = M_PI;
1497                           }
1498                           else if(aSQMParS1 < gp::Resolution())
1499                           {
1500                             //We are walking along tangent zone.
1501                             //It should be continued.
1502                             anAngleS1 = 0.0;
1503                           }
1504                           else
1505                           {
1506                             anAngleS1 = Abs(PrevToParamOnS1.Angle(PrevToCurOnS1));
1507                           }
1508
1509                           if(aSQMCurS2 < gp::Resolution())
1510                           {
1511                             //We came back to the one of previos point.
1512                             //Therefore, we must break;
1513
1514                             anAngleS2 = M_PI;
1515                           }
1516                           else if(aSQMParS2 < gp::Resolution())
1517                           {
1518                             //We are walking along tangent zone.
1519                             //It should be continued;
1520                             anAngleS2 = 0.0;
1521                           }
1522                           else
1523                           {
1524                             anAngleS2 = Abs(PrevToParamOnS2.Angle(PrevToCurOnS2));
1525                           }
1526
1527                           if ((anAngleS1 > MaxAngle) && (anAngleS2 > MaxAngle))
1528                           {
1529                             Arrive = Standard_True;
1530                             break;
1531                           }
1532
1533                           {
1534                             //Check singularity.
1535                             //I.e. check if we are walking along direction, which does not
1536                             //result in comming to any point (i.e. derivative 
1537                             //3D-intersection curve along this direction is equal to 0).
1538                             //A sphere with direction {dU=1, dV=0} from point
1539                             //(U=0, V=M_PI/2) can be considered as example for
1540                             //this case (we cannot find another 3D-point if we go thus).
1541
1542                             //Direction chosen along 1st and 2nd surface correspondingly
1543                             const gp_Vec2d  aDirS1(prevPntOnS1, curPntOnS1),
1544                                             aDirS2(prevPntOnS2, curPntOnS2);
1545
1546                             gp_Pnt aPtemp;
1547                             gp_Vec aDuS1, aDvS1, aDuS2, aDvS2;
1548
1549                             myIntersectionOn2S.Function().AuxillarSurface1()->
1550                                   D1(curPntOnS1.X(), curPntOnS1.Y(), aPtemp, aDuS1, aDvS1);
1551                             myIntersectionOn2S.Function().AuxillarSurface2()->
1552                                   D1(curPntOnS2.X(), curPntOnS2.Y(), aPtemp, aDuS2, aDvS2);
1553
1554                             //Derivative WLine along (it is vector-function indeed)
1555                             //directions chosen
1556                             //(https://en.wikipedia.org/wiki/Directional_derivative#Variation_using_only_direction_of_vector).
1557                             //F1 - on the 1st surface, F2 - on the 2nd surface.
1558                             //x, y, z - coordinates of derivative vector.
1559                             const Standard_Real aF1x =  aDuS1.X()*aDirS1.X() + 
1560                                                         aDvS1.X()*aDirS1.Y();
1561                             const Standard_Real aF1y =  aDuS1.Y()*aDirS1.X() +
1562                                                         aDvS1.Y()*aDirS1.Y();
1563                             const Standard_Real aF1z =  aDuS1.Z()*aDirS1.X() +
1564                                                         aDvS1.Z()*aDirS1.Y();
1565                             const Standard_Real aF2x =  aDuS2.X()*aDirS2.X() +
1566                                                         aDvS2.X()*aDirS2.Y();
1567                             const Standard_Real aF2y =  aDuS2.Y()*aDirS2.X() +
1568                                                         aDvS2.Y()*aDirS2.Y();
1569                             const Standard_Real aF2z =  aDuS2.Z()*aDirS2.X() +
1570                                                         aDvS2.Z()*aDirS2.Y();
1571
1572                             const Standard_Real aF1 = aF1x*aF1x + aF1y*aF1y + aF1z*aF1z;
1573                             const Standard_Real aF2 = aF2x*aF2x + aF2y*aF2y + aF2z*aF2z;
1574
1575                             if((aF1 < gp::Resolution()) && (aF2 < gp::Resolution()))
1576                             {
1577                               //All derivative are equal to 0. Therefore, there is
1578                               //no point in going along direction chosen.
1579                               Arrive = Standard_True;
1580                               break;
1581                             }
1582                           }
1583                         }//if (previoustg) cond.
1584
1585                         ////////////////////////////////////////
1586                         AddAPoint(previousPoint);
1587                         RejectIndex++;
1588
1589                         if(RejectIndex >= RejectIndexMAX)
1590                         {
1591                           Arrive = Standard_True;
1592                           break;
1593                         }
1594
1595                         //
1596
1597                         LevelOfIterWithoutAppend=0;
1598                         Arrive = Standard_True;
1599                       }
1600                       else
1601                       {
1602                         // revert to the last correctly calculated point
1603                         previousPoint = previousPointSave;
1604                         previoustg    = previoustgSave;
1605                         previousd     = previousdSave;
1606                         previousd1    = previousd1Save;
1607                         previousd2    = previousd2Save;
1608                       }
1609                     }
1610
1611                     //
1612                     Standard_Boolean wasExtended = Standard_False;
1613
1614                     if(Arrive && myIntersectionOn2S.IsTangent() && bPrevNotTangent)
1615                     {
1616                       if(ExtendLineInCommonZone(SauvChoixIso, DejaReparti))
1617                       {
1618                         wasExtended = Standard_True;
1619                         Arrive = Standard_False;
1620                         ChoixIso = SauvChoixIso;
1621                       }
1622                     }
1623
1624                     RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
1625
1626                     if(Arrive && 
1627                       myIntersectionOn2S.IsDone() && !myIntersectionOn2S.IsEmpty() &&
1628                       myIntersectionOn2S.IsTangent() && bPrevNotTangent &&
1629                       !wasExtended)
1630                     {
1631                       if(ExtendLineInCommonZone(SauvChoixIso, DejaReparti))
1632                       {
1633                         wasExtended = Standard_True;
1634                         Arrive = Standard_False;
1635                         ChoixIso = SauvChoixIso;
1636                       }
1637                     }
1638                   }//else !TestArret() $
1639                 }//$$ end successful framing on border (!myIntersectionOn2S.IsEmpty())
1640                 else
1641                 {
1642                   //echec framing on border; division of step 
1643                   Arrive = Standard_False;
1644                   NoTestDeflection = Standard_True;
1645                   RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
1646                 }
1647               }//$$$ end framing on border (!close)
1648             }//004 fin TestArret return Arrive = True
1649           } // 006case IntWalk_ArretSurPoint:  end Processing Status = OK  or ArretSurPoint 
1650         } //007  switch(aStatus)
1651       } //008 end processing point  (TEST DEFLECTION)
1652     } //009 end processing line (else if myIntersectionOn2S.IsDone())
1653   }  //010 end if first departure point allows marching  while (!Arrive)
1654
1655   done = Standard_True;
1656 }
1657 // ===========================================================================================================
1658 // function: ExtendLineInCommonZone
1659 // purpose:  Extends already computed line inside tangent zone in the direction given by theChoixIso.
1660 //           Returns Standard_True if the line was extended through tangent zone and the last computed point 
1661 //           is outside the tangent zone (but it is not put into the line). Otherwise returns Standard_False.
1662 // ===========================================================================================================
1663 Standard_Boolean IntWalk_PWalking::ExtendLineInCommonZone(const IntImp_ConstIsoparametric theChoixIso,
1664                                                           const Standard_Boolean          theDirectionFlag) 
1665 {
1666   Standard_Boolean bOutOfTangentZone = Standard_False;
1667   Standard_Boolean bStop = !myIntersectionOn2S.IsTangent();
1668   Standard_Integer dIncKey = 1;
1669   TColStd_Array1OfReal Param(1,4);
1670   IntWalk_StatusDeflection aStatus = IntWalk_OK;
1671   Standard_Integer nbIterWithoutAppend = 0;
1672   Standard_Integer nbEqualPoints = 0;
1673   Standard_Integer parit = 0;
1674   Standard_Integer uvit = 0;
1675   IntSurf_SequenceOfPntOn2S aSeqOfNewPoint;
1676
1677   while (!bStop) {
1678     nbIterWithoutAppend++;
1679
1680     if((nbIterWithoutAppend > 20) || (nbEqualPoints > 20)) {
1681 #ifdef OCCT_DEBUG
1682       cout<<"Infinite loop detected. Stop iterations (IntWalk_PWalking_1.gxx)" << endl;
1683 #endif
1684       bStop = Standard_True;
1685       break;
1686     }
1687     Standard_Real f = 0.;
1688
1689     switch (theChoixIso)
1690     { 
1691     case IntImp_UIsoparametricOnCaro1: f = Abs(previousd1.X()); break;
1692     case IntImp_VIsoparametricOnCaro1: f = Abs(previousd1.Y()); break;
1693     case IntImp_UIsoparametricOnCaro2: f = Abs(previousd2.X()); break;
1694     case IntImp_VIsoparametricOnCaro2: f = Abs(previousd2.Y()); break;
1695     }
1696
1697     if(f<0.1) f=0.1;
1698
1699     previousPoint.Parameters(Param(1),Param(2),Param(3),Param(4));
1700
1701     Standard_Real dP1 = sensCheminement * pasuv[0] * previousd1.X() /f;
1702     Standard_Real dP2 = sensCheminement * pasuv[1] * previousd1.Y() /f;
1703     Standard_Real dP3 = sensCheminement * pasuv[2] * previousd2.X() /f; 
1704     Standard_Real dP4 = sensCheminement * pasuv[3] * previousd2.Y() /f;
1705
1706     if(theChoixIso == IntImp_UIsoparametricOnCaro1 && Abs(dP1) < 1.e-7) dP1 *= (5. * (Standard_Real)dIncKey);
1707     if(theChoixIso == IntImp_VIsoparametricOnCaro1 && Abs(dP2) < 1.e-7) dP2 *= (5. * (Standard_Real)dIncKey);
1708     if(theChoixIso == IntImp_UIsoparametricOnCaro2 && Abs(dP3) < 1.e-7) dP3 *= (5. * (Standard_Real)dIncKey);
1709     if(theChoixIso == IntImp_VIsoparametricOnCaro2 && Abs(dP4) < 1.e-7) dP4 *= (5. * (Standard_Real)dIncKey);
1710
1711     Param(1) += dP1;
1712     Param(2) += dP2;
1713     Param(3) += dP3; 
1714     Param(4) += dP4;
1715     Standard_Real SvParam[4];
1716     IntImp_ConstIsoparametric ChoixIso = theChoixIso;
1717
1718     for(parit = 0; parit < 4; parit++) {
1719       SvParam[parit] = Param(parit+1);
1720     }
1721     math_FunctionSetRoot  Rsnld(myIntersectionOn2S.Function());
1722     ChoixIso = myIntersectionOn2S.Perform(Param,Rsnld, theChoixIso);
1723
1724     if (!myIntersectionOn2S.IsDone()) {
1725       return bOutOfTangentZone;
1726     }
1727     else {
1728       if (myIntersectionOn2S.IsEmpty()) {
1729         return bOutOfTangentZone;
1730       }
1731
1732       aStatus = TestDeflection(ChoixIso);
1733
1734       if(aStatus == IntWalk_OK) {
1735
1736         for(uvit = 0; uvit < 4; uvit++) {
1737           if(pasuv[uvit] < pasInit[uvit]) {
1738             pasuv[uvit] = pasInit[uvit];
1739           }
1740         }
1741       }
1742
1743       switch(aStatus) {
1744       case  IntWalk_ArretSurPointPrecedent:
1745         {
1746           bStop = Standard_True;
1747           bOutOfTangentZone = !myIntersectionOn2S.IsTangent();
1748           break;
1749         }
1750       case IntWalk_PasTropGrand:
1751         {
1752           for(parit = 0; parit < 4; parit++) {
1753             Param(parit+1) = SvParam[parit];
1754           }
1755           Standard_Boolean bDecrease = Standard_False;
1756
1757           for(uvit = 0; uvit < 4; uvit++) {
1758             if(pasSav[uvit] < pasInit[uvit]) { 
1759               pasInit[uvit] -= (pasInit[uvit] - pasSav[uvit]) * 0.1;
1760               bDecrease = Standard_True;
1761             }
1762           }
1763
1764           if(bDecrease) nbIterWithoutAppend--;
1765           break;
1766         }
1767       case IntWalk_PointConfondu:
1768         {
1769           for(uvit = 0; uvit < 4; uvit++) {
1770             if(pasuv[uvit] < pasInit[uvit]) {
1771               pasuv[uvit] += (pasInit[uvit] - pasuv[uvit]) * 0.1;
1772             }
1773           }
1774           break;
1775         }
1776       case IntWalk_OK:
1777       case IntWalk_ArretSurPoint:
1778         {
1779           //
1780           bStop = TestArret(theDirectionFlag, Param, ChoixIso);
1781           //
1782
1783           //
1784           if(!bStop) {
1785             Standard_Real u11,v11,u12,v12; 
1786             myIntersectionOn2S.Point().Parameters(u11,v11,u12,v12); 
1787             Standard_Real u21,v21,u22,v22;
1788             previousPoint.Parameters(u21,v21,u22,v22); 
1789
1790             if(((fabs(u11-u21) < ResoU1) && (fabs(v11-v21) < ResoV1)) ||
1791               ((fabs(u12-u22) < ResoU2) && (fabs(v12-v22) < ResoV2))) {
1792                 nbEqualPoints++;
1793             }
1794             else {
1795               nbEqualPoints = 0;
1796             }
1797           }
1798           //
1799
1800           bStop = bStop || !myIntersectionOn2S.IsTangent();
1801           bOutOfTangentZone = !myIntersectionOn2S.IsTangent();
1802
1803           if(!bStop) {
1804             Standard_Boolean pointisvalid = Standard_False;
1805             Standard_Real u1,v1,u2,v2; 
1806             myIntersectionOn2S.Point().Parameters(u1,v1,u2,v2); 
1807
1808             if(u1 <= UM1  && u2 <= UM2 && v1 <= VM1 && 
1809               v2 <= VM2  && u1 >= Um1 && u2 >= Um2 &&
1810               v1 >= Vm1  && v2 >= Vm2) 
1811               pointisvalid = Standard_True;
1812
1813             if(pointisvalid) {
1814               previousPoint = myIntersectionOn2S.Point();
1815               previoustg = myIntersectionOn2S.IsTangent();
1816
1817               if(!previoustg) {
1818                 previousd  = myIntersectionOn2S.Direction();
1819                 previousd1 = myIntersectionOn2S.DirectionOnS1();
1820                 previousd2 = myIntersectionOn2S.DirectionOnS2();
1821               }
1822               Standard_Boolean bAddPoint = Standard_True;
1823
1824               if(line->NbPoints() >= 1) {
1825                 gp_Pnt pf = line->Value(1).Value();
1826                 gp_Pnt pl = previousPoint.Value(); 
1827
1828                 if(pf.Distance(pl) < Precision::Confusion()) { 
1829                   dIncKey++; 
1830                   if(dIncKey == 5000) return bOutOfTangentZone; 
1831                   else bAddPoint = Standard_False;
1832                 }
1833               }
1834
1835               if(bAddPoint) {
1836                 aSeqOfNewPoint.Append(previousPoint);
1837                 nbIterWithoutAppend = 0;
1838               }
1839             }
1840
1841             if (line->NbPoints() == 2) {
1842               for(uvit = 0; uvit < 4; uvit++) {
1843                 pasSav[uvit] = pasuv[uvit]; 
1844               }
1845             }
1846
1847             if ( !pointisvalid ) {
1848               // decrease step if out of bounds
1849               // otherwise the same calculations will be 
1850               // repeated several times
1851               if ( ( u1 > UM1 ) || ( u1 < Um1 ) )
1852                 pasuv[0] *= 0.5;
1853
1854               if ( ( v1 > VM1 ) || ( v1 < Vm1 ) ) 
1855                 pasuv[1] *= 0.5;
1856
1857               if ( ( u2 > UM2 ) || ( u2 < Um2 ) )
1858                 pasuv[2] *= 0.5;
1859
1860               if ( ( v2 > VM2 ) || ( v2 < Vm2 ) )
1861                 pasuv[3] *= 0.5;
1862             }
1863           } // end if(!bStop)
1864           else { //if(bStop)
1865             if(close && (line->NbPoints() >= 1)) { 
1866
1867               if(!bOutOfTangentZone) {
1868                 aSeqOfNewPoint.Append(line->Value(1)); // line end
1869               }
1870               nbIterWithoutAppend = 0;
1871             }
1872             else {
1873               ChoixIso = myIntersectionOn2S.Perform(Param, Rsnld, theChoixIso);
1874
1875               if(myIntersectionOn2S.IsEmpty()) { 
1876                 bStop = Standard_True;// !myIntersectionOn2S.IsTangent();
1877                 bOutOfTangentZone = Standard_False; // !myIntersectionOn2S.IsTangent();
1878               }
1879               else {
1880                 Standard_Boolean bAddPoint = Standard_True;
1881                 Standard_Boolean pointisvalid = Standard_False;
1882
1883                 previousPoint = myIntersectionOn2S.Point();
1884                 Standard_Real u1,v1,u2,v2; 
1885                 previousPoint.Parameters(u1,v1,u2,v2); 
1886
1887                 if(u1 <= UM1  && u2 <= UM2 && v1 <= VM1 && 
1888                   v2 <= VM2  && u1 >= Um1 && u2 >= Um2 &&
1889                   v1 >= Vm1  && v2 >= Vm2) 
1890                   pointisvalid = Standard_True;
1891
1892                 if(pointisvalid) {
1893
1894                   if(line->NbPoints() >= 1) {
1895                     gp_Pnt pf = line->Value(1).Value();
1896                     gp_Pnt pl = previousPoint.Value(); 
1897
1898                     if(pf.Distance(pl) < Precision::Confusion()) { 
1899                       dIncKey++; 
1900                       if(dIncKey == 5000) return bOutOfTangentZone; 
1901                       else bAddPoint = Standard_False;
1902                     }
1903                   }
1904
1905                   if(bAddPoint && !bOutOfTangentZone) {
1906                     aSeqOfNewPoint.Append(previousPoint);
1907                     nbIterWithoutAppend = 0;
1908                   }
1909                 }
1910               }
1911             }
1912           }
1913           break;
1914         }
1915       default:
1916         {
1917           break;
1918         }
1919       }
1920     }
1921   }
1922   Standard_Boolean bExtendLine = Standard_False;
1923   Standard_Real u1 = 0., v1 = 0., u2 = 0., v2 = 0.; 
1924
1925   Standard_Integer pit = 0;
1926
1927   for(pit = 0; !bExtendLine && (pit < 2); pit++) {
1928     if(pit == 0)
1929       previousPoint.Parameters(u1,v1,u2,v2); 
1930     else {
1931       if(aSeqOfNewPoint.Length() > 0)
1932         aSeqOfNewPoint.Value(aSeqOfNewPoint.Length()).Parameters(u1,v1,u2,v2); 
1933       else
1934         break;
1935     }
1936
1937     if(((u1 - Um1) < ResoU1) ||
1938       ((UM1 - u1) < ResoU1) ||
1939       ((u2 - Um2) < ResoU2) ||
1940       ((UM2 - u2) < ResoU2) ||
1941       ((v1 - Vm1) < ResoV1) ||
1942       ((VM1 - v1) < ResoV1) ||
1943       ((v2 - Vm2) < ResoV2) ||
1944       ((VM2 - v2) < ResoV2))
1945       bExtendLine = Standard_True;
1946   }
1947
1948   if(!bExtendLine) {
1949     //    if(aStatus == IntWalk_OK || aStatus == IntWalk_ArretSurPoint) {
1950     if(aStatus == IntWalk_OK) {
1951       bExtendLine = Standard_True;
1952
1953       if(aSeqOfNewPoint.Length() > 1) {
1954         TColStd_Array1OfReal FirstParams(0, 3), LastParams(0, 3), Resolutions(0, 3);
1955         Resolutions(0) = ResoU1; Resolutions(1) = ResoV1; Resolutions(2) = ResoU2; Resolutions(3) = ResoV2;
1956
1957         aSeqOfNewPoint(1).Parameters(FirstParams.ChangeValue(0), FirstParams.ChangeValue(1),
1958           FirstParams.ChangeValue(2), FirstParams.ChangeValue(3));
1959         aSeqOfNewPoint(aSeqOfNewPoint.Length()).Parameters(LastParams.ChangeValue(0), 
1960           LastParams.ChangeValue(1),
1961           LastParams.ChangeValue(2), 
1962           LastParams.ChangeValue(3)); 
1963         Standard_Integer indexofiso = 0;
1964
1965         if(theChoixIso == IntImp_UIsoparametricOnCaro1) indexofiso = 0;
1966         if(theChoixIso == IntImp_VIsoparametricOnCaro1) indexofiso = 1;
1967         if(theChoixIso == IntImp_UIsoparametricOnCaro2) indexofiso = 2;
1968         if(theChoixIso == IntImp_VIsoparametricOnCaro2) indexofiso = 3;
1969
1970         Standard_Integer afirstindex = (indexofiso < 2) ? 0 : 2;
1971         gp_Vec2d aTangentZoneDir(gp_Pnt2d(FirstParams.Value(afirstindex), FirstParams.Value(afirstindex + 1)),
1972           gp_Pnt2d(LastParams.Value(afirstindex), LastParams.Value(afirstindex + 1)));
1973
1974         gp_Dir2d anIsoDir(0, 1);
1975
1976         if((indexofiso == 1) || (indexofiso == 3))
1977           anIsoDir = gp_Dir2d(1, 0);
1978
1979         if(aTangentZoneDir.SquareMagnitude() > gp::Resolution()) {
1980           Standard_Real piquota = M_PI*0.25;
1981
1982           if(fabs(aTangentZoneDir.Angle(anIsoDir)) > piquota) {
1983             Standard_Integer ii = 1, nextii = 2;
1984             gp_Vec2d d1(0, 0);
1985             Standard_Real asqresol = gp::Resolution();
1986             asqresol *= asqresol;
1987
1988             do {
1989               aSeqOfNewPoint(ii).Parameters(FirstParams.ChangeValue(0), FirstParams.ChangeValue(1),
1990                 FirstParams.ChangeValue(2), FirstParams.ChangeValue(3));
1991               aSeqOfNewPoint(ii + 1).Parameters(LastParams.ChangeValue(0), LastParams.ChangeValue(1),
1992                 LastParams.ChangeValue(2), LastParams.ChangeValue(3));
1993               d1 = gp_Vec2d(gp_Pnt2d(FirstParams.Value(afirstindex),
1994                 FirstParams.Value(afirstindex + 1)),
1995                 gp_Pnt2d(LastParams.Value(afirstindex),
1996                 LastParams.Value(afirstindex + 1)));
1997               ii++;
1998             }
1999             while((d1.SquareMagnitude() < asqresol) &&
2000               (ii < aSeqOfNewPoint.Length()));
2001
2002             nextii = ii;
2003
2004             while(nextii < aSeqOfNewPoint.Length()) {
2005
2006               gp_Vec2d nextd1(0, 0);
2007               Standard_Integer jj = nextii;
2008
2009               do {
2010                 aSeqOfNewPoint(jj).Parameters(FirstParams.ChangeValue(0), FirstParams.ChangeValue(1),
2011                   FirstParams.ChangeValue(2), FirstParams.ChangeValue(3));
2012                 aSeqOfNewPoint(jj + 1).Parameters(LastParams.ChangeValue(0), LastParams.ChangeValue(1),
2013                   LastParams.ChangeValue(2), LastParams.ChangeValue(3));
2014                 nextd1 = gp_Vec2d(gp_Pnt2d(FirstParams.Value(afirstindex),
2015                   FirstParams.Value(afirstindex + 1)),
2016                   gp_Pnt2d(LastParams.Value(afirstindex),
2017                   LastParams.Value(afirstindex + 1)));
2018                 jj++;
2019
2020               }
2021               while((nextd1.SquareMagnitude() < asqresol) &&
2022                 (jj < aSeqOfNewPoint.Length()));
2023               nextii = jj;
2024
2025               if(fabs(d1.Angle(nextd1)) > piquota) {
2026                 bExtendLine = Standard_False;
2027                 break;
2028               }
2029               d1 = nextd1;
2030             }
2031           }
2032           // end if(fabs(aTangentZoneDir.Angle(anIsoDir)
2033         }
2034       }
2035     }
2036   }
2037
2038   if(!bExtendLine) {
2039     return Standard_False;
2040   }
2041   Standard_Integer i = 0;
2042
2043   for(i = 1; i <= aSeqOfNewPoint.Length(); i++) {
2044     AddAPoint(aSeqOfNewPoint.Value(i));
2045   }
2046
2047   return bOutOfTangentZone;
2048 }
2049
2050 //=======================================================================
2051 //function : DistanceMinimizeByGradient
2052 //purpose  : 
2053 //
2054 // ATTENTION!!!
2055 //  theInit should be initialized before function calling.
2056 //=======================================================================
2057 Standard_Boolean IntWalk_PWalking::
2058 DistanceMinimizeByGradient( const Handle(Adaptor3d_HSurface)& theASurf1,
2059                            const Handle(Adaptor3d_HSurface)& theASurf2,
2060                            TColStd_Array1OfReal& theInit,
2061                            const Standard_Real* theStep0)
2062 {
2063   const Standard_Integer aNbIterMAX = 60;
2064   const Standard_Real aTol = 1.0e-14;
2065   const Standard_Real aTolNul = 1.0 / Precision::Infinite();
2066
2067   // I.e. if theU1 = 0.0 then Epsilon(theU1) = DBL_MIN (~1.0e-308).
2068   // Work with this number is impossible: there is a dangerous to 
2069   // obtain Floating-point-overflow. Therefore, we limit this value.
2070   const Standard_Real aMinAddValU1 = Max(Epsilon(theInit(1)), aTolNul);
2071   const Standard_Real aMinAddValV1 = Max(Epsilon(theInit(2)), aTolNul);
2072   const Standard_Real aMinAddValU2 = Max(Epsilon(theInit(3)), aTolNul);
2073   const Standard_Real aMinAddValV2 = Max(Epsilon(theInit(4)), aTolNul);
2074
2075   Handle(Geom_Surface) aS1, aS2;
2076
2077   if (theASurf1->GetType() != GeomAbs_BezierSurface &&
2078       theASurf1->GetType() != GeomAbs_BSplineSurface)
2079       return Standard_True;
2080   if (theASurf2->GetType() != GeomAbs_BezierSurface &&
2081       theASurf2->GetType() != GeomAbs_BSplineSurface)
2082       return Standard_True;
2083
2084   Standard_Boolean aStatus = Standard_False;
2085
2086   gp_Pnt aP1, aP2;
2087   gp_Vec aD1u, aD1v, aD2U, aD2V;
2088
2089   theASurf1->D1(theInit(1), theInit(2), aP1, aD1u, aD1v);
2090   theASurf2->D1(theInit(3), theInit(4), aP2, aD2U, aD2V);
2091
2092   Standard_Real aSQDistPrev = aP1.SquareDistance(aP2);
2093
2094   gp_Vec aP12(aP1, aP2);
2095
2096   Standard_Real aGradFu(-aP12.Dot(aD1u));
2097   Standard_Real aGradFv(-aP12.Dot(aD1v));
2098   Standard_Real aGradFU( aP12.Dot(aD2U));
2099   Standard_Real aGradFV( aP12.Dot(aD2V));
2100
2101   Standard_Real aStepU1 = 1.0e-6, aStepV1 = 1.0e-6,
2102                 aStepU2 = 1.0e-6, aStepV2 = 1.0e-6;
2103     
2104   if (theStep0)
2105   {
2106     aStepU1 = theStep0[0];
2107     aStepV1 = theStep0[1];
2108     aStepU2 = theStep0[2];
2109     aStepV2 = theStep0[3];
2110   }
2111
2112   Standard_Boolean flRepeat = Standard_True;
2113   Standard_Integer aNbIter = aNbIterMAX;
2114
2115   while(flRepeat)
2116   {
2117     Standard_Real anAdd = aGradFu*aStepU1;
2118     const Standard_Real aPARu = theInit(1) - Sign(Max(Abs(anAdd), aMinAddValU1), anAdd);
2119     
2120     anAdd = aGradFv*aStepV1;
2121     const Standard_Real aPARv = theInit(2) - Sign(Max(Abs(anAdd), aMinAddValV1), anAdd);
2122
2123     anAdd = aGradFU*aStepU2;
2124     const Standard_Real aParU = theInit(3) - Sign(Max(Abs(anAdd), aMinAddValU2), anAdd);
2125
2126     anAdd = aGradFV*aStepV2;
2127     const Standard_Real aParV = theInit(4) - Sign(Max(Abs(anAdd), aMinAddValV2), anAdd);
2128
2129     gp_Pnt aPt1, aPt2;
2130
2131     theASurf1->D1(aPARu, aPARv, aPt1, aD1u, aD1v);
2132     theASurf2->D1(aParU, aParV, aPt2, aD2U, aD2V);
2133
2134     Standard_Real aSQDist = aPt1.SquareDistance(aPt2);
2135
2136     if(aSQDist < aSQDistPrev)
2137     {
2138       aSQDistPrev = aSQDist;
2139       theInit(1) = aPARu;
2140       theInit(2) = aPARv;
2141       theInit(3) = aParU;
2142       theInit(4) = aParV;
2143
2144       aStatus = aSQDistPrev < aTol;
2145       aStepU1 *= 1.2;
2146       aStepV1 *= 1.2;
2147       aStepU2 *= 1.2;
2148       aStepV2 *= 1.2;
2149     }
2150     else
2151     {
2152       if(--aNbIter < 0)
2153       {
2154         flRepeat = Standard_False;
2155       }
2156       else
2157       {
2158         theASurf1->D1(theInit(1), theInit(2), aPt1, aD1u, aD1v);
2159         theASurf2->D1(theInit(3), theInit(4), aPt2, aD2U, aD2V);
2160
2161         gp_Vec aPt12(aPt1, aPt2);
2162         aGradFu = -aPt12.Dot(aD1u);
2163         aGradFv = -aPt12.Dot(aD1v);
2164         aGradFU = aPt12.Dot(aD2U);
2165         aGradFV = aPt12.Dot(aD2V);
2166
2167         if (theStep0)
2168         {
2169           aStepU1 = theStep0[0];
2170           aStepV1 = theStep0[1];
2171           aStepU2 = theStep0[2];
2172           aStepV2 = theStep0[3];
2173         }
2174         else
2175         {
2176           aStepU1 = aStepV1 = aStepU2 = aStepV2 = 1.0e-6;
2177         }
2178       }
2179     }
2180   }
2181
2182   return aStatus;
2183 }
2184
2185 //=======================================================================
2186 //function : DistanceMinimizeByExtrema
2187 //purpose  : 
2188 //
2189 // ATTENTION!!!
2190 //  theP0, theU0 and theV0 parameters should be initialized
2191 //    before the function calling.
2192 //=======================================================================
2193 Standard_Boolean IntWalk_PWalking::
2194 DistanceMinimizeByExtrema(const Handle(Adaptor3d_HSurface)& theASurf, 
2195                           const gp_Pnt& theP0,
2196                           Standard_Real& theU0,
2197                           Standard_Real& theV0,
2198                           const Standard_Real* theStep0)
2199 {
2200   const Standard_Real aTol = 1.0e-14;
2201   gp_Pnt aPS;
2202   gp_Vec aD1Su, aD1Sv, aD2Su, aD2Sv, aD2SuvTemp;
2203   Standard_Real aSQDistPrev = RealLast();
2204   Standard_Real aU = theU0, aV = theV0;
2205
2206   const Standard_Real aStep0[2] = { theStep0 ? theStep0[0] : 1.0,
2207                                     theStep0 ? theStep0[1] : 1.0 };
2208
2209   Standard_Integer aNbIter = 10;
2210   do
2211   {
2212     theASurf->D2(aU, aV, aPS, aD1Su, aD1Sv, aD2Su, aD2Sv, aD2SuvTemp);
2213
2214     gp_Vec aVec(theP0, aPS);
2215
2216     Standard_Real aSQDist = aVec.SquareMagnitude();
2217
2218     if(aSQDist >= aSQDistPrev)
2219       break;
2220
2221     aSQDistPrev = aSQDist;
2222     theU0 = aU;
2223     theV0 = aV;
2224     aNbIter--;
2225
2226     if(aSQDistPrev < aTol)
2227       break;
2228
2229     //Functions
2230     const Standard_Real aF1 = aD1Su.Dot(aVec), aF2 = aD1Sv.Dot(aVec);
2231
2232     //Derivatives
2233     const Standard_Real aDf1u = aD2Su.Dot(aVec) + aD1Su.Dot(aD1Su),
2234       aDf1v = aD2Su.Dot(aD1Sv),
2235       aDf2u = aDf1v,
2236       aDf2v = aD2Sv.Dot(aVec) + aD1Sv.Dot(aD1Sv);
2237
2238     const Standard_Real aDet = aDf1u*aDf2v - aDf1v*aDf2u;
2239     aU -= aStep0[0]*(aDf2v*aF1 - aDf1v*aF2) / aDet;
2240     aV += aStep0[1]*(aDf2u*aF1 - aDf1u*aF2) / aDet;
2241   }
2242   while(aNbIter > 0);
2243
2244   return (aSQDistPrev < aTol);
2245 }
2246
2247 //=======================================================================
2248 //function : HandleSingleSingularPoint
2249 //purpose  : 
2250 //=======================================================================
2251 Standard_Boolean IntWalk_PWalking::HandleSingleSingularPoint(const Handle(Adaptor3d_HSurface)& theASurf1,
2252                                                              const Handle(Adaptor3d_HSurface)& theASurf2,
2253                                                              const Standard_Real the3DTol,
2254                                                              TColStd_Array1OfReal &thePnt)
2255 {
2256   // u1, v1, u2, v2 order is used.
2257   Standard_Real aLowBorder[4] = {theASurf1->FirstUParameter(),
2258                                  theASurf1->FirstVParameter(),
2259                                  theASurf2->FirstUParameter(),
2260                                  theASurf2->FirstVParameter()};
2261   Standard_Real aUppBorder[4] = {theASurf1->LastUParameter(),
2262                                  theASurf1->LastVParameter(),
2263                                  theASurf2->LastUParameter(),
2264                                  theASurf2->LastVParameter()};
2265   IntImp_ConstIsoparametric aLockedDir[4] = {IntImp_UIsoparametricOnCaro1,
2266                                              IntImp_VIsoparametricOnCaro1,
2267                                              IntImp_UIsoparametricOnCaro2,
2268                                              IntImp_VIsoparametricOnCaro2};
2269
2270   // Create new intersector with new tolerance.
2271   IntWalk_TheInt2S anInt(theASurf1, theASurf2, the3DTol);
2272   math_FunctionSetRoot aRsnld(anInt.Function());
2273
2274   for (Standard_Integer i = 1; i <= 4; ++i)
2275   {
2276     if ( Abs(thePnt(i) - aLowBorder[i - 1]) < Precision::PConfusion() ||
2277          Abs(thePnt(i) - aUppBorder[i - 1]) < Precision::PConfusion())
2278     {
2279
2280       anInt.Perform(thePnt,aRsnld, aLockedDir[i - 1]);
2281
2282       if (!anInt.IsDone())
2283         continue;
2284
2285       if (anInt.IsEmpty())
2286         continue;
2287
2288       anInt.Point().Parameters(thePnt(1), thePnt(2), thePnt(3), thePnt(4));
2289
2290       Standard_Boolean isInDomain = Standard_True;
2291       for (Standard_Integer j = 1; isInDomain && (j <= 4); ++j)
2292       {
2293         if ((thePnt(j) - aLowBorder[j - 1] + Precision::PConfusion())*
2294             (thePnt(j) - aUppBorder[j - 1] - Precision::PConfusion()) > 0.0)
2295         {
2296           isInDomain = Standard_False;
2297         }
2298       }
2299
2300       if (isInDomain)
2301         return Standard_True;
2302     }
2303   }
2304
2305   return Standard_False;
2306 }
2307
2308 //=======================================================================
2309 //function : SeekPointOnBoundary
2310 //purpose  : 
2311 //=======================================================================
2312 Standard_Boolean IntWalk_PWalking::
2313         SeekPointOnBoundary(const Handle(Adaptor3d_HSurface)& theASurf1,
2314                             const Handle(Adaptor3d_HSurface)& theASurf2,
2315                             const Standard_Real theU1,
2316                             const Standard_Real theV1,
2317                             const Standard_Real theU2,
2318                             const Standard_Real theV2,
2319                             const Standard_Boolean isTheFirst)
2320 {
2321   Standard_Boolean isOK = Standard_False;
2322
2323   // u1, v1, u2, v2 order is used.
2324   const Standard_Real aLowBorder[4] = {theASurf1->FirstUParameter(),
2325                                        theASurf1->FirstVParameter(),
2326                                        theASurf2->FirstUParameter(),
2327                                        theASurf2->FirstVParameter()};
2328   const Standard_Real aUppBorder[4] = {theASurf1->LastUParameter(),
2329                                        theASurf1->LastVParameter(),
2330                                        theASurf2->LastUParameter(),
2331                                        theASurf2->LastVParameter()};
2332
2333   // Tune solution tolerance according with object size.
2334   const Standard_Real aRes1 = Max(Precision::PConfusion() / theASurf1->UResolution(1.0),
2335                                   Precision::PConfusion() / theASurf1->VResolution(1.0));
2336   const Standard_Real aRes2 = Max(Precision::PConfusion() / theASurf2->UResolution(1.0),
2337                                   Precision::PConfusion() / theASurf2->VResolution(1.0));
2338   const Standard_Real a3DTol = Max(aRes1, aRes2);
2339   const Standard_Real aTol = Max(Precision::Confusion(), a3DTol);
2340
2341   // u1, v1, u2, v2 order is used.
2342   TColStd_Array1OfReal aPnt(1,4);
2343   aPnt(1) = theU1; aPnt(2) = theV1; aPnt(3) = theU2; aPnt(4) = theV2;
2344   TColStd_Array1OfReal aSingularPnt(aPnt);
2345
2346   Standard_Integer aNbIter = 20;
2347   Standard_Boolean aStatus = Standard_False;
2348   do
2349   {
2350     aNbIter--;
2351     aStatus = DistanceMinimizeByGradient(theASurf1, theASurf2, aPnt);
2352     if (aStatus && !AdjustToDomain(4, &aPnt(1), &aLowBorder[0], &aUppBorder[0]))
2353       break;
2354
2355     aStatus = DistanceMinimizeByExtrema(theASurf1, theASurf2->Value(aPnt(3), aPnt(4)),
2356                                         aPnt(1), aPnt(2));
2357     if (aStatus && !AdjustToDomain(2, &aPnt(1), &aLowBorder[0], &aUppBorder[0]))
2358       break;
2359
2360     aStatus = DistanceMinimizeByExtrema(theASurf2, theASurf1->Value(aPnt(1), aPnt(2)),
2361                                         aPnt(3), aPnt(4));
2362     if (aStatus && !AdjustToDomain(2, &aPnt(3), &aLowBorder[2], &aUppBorder[2]))
2363       break;
2364   }
2365   while(!aStatus && (aNbIter > 0));
2366
2367   // Handle singular points.
2368   Standard_Boolean aSingularStatus = HandleSingleSingularPoint(theASurf1,
2369                                                                theASurf2,
2370                                                                aTol,
2371                                                                aSingularPnt);
2372   if (aSingularStatus)
2373     aPnt = aSingularPnt;
2374
2375   if (!aStatus && !aSingularStatus)
2376   {
2377     return isOK;
2378   }
2379
2380   gp_Pnt aP1 = theASurf1->Value(aPnt(1), aPnt(2));
2381   gp_Pnt aP2 = theASurf2->Value(aPnt(3), aPnt(4));
2382   const gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));
2383
2384   const Standard_Real aSQDist = aPInt.SquareDistance(aP1);
2385   if (aSQDist > aTol * aTol)
2386   {
2387     return isOK;
2388   }
2389
2390   //Found point is true intersection point
2391   IntSurf_PntOn2S anIP;
2392   anIP.SetValue(aPInt, aPnt(1), aPnt(2), aPnt(3), aPnt(4));
2393
2394   //The main idea of checks below is to define if insertion of
2395   //addition point (on the boundary) does not lead to invalid
2396   //intersection curve (e.g. having a loop).
2397   //
2398   //Loops are detected with rotation angle of the Walking-line (WL).
2399   //If there is hairpin bend then insertion is forbidden.
2400
2401   //There are at least two possible problems:
2402   //  1. There are some cases when two neighbor points of the WL
2403   //      are almost coincident (the distance between them is less
2404   //      than Precision::Confusion). It is impossible to define
2405   //      rotation angle in these cases. Therefore, points with
2406   //      "good" distances should be selected.
2407
2408   //  2. Intersection point on the surface boundary has highest
2409   //      priority in compare with other "middle" points. Therefore,
2410   //      if insertion of new point will result in a bend then some
2411   //      "middle" points should be deleted in order to provide
2412   //      correct insertion.
2413
2414   //Problem test cases:
2415   //  test bugs modalg_5 bug24585_1
2416   //  test boolean bcut_complex G7
2417   //  test bugs moddata_2 bug469
2418
2419   if (isTheFirst)
2420   {
2421     while (line->NbPoints() > 1)
2422     {
2423       const Standard_Integer aNbPnts = line->NbPoints();
2424
2425       Standard_Integer aPInd = 1;
2426       for (; aPInd <= aNbPnts; aPInd++)
2427       {
2428         aP1.SetXYZ(line->Value(aPInd).Value().XYZ());
2429         if (aP1.SquareDistance(aPInt) > Precision::SquareConfusion())
2430         {
2431           break;
2432         }
2433         else if (aPInd == 1)
2434         {
2435           // After insertion, we will obtain
2436           // two coincident points in the line.
2437           // Therefore, insertion is forbidden.
2438           return isOK;
2439         }
2440       }
2441
2442       for (++aPInd; aPInd <= aNbPnts; aPInd++)
2443       {
2444         aP2.SetXYZ(line->Value(aPInd).Value().XYZ());
2445         if (aP1.SquareDistance(aP2) > Precision::SquareConfusion())
2446           break;
2447       }
2448
2449       if (aPInd > aNbPnts)
2450       {
2451         return isOK;
2452       }
2453
2454       const gp_XYZ aDir01(aP1.XYZ() - aPInt.XYZ());
2455       const gp_XYZ aDir12(aP2.XYZ() - aP1.XYZ());
2456
2457       if (aDir01.Dot(aDir12) > 0.0)
2458       {
2459         break;
2460       }
2461
2462       RemoveAPoint(1);
2463     }
2464
2465     line->InsertBefore(1, anIP);
2466     isOK = Standard_True;
2467   }
2468   else
2469   {
2470     while (line->NbPoints() > 1)
2471     {
2472       const Standard_Integer aNbPnts = line->NbPoints();
2473
2474       gp_Pnt aPPrev, aPCurr;
2475       Standard_Integer aPInd = aNbPnts;
2476       for (; aPInd > 0; aPInd--)
2477       {
2478         aPCurr.SetXYZ(line->Value(aPInd).Value().XYZ());
2479         if (aPCurr.SquareDistance(aPInt) > Precision::SquareConfusion())
2480         {
2481           break;
2482         }
2483         else if (aPInd == aNbPnts)
2484         {
2485           // After insertion, we will obtain
2486           // two coincident points in the line.
2487           // Therefore, insertion is forbidden.
2488           return isOK;
2489         }
2490       }
2491
2492       for (--aPInd; aPInd > 0; aPInd--)
2493       {
2494         aPPrev.SetXYZ(line->Value(aPInd).Value().XYZ());
2495         if (aPCurr.SquareDistance(aPPrev) > Precision::SquareConfusion())
2496           break;
2497       }
2498
2499       if (aPInd < 1)
2500       {
2501         return isOK;
2502       }
2503
2504       const gp_XYZ aDirPC(aPCurr.XYZ() - aPPrev.XYZ());
2505       const gp_XYZ aDirCN(aPInt.XYZ() - aPCurr.XYZ());
2506
2507       if (aDirPC.Dot(aDirCN) > 0.0)
2508       {
2509         break;
2510       }
2511
2512       RemoveAPoint(aNbPnts);
2513     }
2514
2515     line->Add(anIP);
2516     isOK = Standard_True;
2517   }
2518
2519   return isOK;
2520 }
2521
2522 //=======================================================================
2523 //function : PutToBoundary
2524 //purpose  : 
2525 //=======================================================================
2526 Standard_Boolean IntWalk_PWalking::
2527 PutToBoundary(const Handle(Adaptor3d_HSurface)& theASurf1,
2528               const Handle(Adaptor3d_HSurface)& theASurf2)
2529 {
2530   const Standard_Real aTolMin = Precision::Confusion();
2531
2532   Standard_Boolean hasBeenAdded = Standard_False;
2533
2534   const Standard_Real aU1bFirst = theASurf1->FirstUParameter();
2535   const Standard_Real aU1bLast = theASurf1->LastUParameter();
2536   const Standard_Real aU2bFirst = theASurf2->FirstUParameter();
2537   const Standard_Real aU2bLast = theASurf2->LastUParameter();
2538   const Standard_Real aV1bFirst = theASurf1->FirstVParameter();
2539   const Standard_Real aV1bLast = theASurf1->LastVParameter();
2540   const Standard_Real aV2bFirst = theASurf2->FirstVParameter();
2541   const Standard_Real aV2bLast = theASurf2->LastVParameter();
2542
2543   Standard_Real aTol = 1.0;
2544   aTol = Min(aTol, aU1bLast - aU1bFirst);
2545   aTol = Min(aTol, aU2bLast - aU2bFirst);
2546   aTol = Min(aTol, aV1bLast - aV1bFirst);
2547   aTol = Min(aTol, aV2bLast - aV2bFirst)*1.0e-3;
2548
2549   if(aTol <= 2.0*aTolMin)
2550     return hasBeenAdded;
2551
2552   Standard_Boolean isNeedAdding = Standard_False;
2553   Standard_Boolean isU1parallel = Standard_False, isV1parallel = Standard_False;
2554   Standard_Boolean isU2parallel = Standard_False, isV2parallel = Standard_False;
2555   IsParallel(line, Standard_True, aTol, isU1parallel, isV1parallel);
2556   IsParallel(line, Standard_False, aTol, isU2parallel, isV2parallel);
2557
2558   Standard_Real u1, v1, u2, v2;
2559   line->Value(1).Parameters(u1, v1, u2, v2);
2560   Standard_Real aDelta = 0.0;
2561
2562   if(!isV1parallel)
2563   {
2564     aDelta = u1 - aU1bFirst;
2565     if((aTolMin < aDelta) && (aDelta < aTol))
2566     {
2567       u1 = aU1bFirst;
2568       isNeedAdding = Standard_True;
2569     }
2570     else
2571     {
2572       aDelta = aU1bLast - u1;
2573       if((aTolMin < aDelta) && (aDelta < aTol))
2574       {
2575         u1 = aU1bLast;
2576         isNeedAdding = Standard_True;
2577       }
2578     }
2579   }
2580
2581   if(!isV2parallel)
2582   {
2583     aDelta = u2 - aU2bFirst;
2584     if((aTolMin < aDelta) && (aDelta < aTol))
2585     {
2586       u2 = aU2bFirst;
2587       isNeedAdding = Standard_True;
2588     }
2589     else
2590     {
2591       aDelta = aU2bLast - u2;
2592       if((aTolMin < aDelta) && (aDelta < aTol))
2593       {
2594         u2 = aU2bLast;
2595         isNeedAdding = Standard_True;
2596       }
2597     }
2598   }
2599
2600   if(!isU1parallel)
2601   {
2602     aDelta = v1 - aV1bFirst;
2603     if((aTolMin < aDelta) && (aDelta < aTol))
2604     {
2605       v1 = aV1bFirst;
2606       isNeedAdding = Standard_True;
2607     }
2608     else
2609     {
2610       aDelta = aV1bLast - v1;
2611       if((aTolMin < aDelta) && (aDelta < aTol))
2612       {
2613         v1 = aV1bLast;
2614         isNeedAdding = Standard_True;
2615       }
2616     }
2617   }
2618
2619   if(!isU2parallel)
2620   {
2621     aDelta = v2 - aV2bFirst;
2622     if((aTolMin < aDelta) && (aDelta < aTol))
2623     {
2624       v2 = aV2bFirst;
2625       isNeedAdding = Standard_True;
2626     }
2627     else
2628     {
2629       aDelta = aV2bLast - v2;
2630       if((aTolMin < aDelta) && (aDelta < aTol))
2631       {
2632         v2 = aV2bLast;
2633         isNeedAdding = Standard_True;
2634       }
2635     }
2636   }
2637
2638   if(isNeedAdding)
2639   {
2640     hasBeenAdded = 
2641       SeekPointOnBoundary(theASurf1, theASurf2, u1, 
2642       v1, u2, v2, Standard_True);
2643   }
2644
2645   const Standard_Integer aNbPnts = line->NbPoints();
2646   isNeedAdding = Standard_False;
2647   line->Value(aNbPnts).Parameters(u1, v1, u2, v2);
2648
2649   if(!isV1parallel)
2650   {
2651     aDelta = u1 - aU1bFirst;
2652     if((aTolMin < aDelta) && (aDelta < aTol))
2653     {
2654       u1 = aU1bFirst;
2655       isNeedAdding = Standard_True;
2656     }
2657     else
2658     {
2659       aDelta = aU1bLast - u1;
2660       if((aTolMin < aDelta) && (aDelta < aTol))
2661       {
2662         u1 = aU1bLast;
2663         isNeedAdding = Standard_True;
2664       }
2665     }
2666   }
2667
2668   if(!isV2parallel)
2669   {
2670     aDelta = u2 - aU2bFirst;
2671     if((aTolMin < aDelta) && (aDelta < aTol))
2672     {
2673       u2 = aU2bFirst;
2674       isNeedAdding = Standard_True;
2675     }
2676     else
2677     {
2678       aDelta = aU2bLast - u2;
2679       if((aTolMin < aDelta) && (aDelta < aTol))
2680       {
2681         u2 = aU2bLast;
2682         isNeedAdding = Standard_True;
2683       }
2684     }
2685   }
2686
2687   if(!isU1parallel)
2688   {
2689     aDelta = v1 - aV1bFirst;
2690     if((aTolMin < aDelta) && (aDelta < aTol))
2691     {
2692       v1 = aV1bFirst;
2693       isNeedAdding = Standard_True;
2694     }
2695     else
2696     {
2697       aDelta = aV1bLast - v1;
2698       if((aTolMin < aDelta) && (aDelta < aTol))
2699       {
2700         v1 = aV1bLast;
2701         isNeedAdding = Standard_True;
2702       }
2703     }
2704   }
2705
2706   if(!isU2parallel)
2707   {
2708     aDelta = v2 - aV2bFirst;
2709     if((aTolMin < aDelta) && (aDelta < aTol))
2710     {
2711       v2 = aV2bFirst;
2712       isNeedAdding = Standard_True;
2713     }
2714     else
2715     {
2716       aDelta = aV2bLast - v2;
2717       if((aTolMin < aDelta) && (aDelta < aTol))
2718       {
2719         v2 = aV2bLast;
2720         isNeedAdding = Standard_True;
2721       }
2722     }
2723   }
2724
2725   if(isNeedAdding)
2726   {
2727     hasBeenAdded = 
2728       SeekPointOnBoundary(theASurf1, theASurf2, u1, 
2729       v1, u2, v2, Standard_False);
2730   }
2731
2732   return hasBeenAdded;
2733 }
2734
2735 //=======================================================================
2736 //function : SeekAdditionalPoints
2737 //purpose  : 
2738 //=======================================================================
2739 Standard_Boolean IntWalk_PWalking::
2740 SeekAdditionalPoints( const Handle(Adaptor3d_HSurface)& theASurf1,
2741                      const Handle(Adaptor3d_HSurface)& theASurf2,
2742                      const Standard_Integer theMinNbPoints)
2743 {
2744   const Standard_Real aTol = 1.0e-14;
2745   Standard_Integer aNbPoints = line->NbPoints();
2746   if(aNbPoints > theMinNbPoints)
2747     return Standard_True;
2748
2749   const Standard_Real aU1bFirst = theASurf1->FirstUParameter();
2750   const Standard_Real aU1bLast = theASurf1->LastUParameter();
2751   const Standard_Real aU2bFirst = theASurf2->FirstUParameter();
2752   const Standard_Real aU2bLast = theASurf2->LastUParameter();
2753   const Standard_Real aV1bFirst = theASurf1->FirstVParameter();
2754   const Standard_Real aV1bLast = theASurf1->LastVParameter();
2755   const Standard_Real aV2bFirst = theASurf2->FirstVParameter();
2756   const Standard_Real aV2bLast = theASurf2->LastVParameter();
2757
2758
2759   Standard_Boolean isPrecise = Standard_False;
2760
2761   TColStd_Array1OfReal aPnt(1, 4);
2762   aPnt.Init(0.0);
2763
2764   Standard_Integer aNbPointsPrev = 0;
2765   while(aNbPoints < theMinNbPoints && (aNbPoints != aNbPointsPrev))
2766   {
2767     aNbPointsPrev = aNbPoints;
2768     for(Standard_Integer fp = 1, lp = 2; fp < aNbPoints; fp = lp + 1)
2769     {
2770       Standard_Real U1f, V1f, U2f, V2f; //first point in 1st and 2nd surafaces
2771       Standard_Real U1l, V1l, U2l, V2l; //last  point in 1st and 2nd surafaces
2772
2773       lp = fp+1;
2774       line->Value(fp).Parameters(U1f, V1f, U2f, V2f);
2775       line->Value(lp).Parameters(U1l, V1l, U2l, V2l);
2776
2777       aPnt(1) = 0.5*(U1f + U1l);
2778       if(aPnt(1) < aU1bFirst)
2779         aPnt(1) = aU1bFirst;
2780       if(aPnt(1) > aU1bLast)
2781         aPnt(1) = aU1bLast;
2782
2783       aPnt(2) = 0.5*(V1f+V1l);
2784       if(aPnt(2) < aV1bFirst)
2785         aPnt(2) = aV1bFirst;
2786       if(aPnt(2) > aV1bLast)
2787         aPnt(2) = aV1bLast;
2788
2789       aPnt(3) = 0.5*(U2f+U2l);
2790       if(aPnt(3) < aU2bFirst)
2791         aPnt(3) = aU2bFirst;
2792       if(aPnt(3) > aU2bLast)
2793         aPnt(3) = aU2bLast;
2794
2795       aPnt(4) = 0.5*(V2f+V2l);
2796       if(aPnt(4) < aV2bFirst)
2797         aPnt(4) = aV2bFirst;
2798       if(aPnt(4) > aV2bLast)
2799         aPnt(4) = aV2bLast;
2800
2801       Standard_Boolean aStatus = Standard_False;
2802       Standard_Integer aNbIter = 5;
2803       do
2804       {
2805         aStatus = DistanceMinimizeByGradient(theASurf1, theASurf2, aPnt);
2806         if(aStatus)
2807         {
2808           break;
2809         }
2810
2811         aStatus = DistanceMinimizeByExtrema(theASurf1, theASurf2->Value(aPnt(3), aPnt(4)), aPnt(1), aPnt(2));
2812         if(aStatus)
2813         {
2814           break;
2815         }
2816
2817         aStatus = DistanceMinimizeByExtrema(theASurf2, theASurf1->Value(aPnt(1), aPnt(2)), aPnt(3), aPnt(4));
2818         if(aStatus)
2819         {
2820           break;
2821         }
2822       }
2823       while(!aStatus && (--aNbIter > 0));
2824
2825       if(aStatus)
2826       {
2827         gp_Pnt  aP1 = theASurf1->Value(aPnt(1), aPnt(2)),
2828           aP2 = theASurf2->Value(aPnt(3), aPnt(4));
2829         gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));
2830
2831         const Standard_Real aSQDist1 = aPInt.SquareDistance(aP1),
2832           aSQDist2 = aPInt.SquareDistance(aP2);
2833
2834         if((aSQDist1 < aTol) && (aSQDist2 < aTol))
2835         {
2836           IntSurf_PntOn2S anIP;
2837           anIP.SetValue(aPInt, aPnt(1), aPnt(2), aPnt(3), aPnt(4));
2838           line->InsertBefore(lp, anIP);
2839
2840           isPrecise = Standard_True;
2841
2842           if(++aNbPoints >= theMinNbPoints)
2843             break;
2844         }
2845         else
2846         {
2847           lp--;
2848         }
2849       }
2850     }
2851   }
2852
2853   return isPrecise;
2854 }
2855
2856 void IntWalk_PWalking::
2857 RepartirOuDiviser(Standard_Boolean& DejaReparti,
2858                   IntImp_ConstIsoparametric& ChoixIso,
2859                   Standard_Boolean& Arrive) 
2860
2861                   // at the neighborhood of a point, there is a fail of marching 
2862                   // it is required to divide the steps to try to continue
2863                   // if the step is too small if we are on border
2864                   // restart in another direction if it was not done, otherwise stop
2865
2866 {
2867   //  Standard_Integer i;
2868   if (Arrive) {    //restart in the other direction
2869     if (!DejaReparti ) {
2870       Arrive = Standard_False;
2871       DejaReparti = Standard_True;
2872       previousPoint = line->Value(1);
2873       previoustg = Standard_False;
2874       previousd1 = firstd1;
2875       previousd2 = firstd2;
2876       previousd = tgdir;
2877       myTangentIdx = line->NbPoints();
2878       tgdir.Reverse();
2879       line->Reverse();
2880
2881       //-- printf("\nIntWalk_PWalking_2.gxx Reverse %3d\n",indextg);
2882       sensCheminement = -1;
2883       tgfirst      = tglast;
2884       tglast       = Standard_False;
2885       ChoixIso     = choixIsoSav;
2886 #if 0
2887       pasuv[0]=pasSav[0];
2888       pasuv[1]=pasSav[1];
2889       pasuv[2]=pasSav[2];
2890       pasuv[3]=pasSav[3];
2891 #else 
2892       Standard_Real u1,v1,u2,v2;
2893       Standard_Real U1,V1,U2,V2;
2894       Standard_Integer nn=line->NbPoints();
2895       if(nn>2) { 
2896         line->Value(nn).Parameters(u1,v1,u2,v2);
2897         line->Value(nn-1).Parameters(U1,V1,U2,V2);
2898         pasuv[0]=Abs(u1-U1);
2899         pasuv[1]=Abs(v1-V1);
2900         pasuv[2]=Abs(u2-U2);
2901         pasuv[3]=Abs(v2-V2);
2902       }
2903 #endif
2904
2905     }
2906   }  
2907   else  {
2908     if (    pasuv[0]*0.5 < ResoU1
2909       &&  pasuv[1]*0.5 < ResoV1
2910       &&  pasuv[2]*0.5 < ResoU2
2911       &&  pasuv[3]*0.5 < ResoV2
2912       ) {
2913         if (!previoustg) {
2914           tglast = Standard_True;      // IS IT ENOUGH ????
2915         }
2916
2917         if (!DejaReparti)
2918         {  //restart in the other direction
2919           DejaReparti = Standard_True;
2920           previousPoint = line->Value(1);
2921           previoustg = Standard_False;
2922           previousd1 = firstd1;
2923           previousd2 = firstd2;
2924           previousd = tgdir;
2925           myTangentIdx = line->NbPoints();
2926           tgdir.Reverse();
2927           line->Reverse();
2928
2929           //-- printf("\nIntWalk_PWalking_2.gxx Reverse %3d\n",indextg);
2930
2931           sensCheminement   = -1;
2932           tgfirst           = tglast;
2933           tglast            = Standard_False;
2934           ChoixIso          = choixIsoSav;
2935
2936 #if 0 
2937           pasuv[0]=pasSav[0];
2938           pasuv[1]=pasSav[1];
2939           pasuv[2]=pasSav[2];
2940           pasuv[3]=pasSav[3];
2941 #else 
2942           Standard_Real u1,v1,u2,v2;
2943           Standard_Real U1,V1,U2,V2;
2944           Standard_Integer nn=line->NbPoints();
2945           if(nn>2) { 
2946             line->Value(nn).Parameters(u1,v1,u2,v2);
2947             line->Value(nn-1).Parameters(U1,V1,U2,V2);
2948             pasuv[0]=Abs(u1-U1);
2949             pasuv[1]=Abs(v1-V1);
2950             pasuv[2]=Abs(u2-U2);
2951             pasuv[3]=Abs(v2-V2);
2952           }
2953 #endif
2954         }
2955         else Arrive = Standard_True;
2956     }
2957     else {
2958       pasuv[0]*=0.5;
2959       pasuv[1]*=0.5;
2960       pasuv[2]*=0.5;
2961       pasuv[3]*=0.5; 
2962     }
2963   }
2964 }
2965
2966 namespace {
2967   //OCC431(apo): modified ->
2968   static const Standard_Real CosRef2D =  Cos(M_PI/9.0),  AngRef2D = M_PI/2.0; 
2969
2970   static const Standard_Real d = 7.0;
2971 }
2972
2973 IntWalk_StatusDeflection  IntWalk_PWalking::TestDeflection(const IntImp_ConstIsoparametric choixIso)
2974
2975 // test if vector is observed by calculating an increase of vector 
2976 //     or the previous point and its tangent, the new calculated point and its  
2977 //     tangent; it is possible to find a cube passing by the 2 points and having as a 
2978 //     derivative the tangents of the intersection
2979 //     calculate the point with parameter 0.5 on cube=p1 
2980 //     calculate the medium point of 2 points of intersection=p2
2981 //   if arrow/2<=||p1p2||<= arrow consider that the vector is observed
2982 //   otherwise adjust the step depending on the ratio ||p1p2||/vector
2983 //   and the previous step 
2984 // test if in  2 tangent planes of surfaces there is no too great angle2d 
2985 // grand : if yes divide the step
2986 // test if there is no change of side
2987 //  
2988 {
2989   if(line->NbPoints() ==1 ) { 
2990     STATIC_BLOCAGE_SUR_PAS_TROP_GRAND=STATIC_PRECEDENT_INFLEXION=0;
2991   }
2992
2993   IntWalk_StatusDeflection aStatus = IntWalk_OK;
2994   Standard_Real FlecheCourante , Ratio = 1.0;
2995
2996   // Caro1 and Caro2
2997   const Handle(Adaptor3d_HSurface)& Caro1 = myIntersectionOn2S.Function().AuxillarSurface1();
2998   const Handle(Adaptor3d_HSurface)& Caro2 = myIntersectionOn2S.Function().AuxillarSurface2();
2999
3000   const IntSurf_PntOn2S& CurrentPoint = myIntersectionOn2S.Point(); 
3001   //==================================================================================
3002   //=========               S t o p   o n   p o i n t                 ============
3003   //================================================================================== 
3004   if (myIntersectionOn2S.IsTangent())  { 
3005     return IntWalk_ArretSurPoint;  
3006   }
3007
3008   const gp_Dir& TgCourante = myIntersectionOn2S.Direction();
3009
3010   const Standard_Real aCosBetweenTangent = TgCourante.Dot(previousd);
3011
3012   //==================================================================================
3013   //=========   R i s k   o f    i n f l e x i o n   p o i n t  ============
3014   //==================================================================================  
3015   if (aCosBetweenTangent < 0) {
3016     //------------------------------------------------------------
3017     //-- Risk of inflexion point : Divide the step by 2
3018     //-- Initialize STATIC_PRECEDENT_INFLEXION so that 
3019     //-- at the next call to return Pas_OK if there is no 
3020     //-- more risk of the point of inflexion
3021     //------------------------------------------------------------
3022
3023     pasuv[0]*=0.5;
3024     pasuv[1]*=0.5;
3025     pasuv[2]*=0.5;
3026     pasuv[3]*=0.5;
3027     STATIC_PRECEDENT_INFLEXION+=3; 
3028     if (pasuv[0] < ResoU1 && pasuv[1] <ResoV1 && pasuv[2] <ResoU2 && pasuv[3] < ResoV2)
3029       return IntWalk_ArretSurPointPrecedent;
3030     else 
3031       return IntWalk_PasTropGrand;
3032   }
3033   else {
3034     if(STATIC_PRECEDENT_INFLEXION  > 0) { 
3035       STATIC_PRECEDENT_INFLEXION -- ;
3036       return IntWalk_OK;
3037     }
3038   }
3039
3040   //==================================================================================
3041   //=========  D e t e c t    c o n f u s e d    P o in t s       ===========
3042   //==================================================================================
3043
3044   const Standard_Real aSqDist = previousPoint.Value().
3045     SquareDistance(CurrentPoint.Value());
3046
3047
3048   if (aSqDist < Precision::SquareConfusion()) { 
3049     pasInit[0] = Max(pasInit[0], 5.0*ResoU1);
3050     pasInit[1] = Max(pasInit[1], 5.0*ResoV1);
3051     pasInit[2] = Max(pasInit[2], 5.0*ResoU2);
3052     pasInit[3] = Max(pasInit[3], 5.0*ResoV2);
3053
3054     for(Standard_Integer i = 0; i < 4; i++)
3055     {
3056       pasuv[i] = Max(pasuv[i], Min(1.5*pasuv[i], pasInit[i]));
3057     }
3058     //Compute local resolution: for OCC26717
3059     if (Abs(pasuv[choixIso] - pasInit[choixIso]) <= Precision::Confusion())
3060     {
3061       Standard_Real CurU, CurV;
3062       if (choixIso == IntImp_UIsoparametricOnCaro1 ||
3063           choixIso == IntImp_VIsoparametricOnCaro1)
3064         previousPoint.ParametersOnS1(CurU, CurV);
3065       else
3066         previousPoint.ParametersOnS2(CurU, CurV);
3067       gp_Pnt CurPnt = (choixIso == IntImp_UIsoparametricOnCaro1 ||
3068                        choixIso == IntImp_VIsoparametricOnCaro1)?
3069         Adaptor3d_HSurfaceTool::Value(Caro1, CurU, CurV) :
3070         Adaptor3d_HSurfaceTool::Value(Caro2, CurU, CurV);
3071       gp_Pnt OffsetPnt;
3072       switch(choixIso)
3073       {
3074       case IntImp_UIsoparametricOnCaro1:
3075         OffsetPnt =
3076           Adaptor3d_HSurfaceTool::Value(Caro1,
3077                                         CurU + sensCheminement*pasuv[0],
3078                                         CurV);
3079         break;
3080       case IntImp_VIsoparametricOnCaro1:
3081         OffsetPnt =
3082           Adaptor3d_HSurfaceTool::Value(Caro1,
3083                                         CurU,
3084                                         CurV + sensCheminement*pasuv[1]);
3085         break;
3086       case IntImp_UIsoparametricOnCaro2:
3087         OffsetPnt =
3088           Adaptor3d_HSurfaceTool::Value(Caro2,
3089                                         CurU + sensCheminement*pasuv[2],
3090                                         CurV);
3091         break;
3092       case IntImp_VIsoparametricOnCaro2:
3093         OffsetPnt =
3094           Adaptor3d_HSurfaceTool::Value(Caro2,
3095                                         CurU,
3096                                         CurV + sensCheminement*pasuv[3]);
3097         break;
3098       default:break;
3099       }
3100       Standard_Real RefDist = CurPnt.Distance(OffsetPnt);
3101       Standard_Real LocalResol = 0.;
3102       if (RefDist > gp::Resolution())
3103         LocalResol = pasuv[choixIso] * tolconf / RefDist;
3104       if (pasuv[choixIso] < 2*LocalResol)
3105         pasuv[choixIso] = pasInit[choixIso] = 2*LocalResol;
3106     }
3107     ////////////////////////////////////////
3108     aStatus = IntWalk_PointConfondu;
3109   }
3110
3111   //==================================================================================
3112   Standard_Real Up1,Vp1,Uc1,Vc1,Du1,Dv1,AbsDu1,AbsDu2,AbsDv1,AbsDv2;
3113   Standard_Real Up2,Vp2,Uc2,Vc2,Du2,Dv2;
3114
3115   previousPoint.Parameters(Up1,Vp1,Up2,Vp2);
3116   CurrentPoint.Parameters(Uc1,Vc1,Uc2,Vc2);               
3117
3118   Du1 = Uc1 - Up1;   Dv1 = Vc1 - Vp1;
3119   Du2 = Uc2 - Up2;   Dv2 = Vc2 - Vp2;
3120
3121   AbsDu1 = Abs(Du1);
3122   AbsDu2 = Abs(Du2);
3123   AbsDv1 = Abs(Dv1);
3124   AbsDv2 = Abs(Dv2);
3125   //=================================================================================
3126   //====   S t e p   o f   p  r o g r e s s i o n (between previous and Current)   =======
3127   //=================================================================================
3128   if (   AbsDu1 < ResoU1 && AbsDv1 < ResoV1 
3129     && AbsDu2 < ResoU2 && AbsDv2 < ResoV2) {
3130       pasuv[0] = ResoU1; pasuv[1] = ResoV1; pasuv[2] = ResoU2; pasuv[3] = ResoV2;
3131       return(IntWalk_ArretSurPointPrecedent);
3132   }
3133   //==================================================================================
3134
3135   Standard_Real tolArea = 100.0;
3136   if (ResoU1 < Precision::PConfusion() ||
3137     ResoV1 < Precision::PConfusion() ||
3138     ResoU2 < Precision::PConfusion() ||
3139     ResoV2 < Precision::PConfusion() )
3140     tolArea =  tolArea*2.0;
3141
3142   Standard_Real Cosi1, CosRef1, Ang1, AngRef1, ResoUV1, Duv1, d1, tolCoeff1;   
3143   Standard_Real Cosi2, CosRef2, Ang2, AngRef2, ResoUV2, Duv2, d2, tolCoeff2;   
3144   Cosi1 = Du1*previousd1.X() + Dv1*previousd1.Y();
3145   Cosi2 = Du2*previousd2.X() + Dv2*previousd2.Y();
3146   Duv1 = Du1*Du1 + Dv1*Dv1;
3147   Duv2 = Du2*Du2 + Dv2*Dv2;
3148   ResoUV1 = ResoU1*ResoU1 + ResoV1*ResoV1;
3149   ResoUV2 = ResoU2*ResoU2 + ResoV2*ResoV2;
3150   //
3151   //modified by NIZNHY-PKV Wed Nov 13 12:25:44 2002 f
3152   //
3153   Standard_Real aMinDiv2=Precision::Confusion();
3154   aMinDiv2=aMinDiv2*aMinDiv2;
3155   //
3156   d1=d;
3157   if (Duv1>aMinDiv2)  {
3158     d1 = Abs(ResoUV1/Duv1);
3159     d1 = Min(Sqrt(d1)*tolArea, d);  
3160   } 
3161   //d1 = Abs(ResoUV1/Duv1); 
3162   //d1 = Min(Sqrt(d1)*tolArea,d);  
3163   //modified by NIZNHY-PKV Wed Nov 13 12:34:30 2002 t
3164   tolCoeff1 = Exp(d1);
3165   //
3166   //modified by NIZNHY-PKV Wed Nov 13 12:34:43 2002 f
3167   d2=d;
3168   if (Duv2>aMinDiv2) {
3169     d2 = Abs(ResoUV2/Duv2); 
3170     d2 = Min(Sqrt(d2)*tolArea,d); 
3171   }
3172   //d2 = Abs(ResoUV2/Duv2); 
3173   //d2 = Min(Sqrt(d2)*tolArea,d);  
3174   //modified by NIZNHY-PKV Wed Nov 13 12:34:53 2002 t
3175   tolCoeff2 = Exp(d2);
3176   CosRef1 = CosRef2D/tolCoeff1;
3177   CosRef2 = CosRef2D/tolCoeff2;
3178   //
3179   //==================================================================================
3180   //== The points are not confused :                                           ==
3181   //== D e t e c t    t h e   S t o p   a  t   p r e v i o u s  p o i n t ==
3182   //==                           N o t    T o o    G r e a t (angle in space UV)    ==
3183   //==                           C h a n g e    o f    s i d e                ==
3184   //==================================================================================
3185   if (aStatus != IntWalk_PointConfondu) {
3186     if(Cosi1*Cosi1 < CosRef1*Duv1 || Cosi2*Cosi2 < CosRef2*Duv2) {
3187       pasuv[0]*=0.5;  pasuv[1]*=0.5;  pasuv[2]*=0.5;  pasuv[3]*=0.5;
3188       if (pasuv[0]<ResoU1 && pasuv[1]<ResoV1 && pasuv[2]<ResoU2 && pasuv[3]<ResoV2) { 
3189         return(IntWalk_ArretSurPointPrecedent);
3190       }
3191       else {
3192         pasuv[0]*=0.5; pasuv[1]*=0.5; pasuv[2]*=0.5; pasuv[3]*=0.5;
3193         return(IntWalk_PasTropGrand);
3194       }
3195     }
3196     const gp_Dir2d& Tg2dcourante1 = myIntersectionOn2S.DirectionOnS1();
3197     const gp_Dir2d& Tg2dcourante2 = myIntersectionOn2S.DirectionOnS2();
3198     Cosi1 = Du1*Tg2dcourante1.X() + Dv1*Tg2dcourante1.Y();
3199     Cosi2 = Du2*Tg2dcourante2.X() + Dv2*Tg2dcourante2.Y();
3200     Ang1 = Abs(previousd1.Angle(Tg2dcourante1));  
3201     Ang2 = Abs(previousd2.Angle(Tg2dcourante2));  
3202     AngRef1 = AngRef2D*tolCoeff1;
3203     AngRef2 = AngRef2D*tolCoeff2;
3204     //-------------------------------------------------------
3205     //-- Test : Angle too great in space UV       -----
3206     //--        Change of  side                      -----
3207     //-------------------------------------------------------
3208     if(Cosi1*Cosi1 < CosRef1*Duv1 || Cosi2*Cosi2 < CosRef2*Duv2 || Ang1 > AngRef1 || Ang2 > AngRef2) {
3209       pasuv[0]*=0.5;  pasuv[1]*=0.5;  pasuv[2]*=0.5;  pasuv[3]*=0.5;
3210       if (pasuv[0]<ResoU1 && pasuv[1]<ResoV1 && pasuv[2]<ResoU2 && pasuv[3]<ResoV2) 
3211         return(IntWalk_ArretSurPoint);
3212       else 
3213         return(IntWalk_PasTropGrand);
3214     }
3215   }
3216   //<-OCC431(apo)
3217   //==================================================================================
3218   //== D e t e c t i o n   o f    :  Step Too Small 
3219   //==                               STEP TOO Great 
3220   //==================================================================================
3221
3222   //---------------------------------------
3223   //-- Estimate of the vector           --
3224   //---------------------------------------
3225   FlecheCourante =
3226     Sqrt(Abs((previousd.XYZ()-TgCourante.XYZ()).SquareModulus()*aSqDist))/8.;
3227
3228   if ( FlecheCourante<= fleche*0.5) {     //-- Current step too small
3229     if(FlecheCourante>1e-16) { 
3230       Ratio = 0.5*(fleche/FlecheCourante);
3231     }
3232     else { 
3233       Ratio = 10.0;
3234     }
3235     Standard_Real pasSu1 = pasuv[0];
3236     Standard_Real pasSv1 = pasuv[1];
3237     Standard_Real pasSu2 = pasuv[2];
3238     Standard_Real pasSv2 = pasuv[3];
3239
3240     //-- In  case if 
3241     //-- a point at U+DeltaU is required, ....
3242     //-- return a point at U + Epsilon
3243     //-- Epsilon << DeltaU.
3244
3245     if(pasuv[0]< AbsDu1) pasuv[0] = AbsDu1;
3246     if(pasuv[1]< AbsDv1) pasuv[1] = AbsDv1;
3247     if(pasuv[2]< AbsDu2) pasuv[2] = AbsDu2;
3248     if(pasuv[3]< AbsDv2) pasuv[3] = AbsDv2;
3249
3250     if(pasuv[0]<ResoU1) pasuv[0]=ResoU1;
3251     if(pasuv[1]<ResoV1) pasuv[1]=ResoV1;
3252     if(pasuv[2]<ResoU2) pasuv[2]=ResoU2;
3253     if(pasuv[3]<ResoV2) pasuv[3]=ResoV2;
3254     //-- if(Ratio>10.0 ) { Ratio=10.0; } 
3255     Standard_Real R1,R = pasInit[0]/pasuv[0];
3256     R1= pasInit[1]/pasuv[1];     if(R1<R) R=R1;
3257     R1= pasInit[2]/pasuv[2];     if(R1<R) R=R1;
3258     R1= pasInit[3]/pasuv[3];     if(R1<R) R=R1;
3259     if(Ratio > R) Ratio=R;
3260     pasuv[0] = Min(Ratio*pasuv[0],pasInit[0]);
3261     pasuv[1] = Min(Ratio*pasuv[1],pasInit[1]);
3262     pasuv[2] = Min(Ratio*pasuv[2],pasInit[2]);
3263     pasuv[3] = Min(Ratio*pasuv[3],pasInit[3]);
3264     if (pasuv[0] != pasSu1 || pasuv[2] != pasSu2|| 
3265       pasuv[1] != pasSv1 || pasuv[3] != pasSv2) {
3266         if(++STATIC_BLOCAGE_SUR_PAS_TROP_GRAND > 5) {
3267           STATIC_BLOCAGE_SUR_PAS_TROP_GRAND = 0;
3268           return IntWalk_PasTropGrand; 
3269         }
3270     }
3271     if(aStatus == IntWalk_OK) {
3272       STATIC_BLOCAGE_SUR_PAS_TROP_GRAND=0;
3273       //-- Try to increase the step
3274     }
3275     return aStatus;
3276   }
3277   else {                                //-- CurrentVector > vector*0.5 
3278     if (FlecheCourante > fleche) {      //-- Current step too Great
3279       Ratio = fleche/FlecheCourante; 
3280       pasuv[0] = Ratio*pasuv[0];
3281       pasuv[1] = Ratio*pasuv[1];
3282       pasuv[2] = Ratio*pasuv[2];
3283       pasuv[3] = Ratio*pasuv[3];
3284       //if(++STATIC_BLOCAGE_SUR_PAS_TROP_GRAND > 5) {
3285       //        STATIC_BLOCAGE_SUR_PAS_TROP_GRAND = 0;
3286       return IntWalk_PasTropGrand; 
3287       //}
3288     }
3289     else {                             //-- vector/2  <  CurrentVector <= vector   
3290       Ratio = 0.75 * (fleche / FlecheCourante);
3291     }
3292   }
3293
3294   if(aStatus != IntWalk_PointConfondu)
3295   {
3296     //Here, aCosBetweenTangent >= 0.0 definitely.
3297
3298     /*
3299     Brief algorithm description.
3300     We have two (not-coincindent) intersection point (P1 and P2). In every point,
3301     vector of tangent (to the intersection curve) is known (vectors T1 and T2).
3302     Basing on these data, we create osculating circle.
3303
3304                                     * - arc of osculating circle
3305                                 *      * 
3306                            P1 x----------x P2
3307                              /            \
3308                             /              \
3309                           Vec(T1)         Vec(T2)
3310
3311     Let me draw your attention to the following facts:
3312       1. Vectors T1 and T2 direct FROM (not TO) points P1 and P2. Therefore,
3313         one of previously computed vector should be reversed.
3314
3315     In this case, the absolute (!) value of the deflection between the arc of
3316     the osculating circle and the P1P2 segment can be computed as follows:
3317           e = d*(1-sin(B/2))/(2*cos(B/2)),      (1)
3318     where d is the length of P1P2 segment, B is the angle between vectors T1 and T2.
3319     At that,
3320               pi/2 <= B <= pi,
3321               cos(B/2) >= 0,
3322               sin(B/2) > 0,
3323               sin(B) > 0,
3324               cos(B) < 0.
3325
3326     Later, the normal state of algorithm work is (as we apply) 
3327               tolconf/2 <= e <= tolconf.
3328     In this case, we keep previous step.
3329
3330     If e < tolconf/2 then the local curvature of the intersection curve is small.
3331     As result, the step should be increased.
3332
3333     If e > tolconf then the step is too big. Therefore, we should decrease one.
3334
3335     Condition (1) is equivalent to
3336               sin(B/2) = 1 - 2/(1+(d/(2*e))^2) = Fs(e),
3337               cos(B) = 1 - 2*Fs(e)^2 = Fd(e),
3338     where Fs(e)and Fd(e) are some function with parameter "deflection".
3339
3340     Let mean that Fs(e) is decreasing function. Fd(e) is increasing function,
3341     in the range, where Fs(e) > 0.0 (i.e. when e < d/2).
3342
3343     Now, let substitute required deflection (tolconf or tolconf/2) to z. Then
3344     it is necessary to check if e < z or if e > z.
3345
3346     In this case, it is enough to comapare Fs(e) and Fs(z).
3347     At that Fs(e) > 0 because sin(B/2) > 0 always.
3348
3349     Therefore, if Fs(z) < 0.0 then Fs(e) > Fs(z) ==> e < z definitely.
3350     If Fs(z) > 0.0 then we can compare Fs(z)^2 and Fs(e)^2 or, in substance,
3351     values Fd(e) and Fd(z). If Fd(e) > Fd(z) then e > z and vice versa.
3352     */
3353
3354     //Fd(e) is already known (Fd(e) == -aCosBetweenTangent)
3355
3356     const Standard_Real anInvSqAbsArcDeflMax = 0.25*aSqDist/(tolconf*tolconf);
3357     const Standard_Real aSinB2Max = 1.0 - 2.0/(1.0 + anInvSqAbsArcDeflMax);
3358
3359     if(aSinB2Max >= 0.0 && (aCosBetweenTangent <= 2.0 * aSinB2Max * aSinB2Max - 1.0))
3360     {//Real deflection is greater or equal than tolconf
3361       aStatus = IntWalk_PasTropGrand;
3362     }
3363     else
3364     {//Real deflection is less than tolconf
3365       const Standard_Real anInvSqAbsArcDeflMin = 4.0*anInvSqAbsArcDeflMax;
3366       const Standard_Real aSinB2Min = 1.0 - 2.0/(1.0 + anInvSqAbsArcDeflMin);
3367
3368       if((aSinB2Min < 0.0) || (aCosBetweenTangent >= 2.0 * aSinB2Min * aSinB2Min - 1.0))
3369       {//Real deflection is less than tolconf/2.0
3370         aStatus = IntWalk_StepTooSmall;
3371       }
3372     }
3373
3374     if(aStatus == IntWalk_PasTropGrand)
3375     {
3376       pasuv[0]*=0.5; pasuv[1]*=0.5; pasuv[2]*=0.5; pasuv[3]*=0.5;
3377       return aStatus;
3378     }
3379
3380     if(aStatus == IntWalk_StepTooSmall)
3381     {
3382       pasuv[0] = Max(pasuv[0], AbsDu1);
3383       pasuv[1] = Max(pasuv[1], AbsDv1);
3384       pasuv[2] = Max(pasuv[2], AbsDu2);
3385       pasuv[3] = Max(pasuv[3], AbsDv2);
3386
3387       pasInit[0] = Max(pasInit[0], AbsDu1);
3388       pasInit[1] = Max(pasInit[1], AbsDv1);
3389       pasInit[2] = Max(pasInit[2], AbsDu2);
3390       pasInit[3] = Max(pasInit[3], AbsDv2);
3391
3392       return aStatus;
3393     }
3394   }
3395
3396   pasuv[0] = Max(myStepMin[0],Min(Min(Ratio*AbsDu1,pasuv[0]),pasInit[0]));
3397   pasuv[1] = Max(myStepMin[1],Min(Min(Ratio*AbsDv1,pasuv[1]),pasInit[1]));
3398   pasuv[2] = Max(myStepMin[2],Min(Min(Ratio*AbsDu2,pasuv[2]),pasInit[2]));
3399   pasuv[3] = Max(myStepMin[3],Min(Min(Ratio*AbsDv2,pasuv[3]),pasInit[3]));
3400
3401   if(aStatus == IntWalk_OK) STATIC_BLOCAGE_SUR_PAS_TROP_GRAND=0;
3402   return aStatus;
3403 }
3404
3405 Standard_Boolean IntWalk_PWalking::
3406 TestArret(const Standard_Boolean DejaReparti,
3407           TColStd_Array1OfReal& Param,
3408           IntImp_ConstIsoparametric&  ChoixIso)
3409
3410           //
3411           // test if the point of intersection set by these parameters remains in the 
3412           // natural domain of each square.
3413           // if the point outpasses reframe to find the best iso (border)
3414           // that intersects easiest the other square
3415           // otherwise test if closed line is present  
3416           // 
3417 {
3418   Standard_Real Uvd[4],Uvf[4],Epsuv[4],Duv[4],Uvp[4],dv,dv2,ParC[4];
3419   Standard_Real DPc,DPb;
3420   Standard_Integer i = 0, k = 0;
3421   Epsuv[0] = ResoU1;
3422   Epsuv[1] = ResoV1;
3423   Epsuv[2] = ResoU2;
3424   Epsuv[3] = ResoV2;
3425   previousPoint.Parameters(Uvp[0],Uvp[1],Uvp[2],Uvp[3]);
3426
3427   Standard_Real SolParam[4];
3428   myIntersectionOn2S.Point().Parameters(SolParam[0],SolParam[1],SolParam[2],SolParam[3]);
3429
3430   Standard_Boolean Trouve = Standard_False;
3431
3432   Uvd[0]=Um1;   Uvf[0]=UM1;   Uvd[1]=Vm1;   Uvf[1]=VM1;
3433   Uvd[2]=Um2;   Uvf[2]=UM2;   Uvd[3]=Vm2;   Uvf[3]=VM2;
3434
3435   Standard_Integer im1;
3436   for ( i = 1,im1 = 0;i<=4;i++,im1++) {
3437     switch(i) { 
3438     case 1: k=2; break;
3439     case 2: k=1; break;
3440     case 3: k=4; break;
3441     case 4: k=3; break;
3442     }
3443     if (Param(i) < (Uvd[im1]-Epsuv[im1]) ||
3444       SolParam[im1] < (Uvd[im1]-Epsuv[im1]))     //--     Current -----  Bound Inf -----  Previous
3445     {
3446       Trouve    = Standard_True;                   //-- 
3447       DPc       = Uvp[im1]-Param(i);               //--     Previous  - Current
3448       DPb       = Uvp[im1]-Uvd[im1];               //--     Previous  - Bound Inf
3449       ParC[im1] = Uvd[im1];                        //--     ParamCorrige
3450       dv        = Param(k)-Uvp[k-1];               //--     Current   - Previous (other Direction)
3451       dv2       = dv*dv;         
3452       if(dv2>RealEpsilon()) {                       //--    Progress at the other Direction ?
3453         Duv[im1]  = DPc*DPb + dv2;
3454         Duv[im1]  = Duv[im1]*Duv[im1]/(DPc*DPc+dv2)/(DPb*DPb+dv2);
3455       }
3456       else {
3457         Duv[im1]=-1.0;                              //--    If no progress, do not change  
3458       }                                             //--    the choice of iso 
3459     }   
3460     else if (Param(i) > (Uvf[im1] + Epsuv[im1]) ||
3461       SolParam[im1] > (Uvf[im1] + Epsuv[im1]))//--    Previous -----  Bound Sup -----  Current
3462     {
3463       Trouve    = Standard_True;                    //-- 
3464       DPc       = Param(i)-Uvp[im1];                //--     Current   - Previous
3465       DPb       = Uvf[im1]-Uvp[im1];                //--     Bound Sup - Previous 
3466       ParC[im1] = Uvf[im1];                         //--     Param Corrige
3467       dv        = Param(k)-Uvp[k-1];                //--     Current   - Previous (other Direction)
3468       dv2       = dv*dv;
3469       if(dv2>RealEpsilon()) {                       //--     Progress in other Direction ?
3470         Duv[im1]  =  DPc*DPb + dv2;
3471         Duv[im1]  = Duv[im1]*Duv[im1]/(DPc*DPc+dv2)/(DPb*DPb+dv2);
3472       }
3473       else {
3474         Duv[im1]=-1.0;                              //--    If no progress, do not change 
3475       }                                             //--    the choice of iso 
3476     }
3477     else { 
3478       Duv[im1]= -1.;
3479       ParC[im1]=Param(i);
3480     }
3481   }
3482
3483   if (Trouve) {
3484     //--------------------------------------------------
3485     //-- One of Parameters u1,v1,u2,v2 is outside of  --
3486     //-- the natural limits.                          -- 
3487     //-- Find the best direction of                   -- 
3488     //-- progress and reframe the parameters.        --
3489     //--------------------------------------------------
3490     Standard_Real ddv = -1.0;
3491     k=-1;
3492     for (i=0;i<=3;i++) {
3493       Param(i+1) = ParC[i];
3494       if(Duv[i]>ddv) { 
3495         ddv = Duv[i];
3496         k=i;
3497       }
3498     }
3499     if(k!=-1) { 
3500       ChoixIso   = ChoixRef[k];
3501     }
3502     else { 
3503       if((ParC[0]<=Uvd[0]+Epsuv[0]) || (ParC[0]>=Uvf[0]-Epsuv[0])) {
3504         ChoixIso = IntImp_UIsoparametricOnCaro1;
3505       }
3506       else if((ParC[1]<=Uvd[1]+Epsuv[1]) || (ParC[1]>=Uvf[1]-Epsuv[1])) {
3507         ChoixIso = IntImp_VIsoparametricOnCaro1;
3508       }
3509       else if((ParC[2]<=Uvd[2]+Epsuv[2]) || (ParC[2]>=Uvf[2]-Epsuv[2])) {
3510         ChoixIso = IntImp_UIsoparametricOnCaro2;
3511       }
3512       else if((ParC[3]<=Uvd[3]+Epsuv[3]) || (ParC[3]>=Uvf[3]-Epsuv[3])) {
3513         ChoixIso = IntImp_VIsoparametricOnCaro2;
3514       }
3515     }
3516     close = Standard_False;
3517     return Standard_True;
3518   }
3519   else 
3520   {  
3521     if (!DejaReparti) { // find if line closed
3522
3523       Standard_Real u,v;
3524       const IntSurf_PntOn2S& POn2S1=line->Value(1);
3525       //On S1
3526       POn2S1.ParametersOnS1(u,v);
3527       gp_Pnt2d P1uvS1(u,v);
3528       previousPoint.ParametersOnS1(u,v);
3529       gp_Pnt2d PrevuvS1(u,v);
3530       myIntersectionOn2S.Point().ParametersOnS1(u,v);
3531       gp_Pnt2d myIntersuvS1(u,v);
3532       Standard_Boolean close2dS1 = (P1uvS1.XY()-PrevuvS1.XY())*
3533         (P1uvS1.XY()-myIntersuvS1.XY()) < 0.0;
3534       //On S2
3535       POn2S1.ParametersOnS2(u,v);
3536       gp_Pnt2d P1uvS2(u,v);
3537       previousPoint.ParametersOnS2(u,v);
3538       gp_Pnt2d PrevuvS2(u,v);
3539       myIntersectionOn2S.Point().ParametersOnS2(u,v);
3540       gp_Pnt2d myIntersuvS2(u,v);
3541       Standard_Boolean close2dS2 = (P1uvS2.XY()-PrevuvS2.XY())*
3542         (P1uvS2.XY()-myIntersuvS2.XY()) < 0.0;
3543
3544       close = close2dS1 && close2dS2;
3545       return close;
3546     }
3547     else return Standard_False;
3548   }
3549 }
3550