987bf634a5e16bebb4c620c53cc0f30931dbcfaa
[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           break;
2431       }
2432
2433       for (++aPInd; aPInd <= aNbPnts; aPInd++)
2434       {
2435         aP2.SetXYZ(line->Value(aPInd).Value().XYZ());
2436         if (aP1.SquareDistance(aP2) > Precision::SquareConfusion())
2437           break;
2438       }
2439
2440       if (aPInd > aNbPnts)
2441       {
2442         return isOK;
2443       }
2444
2445       const gp_XYZ aDir01(aP1.XYZ() - aPInt.XYZ());
2446       const gp_XYZ aDir12(aP2.XYZ() - aP1.XYZ());
2447
2448       if (aDir01.Dot(aDir12) > 0.0)
2449       {
2450         break;
2451       }
2452
2453       RemoveAPoint(1);
2454     }
2455
2456     line->InsertBefore(1, anIP);
2457     isOK = Standard_True;
2458   }
2459   else
2460   {
2461     while (line->NbPoints() > 1)
2462     {
2463       const Standard_Integer aNbPnts = line->NbPoints();
2464
2465       gp_Pnt aPPrev, aPCurr;
2466       Standard_Integer aPInd = aNbPnts;
2467       for (; aPInd > 0; aPInd--)
2468       {
2469         aPCurr.SetXYZ(line->Value(aPInd).Value().XYZ());
2470         if (aPCurr.SquareDistance(aPInt) > Precision::SquareConfusion())
2471           break;
2472       }
2473
2474       for (--aPInd; aPInd > 0; aPInd--)
2475       {
2476         aPPrev.SetXYZ(line->Value(aPInd).Value().XYZ());
2477         if (aPCurr.SquareDistance(aPPrev) > Precision::SquareConfusion())
2478           break;
2479       }
2480
2481       if (aPInd < 1)
2482       {
2483         return isOK;
2484       }
2485
2486       const gp_XYZ aDirPC(aPCurr.XYZ() - aPPrev.XYZ());
2487       const gp_XYZ aDirCN(aPInt.XYZ() - aPCurr.XYZ());
2488
2489       if (aDirPC.Dot(aDirCN) > 0.0)
2490       {
2491         break;
2492       }
2493
2494       RemoveAPoint(aNbPnts);
2495     }
2496
2497     line->Add(anIP);
2498     isOK = Standard_True;
2499   }
2500
2501   return isOK;
2502 }
2503
2504 //=======================================================================
2505 //function : PutToBoundary
2506 //purpose  : 
2507 //=======================================================================
2508 Standard_Boolean IntWalk_PWalking::
2509 PutToBoundary(const Handle(Adaptor3d_HSurface)& theASurf1,
2510               const Handle(Adaptor3d_HSurface)& theASurf2)
2511 {
2512   const Standard_Real aTolMin = Precision::Confusion();
2513
2514   Standard_Boolean hasBeenAdded = Standard_False;
2515
2516   const Standard_Real aU1bFirst = theASurf1->FirstUParameter();
2517   const Standard_Real aU1bLast = theASurf1->LastUParameter();
2518   const Standard_Real aU2bFirst = theASurf2->FirstUParameter();
2519   const Standard_Real aU2bLast = theASurf2->LastUParameter();
2520   const Standard_Real aV1bFirst = theASurf1->FirstVParameter();
2521   const Standard_Real aV1bLast = theASurf1->LastVParameter();
2522   const Standard_Real aV2bFirst = theASurf2->FirstVParameter();
2523   const Standard_Real aV2bLast = theASurf2->LastVParameter();
2524
2525   Standard_Real aTol = 1.0;
2526   aTol = Min(aTol, aU1bLast - aU1bFirst);
2527   aTol = Min(aTol, aU2bLast - aU2bFirst);
2528   aTol = Min(aTol, aV1bLast - aV1bFirst);
2529   aTol = Min(aTol, aV2bLast - aV2bFirst)*1.0e-3;
2530
2531   if(aTol <= 2.0*aTolMin)
2532     return hasBeenAdded;
2533
2534   Standard_Boolean isNeedAdding = Standard_False;
2535   Standard_Boolean isU1parallel = Standard_False, isV1parallel = Standard_False;
2536   Standard_Boolean isU2parallel = Standard_False, isV2parallel = Standard_False;
2537   IsParallel(line, Standard_True, aTol, isU1parallel, isV1parallel);
2538   IsParallel(line, Standard_False, aTol, isU2parallel, isV2parallel);
2539
2540   Standard_Real u1, v1, u2, v2;
2541   line->Value(1).Parameters(u1, v1, u2, v2);
2542   Standard_Real aDelta = 0.0;
2543
2544   if(!isV1parallel)
2545   {
2546     aDelta = u1 - aU1bFirst;
2547     if((aTolMin < aDelta) && (aDelta < aTol))
2548     {
2549       u1 = aU1bFirst;
2550       isNeedAdding = Standard_True;
2551     }
2552     else
2553     {
2554       aDelta = aU1bLast - u1;
2555       if((aTolMin < aDelta) && (aDelta < aTol))
2556       {
2557         u1 = aU1bLast;
2558         isNeedAdding = Standard_True;
2559       }
2560     }
2561   }
2562
2563   if(!isV2parallel)
2564   {
2565     aDelta = u2 - aU2bFirst;
2566     if((aTolMin < aDelta) && (aDelta < aTol))
2567     {
2568       u2 = aU2bFirst;
2569       isNeedAdding = Standard_True;
2570     }
2571     else
2572     {
2573       aDelta = aU2bLast - u2;
2574       if((aTolMin < aDelta) && (aDelta < aTol))
2575       {
2576         u2 = aU2bLast;
2577         isNeedAdding = Standard_True;
2578       }
2579     }
2580   }
2581
2582   if(!isU1parallel)
2583   {
2584     aDelta = v1 - aV1bFirst;
2585     if((aTolMin < aDelta) && (aDelta < aTol))
2586     {
2587       v1 = aV1bFirst;
2588       isNeedAdding = Standard_True;
2589     }
2590     else
2591     {
2592       aDelta = aV1bLast - v1;
2593       if((aTolMin < aDelta) && (aDelta < aTol))
2594       {
2595         v1 = aV1bLast;
2596         isNeedAdding = Standard_True;
2597       }
2598     }
2599   }
2600
2601   if(!isU2parallel)
2602   {
2603     aDelta = v2 - aV2bFirst;
2604     if((aTolMin < aDelta) && (aDelta < aTol))
2605     {
2606       v2 = aV2bFirst;
2607       isNeedAdding = Standard_True;
2608     }
2609     else
2610     {
2611       aDelta = aV2bLast - v2;
2612       if((aTolMin < aDelta) && (aDelta < aTol))
2613       {
2614         v2 = aV2bLast;
2615         isNeedAdding = Standard_True;
2616       }
2617     }
2618   }
2619
2620   if(isNeedAdding)
2621   {
2622     hasBeenAdded = 
2623       SeekPointOnBoundary(theASurf1, theASurf2, u1, 
2624       v1, u2, v2, Standard_True);
2625   }
2626
2627   const Standard_Integer aNbPnts = line->NbPoints();
2628   isNeedAdding = Standard_False;
2629   line->Value(aNbPnts).Parameters(u1, v1, u2, v2);
2630
2631   if(!isV1parallel)
2632   {
2633     aDelta = u1 - aU1bFirst;
2634     if((aTolMin < aDelta) && (aDelta < aTol))
2635     {
2636       u1 = aU1bFirst;
2637       isNeedAdding = Standard_True;
2638     }
2639     else
2640     {
2641       aDelta = aU1bLast - u1;
2642       if((aTolMin < aDelta) && (aDelta < aTol))
2643       {
2644         u1 = aU1bLast;
2645         isNeedAdding = Standard_True;
2646       }
2647     }
2648   }
2649
2650   if(!isV2parallel)
2651   {
2652     aDelta = u2 - aU2bFirst;
2653     if((aTolMin < aDelta) && (aDelta < aTol))
2654     {
2655       u2 = aU2bFirst;
2656       isNeedAdding = Standard_True;
2657     }
2658     else
2659     {
2660       aDelta = aU2bLast - u2;
2661       if((aTolMin < aDelta) && (aDelta < aTol))
2662       {
2663         u2 = aU2bLast;
2664         isNeedAdding = Standard_True;
2665       }
2666     }
2667   }
2668
2669   if(!isU1parallel)
2670   {
2671     aDelta = v1 - aV1bFirst;
2672     if((aTolMin < aDelta) && (aDelta < aTol))
2673     {
2674       v1 = aV1bFirst;
2675       isNeedAdding = Standard_True;
2676     }
2677     else
2678     {
2679       aDelta = aV1bLast - v1;
2680       if((aTolMin < aDelta) && (aDelta < aTol))
2681       {
2682         v1 = aV1bLast;
2683         isNeedAdding = Standard_True;
2684       }
2685     }
2686   }
2687
2688   if(!isU2parallel)
2689   {
2690     aDelta = v2 - aV2bFirst;
2691     if((aTolMin < aDelta) && (aDelta < aTol))
2692     {
2693       v2 = aV2bFirst;
2694       isNeedAdding = Standard_True;
2695     }
2696     else
2697     {
2698       aDelta = aV2bLast - v2;
2699       if((aTolMin < aDelta) && (aDelta < aTol))
2700       {
2701         v2 = aV2bLast;
2702         isNeedAdding = Standard_True;
2703       }
2704     }
2705   }
2706
2707   if(isNeedAdding)
2708   {
2709     hasBeenAdded = 
2710       SeekPointOnBoundary(theASurf1, theASurf2, u1, 
2711       v1, u2, v2, Standard_False);
2712   }
2713
2714   return hasBeenAdded;
2715 }
2716
2717 //=======================================================================
2718 //function : SeekAdditionalPoints
2719 //purpose  : 
2720 //=======================================================================
2721 Standard_Boolean IntWalk_PWalking::
2722 SeekAdditionalPoints( const Handle(Adaptor3d_HSurface)& theASurf1,
2723                      const Handle(Adaptor3d_HSurface)& theASurf2,
2724                      const Standard_Integer theMinNbPoints)
2725 {
2726   const Standard_Real aTol = 1.0e-14;
2727   Standard_Integer aNbPoints = line->NbPoints();
2728   if(aNbPoints > theMinNbPoints)
2729     return Standard_True;
2730
2731   const Standard_Real aU1bFirst = theASurf1->FirstUParameter();
2732   const Standard_Real aU1bLast = theASurf1->LastUParameter();
2733   const Standard_Real aU2bFirst = theASurf2->FirstUParameter();
2734   const Standard_Real aU2bLast = theASurf2->LastUParameter();
2735   const Standard_Real aV1bFirst = theASurf1->FirstVParameter();
2736   const Standard_Real aV1bLast = theASurf1->LastVParameter();
2737   const Standard_Real aV2bFirst = theASurf2->FirstVParameter();
2738   const Standard_Real aV2bLast = theASurf2->LastVParameter();
2739
2740
2741   Standard_Boolean isPrecise = Standard_False;
2742
2743   TColStd_Array1OfReal aPnt(1, 4);
2744   aPnt.Init(0.0);
2745
2746   Standard_Integer aNbPointsPrev = 0;
2747   while(aNbPoints < theMinNbPoints && (aNbPoints != aNbPointsPrev))
2748   {
2749     aNbPointsPrev = aNbPoints;
2750     for(Standard_Integer fp = 1, lp = 2; fp < aNbPoints; fp = lp + 1)
2751     {
2752       Standard_Real U1f, V1f, U2f, V2f; //first point in 1st and 2nd surafaces
2753       Standard_Real U1l, V1l, U2l, V2l; //last  point in 1st and 2nd surafaces
2754
2755       lp = fp+1;
2756       line->Value(fp).Parameters(U1f, V1f, U2f, V2f);
2757       line->Value(lp).Parameters(U1l, V1l, U2l, V2l);
2758
2759       aPnt(1) = 0.5*(U1f + U1l);
2760       if(aPnt(1) < aU1bFirst)
2761         aPnt(1) = aU1bFirst;
2762       if(aPnt(1) > aU1bLast)
2763         aPnt(1) = aU1bLast;
2764
2765       aPnt(2) = 0.5*(V1f+V1l);
2766       if(aPnt(2) < aV1bFirst)
2767         aPnt(2) = aV1bFirst;
2768       if(aPnt(2) > aV1bLast)
2769         aPnt(2) = aV1bLast;
2770
2771       aPnt(3) = 0.5*(U2f+U2l);
2772       if(aPnt(3) < aU2bFirst)
2773         aPnt(3) = aU2bFirst;
2774       if(aPnt(3) > aU2bLast)
2775         aPnt(3) = aU2bLast;
2776
2777       aPnt(4) = 0.5*(V2f+V2l);
2778       if(aPnt(4) < aV2bFirst)
2779         aPnt(4) = aV2bFirst;
2780       if(aPnt(4) > aV2bLast)
2781         aPnt(4) = aV2bLast;
2782
2783       Standard_Boolean aStatus = Standard_False;
2784       Standard_Integer aNbIter = 5;
2785       do
2786       {
2787         aStatus = DistanceMinimizeByGradient(theASurf1, theASurf2, aPnt);
2788         if(aStatus)
2789         {
2790           break;
2791         }
2792
2793         aStatus = DistanceMinimizeByExtrema(theASurf1, theASurf2->Value(aPnt(3), aPnt(4)), aPnt(1), aPnt(2));
2794         if(aStatus)
2795         {
2796           break;
2797         }
2798
2799         aStatus = DistanceMinimizeByExtrema(theASurf2, theASurf1->Value(aPnt(1), aPnt(2)), aPnt(3), aPnt(4));
2800         if(aStatus)
2801         {
2802           break;
2803         }
2804       }
2805       while(!aStatus && (--aNbIter > 0));
2806
2807       if(aStatus)
2808       {
2809         gp_Pnt  aP1 = theASurf1->Value(aPnt(1), aPnt(2)),
2810           aP2 = theASurf2->Value(aPnt(3), aPnt(4));
2811         gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));
2812
2813         const Standard_Real aSQDist1 = aPInt.SquareDistance(aP1),
2814           aSQDist2 = aPInt.SquareDistance(aP2);
2815
2816         if((aSQDist1 < aTol) && (aSQDist2 < aTol))
2817         {
2818           IntSurf_PntOn2S anIP;
2819           anIP.SetValue(aPInt, aPnt(1), aPnt(2), aPnt(3), aPnt(4));
2820           line->InsertBefore(lp, anIP);
2821
2822           isPrecise = Standard_True;
2823
2824           if(++aNbPoints >= theMinNbPoints)
2825             break;
2826         }
2827         else
2828         {
2829           lp--;
2830         }
2831       }
2832     }
2833   }
2834
2835   return isPrecise;
2836 }
2837
2838 void IntWalk_PWalking::
2839 RepartirOuDiviser(Standard_Boolean& DejaReparti,
2840                   IntImp_ConstIsoparametric& ChoixIso,
2841                   Standard_Boolean& Arrive) 
2842
2843                   // at the neighborhood of a point, there is a fail of marching 
2844                   // it is required to divide the steps to try to continue
2845                   // if the step is too small if we are on border
2846                   // restart in another direction if it was not done, otherwise stop
2847
2848 {
2849   //  Standard_Integer i;
2850   if (Arrive) {    //restart in the other direction
2851     if (!DejaReparti ) {
2852       Arrive = Standard_False;
2853       DejaReparti = Standard_True;
2854       previousPoint = line->Value(1);
2855       previoustg = Standard_False;
2856       previousd1 = firstd1;
2857       previousd2 = firstd2;
2858       previousd = tgdir;
2859       myTangentIdx = line->NbPoints();
2860       tgdir.Reverse();
2861       line->Reverse();
2862
2863       //-- printf("\nIntWalk_PWalking_2.gxx Reverse %3d\n",indextg);
2864       sensCheminement = -1;
2865       tgfirst      = tglast;
2866       tglast       = Standard_False;
2867       ChoixIso     = choixIsoSav;
2868 #if 0
2869       pasuv[0]=pasSav[0];
2870       pasuv[1]=pasSav[1];
2871       pasuv[2]=pasSav[2];
2872       pasuv[3]=pasSav[3];
2873 #else 
2874       Standard_Real u1,v1,u2,v2;
2875       Standard_Real U1,V1,U2,V2;
2876       Standard_Integer nn=line->NbPoints();
2877       if(nn>2) { 
2878         line->Value(nn).Parameters(u1,v1,u2,v2);
2879         line->Value(nn-1).Parameters(U1,V1,U2,V2);
2880         pasuv[0]=Abs(u1-U1);
2881         pasuv[1]=Abs(v1-V1);
2882         pasuv[2]=Abs(u2-U2);
2883         pasuv[3]=Abs(v2-V2);
2884       }
2885 #endif
2886
2887     }
2888   }  
2889   else  {
2890     if (    pasuv[0]*0.5 < ResoU1
2891       &&  pasuv[1]*0.5 < ResoV1
2892       &&  pasuv[2]*0.5 < ResoU2
2893       &&  pasuv[3]*0.5 < ResoV2
2894       ) {
2895         if (!previoustg) {
2896           tglast = Standard_True;      // IS IT ENOUGH ????
2897         }
2898
2899         if (!DejaReparti)
2900         {  //restart in the other direction
2901           DejaReparti = Standard_True;
2902           previousPoint = line->Value(1);
2903           previoustg = Standard_False;
2904           previousd1 = firstd1;
2905           previousd2 = firstd2;
2906           previousd = tgdir;
2907           myTangentIdx = line->NbPoints();
2908           tgdir.Reverse();
2909           line->Reverse();
2910
2911           //-- printf("\nIntWalk_PWalking_2.gxx Reverse %3d\n",indextg);
2912
2913           sensCheminement   = -1;
2914           tgfirst           = tglast;
2915           tglast            = Standard_False;
2916           ChoixIso          = choixIsoSav;
2917
2918 #if 0 
2919           pasuv[0]=pasSav[0];
2920           pasuv[1]=pasSav[1];
2921           pasuv[2]=pasSav[2];
2922           pasuv[3]=pasSav[3];
2923 #else 
2924           Standard_Real u1,v1,u2,v2;
2925           Standard_Real U1,V1,U2,V2;
2926           Standard_Integer nn=line->NbPoints();
2927           if(nn>2) { 
2928             line->Value(nn).Parameters(u1,v1,u2,v2);
2929             line->Value(nn-1).Parameters(U1,V1,U2,V2);
2930             pasuv[0]=Abs(u1-U1);
2931             pasuv[1]=Abs(v1-V1);
2932             pasuv[2]=Abs(u2-U2);
2933             pasuv[3]=Abs(v2-V2);
2934           }
2935 #endif
2936         }
2937         else Arrive = Standard_True;
2938     }
2939     else {
2940       pasuv[0]*=0.5;
2941       pasuv[1]*=0.5;
2942       pasuv[2]*=0.5;
2943       pasuv[3]*=0.5; 
2944     }
2945   }
2946 }
2947
2948 namespace {
2949   //OCC431(apo): modified ->
2950   static const Standard_Real CosRef2D =  Cos(M_PI/9.0),  AngRef2D = M_PI/2.0; 
2951
2952   static const Standard_Real d = 7.0;
2953 }
2954
2955 IntWalk_StatusDeflection  IntWalk_PWalking::TestDeflection(const IntImp_ConstIsoparametric choixIso)
2956
2957 // test if vector is observed by calculating an increase of vector 
2958 //     or the previous point and its tangent, the new calculated point and its  
2959 //     tangent; it is possible to find a cube passing by the 2 points and having as a 
2960 //     derivative the tangents of the intersection
2961 //     calculate the point with parameter 0.5 on cube=p1 
2962 //     calculate the medium point of 2 points of intersection=p2
2963 //   if arrow/2<=||p1p2||<= arrow consider that the vector is observed
2964 //   otherwise adjust the step depending on the ratio ||p1p2||/vector
2965 //   and the previous step 
2966 // test if in  2 tangent planes of surfaces there is no too great angle2d 
2967 // grand : if yes divide the step
2968 // test if there is no change of side
2969 //  
2970 {
2971   if(line->NbPoints() ==1 ) { 
2972     STATIC_BLOCAGE_SUR_PAS_TROP_GRAND=STATIC_PRECEDENT_INFLEXION=0;
2973   }
2974
2975   IntWalk_StatusDeflection aStatus = IntWalk_OK;
2976   Standard_Real FlecheCourante , Ratio = 1.0;
2977
2978   // Caro1 and Caro2
2979   const Handle(Adaptor3d_HSurface)& Caro1 = myIntersectionOn2S.Function().AuxillarSurface1();
2980   const Handle(Adaptor3d_HSurface)& Caro2 = myIntersectionOn2S.Function().AuxillarSurface2();
2981
2982   const IntSurf_PntOn2S& CurrentPoint = myIntersectionOn2S.Point(); 
2983   //==================================================================================
2984   //=========               S t o p   o n   p o i n t                 ============
2985   //================================================================================== 
2986   if (myIntersectionOn2S.IsTangent())  { 
2987     return IntWalk_ArretSurPoint;  
2988   }
2989
2990   const gp_Dir& TgCourante = myIntersectionOn2S.Direction();
2991
2992   const Standard_Real aCosBetweenTangent = TgCourante.Dot(previousd);
2993
2994   //==================================================================================
2995   //=========   R i s k   o f    i n f l e x i o n   p o i n t  ============
2996   //==================================================================================  
2997   if (aCosBetweenTangent < 0) {
2998     //------------------------------------------------------------
2999     //-- Risk of inflexion point : Divide the step by 2
3000     //-- Initialize STATIC_PRECEDENT_INFLEXION so that 
3001     //-- at the next call to return Pas_OK if there is no 
3002     //-- more risk of the point of inflexion
3003     //------------------------------------------------------------
3004
3005     pasuv[0]*=0.5;
3006     pasuv[1]*=0.5;
3007     pasuv[2]*=0.5;
3008     pasuv[3]*=0.5;
3009     STATIC_PRECEDENT_INFLEXION+=3; 
3010     if (pasuv[0] < ResoU1 && pasuv[1] <ResoV1 && pasuv[2] <ResoU2 && pasuv[3] < ResoV2)
3011       return IntWalk_ArretSurPointPrecedent;
3012     else 
3013       return IntWalk_PasTropGrand;
3014   }
3015   else {
3016     if(STATIC_PRECEDENT_INFLEXION  > 0) { 
3017       STATIC_PRECEDENT_INFLEXION -- ;
3018       return IntWalk_OK;
3019     }
3020   }
3021
3022   //==================================================================================
3023   //=========  D e t e c t    c o n f u s e d    P o in t s       ===========
3024   //==================================================================================
3025
3026   const Standard_Real aSqDist = previousPoint.Value().
3027     SquareDistance(CurrentPoint.Value());
3028
3029
3030   if (aSqDist < Precision::SquareConfusion()) { 
3031     pasInit[0] = Max(pasInit[0], 5.0*ResoU1);
3032     pasInit[1] = Max(pasInit[1], 5.0*ResoV1);
3033     pasInit[2] = Max(pasInit[2], 5.0*ResoU2);
3034     pasInit[3] = Max(pasInit[3], 5.0*ResoV2);
3035
3036     for(Standard_Integer i = 0; i < 4; i++)
3037     {
3038       pasuv[i] = Max(pasuv[i], Min(1.5*pasuv[i], pasInit[i]));
3039     }
3040     //Compute local resolution: for OCC26717
3041     if (Abs(pasuv[choixIso] - pasInit[choixIso]) <= Precision::Confusion())
3042     {
3043       Standard_Real CurU, CurV;
3044       if (choixIso == IntImp_UIsoparametricOnCaro1 ||
3045           choixIso == IntImp_VIsoparametricOnCaro1)
3046         previousPoint.ParametersOnS1(CurU, CurV);
3047       else
3048         previousPoint.ParametersOnS2(CurU, CurV);
3049       gp_Pnt CurPnt = (choixIso == IntImp_UIsoparametricOnCaro1 ||
3050                        choixIso == IntImp_VIsoparametricOnCaro1)?
3051         Adaptor3d_HSurfaceTool::Value(Caro1, CurU, CurV) :
3052         Adaptor3d_HSurfaceTool::Value(Caro2, CurU, CurV);
3053       gp_Pnt OffsetPnt;
3054       switch(choixIso)
3055       {
3056       case IntImp_UIsoparametricOnCaro1:
3057         OffsetPnt =
3058           Adaptor3d_HSurfaceTool::Value(Caro1,
3059                                         CurU + sensCheminement*pasuv[0],
3060                                         CurV);
3061         break;
3062       case IntImp_VIsoparametricOnCaro1:
3063         OffsetPnt =
3064           Adaptor3d_HSurfaceTool::Value(Caro1,
3065                                         CurU,
3066                                         CurV + sensCheminement*pasuv[1]);
3067         break;
3068       case IntImp_UIsoparametricOnCaro2:
3069         OffsetPnt =
3070           Adaptor3d_HSurfaceTool::Value(Caro2,
3071                                         CurU + sensCheminement*pasuv[2],
3072                                         CurV);
3073         break;
3074       case IntImp_VIsoparametricOnCaro2:
3075         OffsetPnt =
3076           Adaptor3d_HSurfaceTool::Value(Caro2,
3077                                         CurU,
3078                                         CurV + sensCheminement*pasuv[3]);
3079         break;
3080       default:break;
3081       }
3082       Standard_Real RefDist = CurPnt.Distance(OffsetPnt);
3083       Standard_Real LocalResol = 0.;
3084       if (RefDist > gp::Resolution())
3085         LocalResol = pasuv[choixIso] * tolconf / RefDist;
3086       if (pasuv[choixIso] < 2*LocalResol)
3087         pasuv[choixIso] = pasInit[choixIso] = 2*LocalResol;
3088     }
3089     ////////////////////////////////////////
3090     aStatus = IntWalk_PointConfondu;
3091   }
3092
3093   //==================================================================================
3094   Standard_Real Up1,Vp1,Uc1,Vc1,Du1,Dv1,AbsDu1,AbsDu2,AbsDv1,AbsDv2;
3095   Standard_Real Up2,Vp2,Uc2,Vc2,Du2,Dv2;
3096
3097   previousPoint.Parameters(Up1,Vp1,Up2,Vp2);
3098   CurrentPoint.Parameters(Uc1,Vc1,Uc2,Vc2);               
3099
3100   Du1 = Uc1 - Up1;   Dv1 = Vc1 - Vp1;
3101   Du2 = Uc2 - Up2;   Dv2 = Vc2 - Vp2;
3102
3103   AbsDu1 = Abs(Du1);
3104   AbsDu2 = Abs(Du2);
3105   AbsDv1 = Abs(Dv1);
3106   AbsDv2 = Abs(Dv2);
3107   //=================================================================================
3108   //====   S t e p   o f   p  r o g r e s s i o n (between previous and Current)   =======
3109   //=================================================================================
3110   if (   AbsDu1 < ResoU1 && AbsDv1 < ResoV1 
3111     && AbsDu2 < ResoU2 && AbsDv2 < ResoV2) {
3112       pasuv[0] = ResoU1; pasuv[1] = ResoV1; pasuv[2] = ResoU2; pasuv[3] = ResoV2;
3113       return(IntWalk_ArretSurPointPrecedent);
3114   }
3115   //==================================================================================
3116
3117   Standard_Real tolArea = 100.0;
3118   if (ResoU1 < Precision::PConfusion() ||
3119     ResoV1 < Precision::PConfusion() ||
3120     ResoU2 < Precision::PConfusion() ||
3121     ResoV2 < Precision::PConfusion() )
3122     tolArea =  tolArea*2.0;
3123
3124   Standard_Real Cosi1, CosRef1, Ang1, AngRef1, ResoUV1, Duv1, d1, tolCoeff1;   
3125   Standard_Real Cosi2, CosRef2, Ang2, AngRef2, ResoUV2, Duv2, d2, tolCoeff2;   
3126   Cosi1 = Du1*previousd1.X() + Dv1*previousd1.Y();
3127   Cosi2 = Du2*previousd2.X() + Dv2*previousd2.Y();
3128   Duv1 = Du1*Du1 + Dv1*Dv1;
3129   Duv2 = Du2*Du2 + Dv2*Dv2;
3130   ResoUV1 = ResoU1*ResoU1 + ResoV1*ResoV1;
3131   ResoUV2 = ResoU2*ResoU2 + ResoV2*ResoV2;
3132   //
3133   //modified by NIZNHY-PKV Wed Nov 13 12:25:44 2002 f
3134   //
3135   Standard_Real aMinDiv2=Precision::Confusion();
3136   aMinDiv2=aMinDiv2*aMinDiv2;
3137   //
3138   d1=d;
3139   if (Duv1>aMinDiv2)  {
3140     d1 = Abs(ResoUV1/Duv1);
3141     d1 = Min(Sqrt(d1)*tolArea, d);  
3142   } 
3143   //d1 = Abs(ResoUV1/Duv1); 
3144   //d1 = Min(Sqrt(d1)*tolArea,d);  
3145   //modified by NIZNHY-PKV Wed Nov 13 12:34:30 2002 t
3146   tolCoeff1 = Exp(d1);
3147   //
3148   //modified by NIZNHY-PKV Wed Nov 13 12:34:43 2002 f
3149   d2=d;
3150   if (Duv2>aMinDiv2) {
3151     d2 = Abs(ResoUV2/Duv2); 
3152     d2 = Min(Sqrt(d2)*tolArea,d); 
3153   }
3154   //d2 = Abs(ResoUV2/Duv2); 
3155   //d2 = Min(Sqrt(d2)*tolArea,d);  
3156   //modified by NIZNHY-PKV Wed Nov 13 12:34:53 2002 t
3157   tolCoeff2 = Exp(d2);
3158   CosRef1 = CosRef2D/tolCoeff1;
3159   CosRef2 = CosRef2D/tolCoeff2;
3160   //
3161   //==================================================================================
3162   //== The points are not confused :                                           ==
3163   //== 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 ==
3164   //==                           N o t    T o o    G r e a t (angle in space UV)    ==
3165   //==                           C h a n g e    o f    s i d e                ==
3166   //==================================================================================
3167   if (aStatus != IntWalk_PointConfondu) {
3168     if(Cosi1*Cosi1 < CosRef1*Duv1 || Cosi2*Cosi2 < CosRef2*Duv2) {
3169       pasuv[0]*=0.5;  pasuv[1]*=0.5;  pasuv[2]*=0.5;  pasuv[3]*=0.5;
3170       if (pasuv[0]<ResoU1 && pasuv[1]<ResoV1 && pasuv[2]<ResoU2 && pasuv[3]<ResoV2) { 
3171         return(IntWalk_ArretSurPointPrecedent);
3172       }
3173       else {
3174         pasuv[0]*=0.5; pasuv[1]*=0.5; pasuv[2]*=0.5; pasuv[3]*=0.5;
3175         return(IntWalk_PasTropGrand);
3176       }
3177     }
3178     const gp_Dir2d& Tg2dcourante1 = myIntersectionOn2S.DirectionOnS1();
3179     const gp_Dir2d& Tg2dcourante2 = myIntersectionOn2S.DirectionOnS2();
3180     Cosi1 = Du1*Tg2dcourante1.X() + Dv1*Tg2dcourante1.Y();
3181     Cosi2 = Du2*Tg2dcourante2.X() + Dv2*Tg2dcourante2.Y();
3182     Ang1 = Abs(previousd1.Angle(Tg2dcourante1));  
3183     Ang2 = Abs(previousd2.Angle(Tg2dcourante2));  
3184     AngRef1 = AngRef2D*tolCoeff1;
3185     AngRef2 = AngRef2D*tolCoeff2;
3186     //-------------------------------------------------------
3187     //-- Test : Angle too great in space UV       -----
3188     //--        Change of  side                      -----
3189     //-------------------------------------------------------
3190     if(Cosi1*Cosi1 < CosRef1*Duv1 || Cosi2*Cosi2 < CosRef2*Duv2 || Ang1 > AngRef1 || Ang2 > AngRef2) {
3191       pasuv[0]*=0.5;  pasuv[1]*=0.5;  pasuv[2]*=0.5;  pasuv[3]*=0.5;
3192       if (pasuv[0]<ResoU1 && pasuv[1]<ResoV1 && pasuv[2]<ResoU2 && pasuv[3]<ResoV2) 
3193         return(IntWalk_ArretSurPoint);
3194       else 
3195         return(IntWalk_PasTropGrand);
3196     }
3197   }
3198   //<-OCC431(apo)
3199   //==================================================================================
3200   //== D e t e c t i o n   o f    :  Step Too Small 
3201   //==                               STEP TOO Great 
3202   //==================================================================================
3203
3204   //---------------------------------------
3205   //-- Estimate of the vector           --
3206   //---------------------------------------
3207   FlecheCourante =
3208     Sqrt(Abs((previousd.XYZ()-TgCourante.XYZ()).SquareModulus()*aSqDist))/8.;
3209
3210   if ( FlecheCourante<= fleche*0.5) {     //-- Current step too small
3211     if(FlecheCourante>1e-16) { 
3212       Ratio = 0.5*(fleche/FlecheCourante);
3213     }
3214     else { 
3215       Ratio = 10.0;
3216     }
3217     Standard_Real pasSu1 = pasuv[0];
3218     Standard_Real pasSv1 = pasuv[1];
3219     Standard_Real pasSu2 = pasuv[2];
3220     Standard_Real pasSv2 = pasuv[3];
3221
3222     //-- In  case if 
3223     //-- a point at U+DeltaU is required, ....
3224     //-- return a point at U + Epsilon
3225     //-- Epsilon << DeltaU.
3226
3227     if(pasuv[0]< AbsDu1) pasuv[0] = AbsDu1;
3228     if(pasuv[1]< AbsDv1) pasuv[1] = AbsDv1;
3229     if(pasuv[2]< AbsDu2) pasuv[2] = AbsDu2;
3230     if(pasuv[3]< AbsDv2) pasuv[3] = AbsDv2;
3231
3232     if(pasuv[0]<ResoU1) pasuv[0]=ResoU1;
3233     if(pasuv[1]<ResoV1) pasuv[1]=ResoV1;
3234     if(pasuv[2]<ResoU2) pasuv[2]=ResoU2;
3235     if(pasuv[3]<ResoV2) pasuv[3]=ResoV2;
3236     //-- if(Ratio>10.0 ) { Ratio=10.0; } 
3237     Standard_Real R1,R = pasInit[0]/pasuv[0];
3238     R1= pasInit[1]/pasuv[1];     if(R1<R) R=R1;
3239     R1= pasInit[2]/pasuv[2];     if(R1<R) R=R1;
3240     R1= pasInit[3]/pasuv[3];     if(R1<R) R=R1;
3241     if(Ratio > R) Ratio=R;
3242     pasuv[0] = Min(Ratio*pasuv[0],pasInit[0]);
3243     pasuv[1] = Min(Ratio*pasuv[1],pasInit[1]);
3244     pasuv[2] = Min(Ratio*pasuv[2],pasInit[2]);
3245     pasuv[3] = Min(Ratio*pasuv[3],pasInit[3]);
3246     if (pasuv[0] != pasSu1 || pasuv[2] != pasSu2|| 
3247       pasuv[1] != pasSv1 || pasuv[3] != pasSv2) {
3248         if(++STATIC_BLOCAGE_SUR_PAS_TROP_GRAND > 5) {
3249           STATIC_BLOCAGE_SUR_PAS_TROP_GRAND = 0;
3250           return IntWalk_PasTropGrand; 
3251         }
3252     }
3253     if(aStatus == IntWalk_OK) {
3254       STATIC_BLOCAGE_SUR_PAS_TROP_GRAND=0;
3255       //-- Try to increase the step
3256     }
3257     return aStatus;
3258   }
3259   else {                                //-- CurrentVector > vector*0.5 
3260     if (FlecheCourante > fleche) {      //-- Current step too Great
3261       Ratio = fleche/FlecheCourante; 
3262       pasuv[0] = Ratio*pasuv[0];
3263       pasuv[1] = Ratio*pasuv[1];
3264       pasuv[2] = Ratio*pasuv[2];
3265       pasuv[3] = Ratio*pasuv[3];
3266       //if(++STATIC_BLOCAGE_SUR_PAS_TROP_GRAND > 5) {
3267       //        STATIC_BLOCAGE_SUR_PAS_TROP_GRAND = 0;
3268       return IntWalk_PasTropGrand; 
3269       //}
3270     }
3271     else {                             //-- vector/2  <  CurrentVector <= vector   
3272       Ratio = 0.75 * (fleche / FlecheCourante);
3273     }
3274   }
3275
3276   if(aStatus != IntWalk_PointConfondu)
3277   {
3278     //Here, aCosBetweenTangent >= 0.0 definitely.
3279
3280     /*
3281     Brief algorithm description.
3282     We have two (not-coincindent) intersection point (P1 and P2). In every point,
3283     vector of tangent (to the intersection curve) is known (vectors T1 and T2).
3284     Basing on these data, we create osculating circle.
3285
3286                                     * - arc of osculating circle
3287                                 *      * 
3288                            P1 x----------x P2
3289                              /            \
3290                             /              \
3291                           Vec(T1)         Vec(T2)
3292
3293     Let me draw your attention to the following facts:
3294       1. Vectors T1 and T2 direct FROM (not TO) points P1 and P2. Therefore,
3295         one of previously computed vector should be reversed.
3296
3297     In this case, the absolute (!) value of the deflection between the arc of
3298     the osculating circle and the P1P2 segment can be computed as follows:
3299           e = d*(1-sin(B/2))/(2*cos(B/2)),      (1)
3300     where d is the length of P1P2 segment, B is the angle between vectors T1 and T2.
3301     At that,
3302               pi/2 <= B <= pi,
3303               cos(B/2) >= 0,
3304               sin(B/2) > 0,
3305               sin(B) > 0,
3306               cos(B) < 0.
3307
3308     Later, the normal state of algorithm work is (as we apply) 
3309               tolconf/2 <= e <= tolconf.
3310     In this case, we keep previous step.
3311
3312     If e < tolconf/2 then the local curvature of the intersection curve is small.
3313     As result, the step should be increased.
3314
3315     If e > tolconf then the step is too big. Therefore, we should decrease one.
3316
3317     Condition (1) is equivalent to
3318               sin(B/2) = 1 - 2/(1+(d/(2*e))^2) = Fs(e),
3319               cos(B) = 1 - 2*Fs(e)^2 = Fd(e),
3320     where Fs(e)and Fd(e) are some function with parameter "deflection".
3321
3322     Let mean that Fs(e) is decreasing function. Fd(e) is increasing function,
3323     in the range, where Fs(e) > 0.0 (i.e. when e < d/2).
3324
3325     Now, let substitute required deflection (tolconf or tolconf/2) to z. Then
3326     it is necessary to check if e < z or if e > z.
3327
3328     In this case, it is enough to comapare Fs(e) and Fs(z).
3329     At that Fs(e) > 0 because sin(B/2) > 0 always.
3330
3331     Therefore, if Fs(z) < 0.0 then Fs(e) > Fs(z) ==> e < z definitely.
3332     If Fs(z) > 0.0 then we can compare Fs(z)^2 and Fs(e)^2 or, in substance,
3333     values Fd(e) and Fd(z). If Fd(e) > Fd(z) then e > z and vice versa.
3334     */
3335
3336     //Fd(e) is already known (Fd(e) == -aCosBetweenTangent)
3337
3338     const Standard_Real anInvSqAbsArcDeflMax = 0.25*aSqDist/(tolconf*tolconf);
3339     const Standard_Real aSinB2Max = 1.0 - 2.0/(1.0 + anInvSqAbsArcDeflMax);
3340
3341     if(aSinB2Max >= 0.0 && (aCosBetweenTangent <= 2.0 * aSinB2Max * aSinB2Max - 1.0))
3342     {//Real deflection is greater or equal than tolconf
3343       aStatus = IntWalk_PasTropGrand;
3344     }
3345     else
3346     {//Real deflection is less than tolconf
3347       const Standard_Real anInvSqAbsArcDeflMin = 4.0*anInvSqAbsArcDeflMax;
3348       const Standard_Real aSinB2Min = 1.0 - 2.0/(1.0 + anInvSqAbsArcDeflMin);
3349
3350       if((aSinB2Min < 0.0) || (aCosBetweenTangent >= 2.0 * aSinB2Min * aSinB2Min - 1.0))
3351       {//Real deflection is less than tolconf/2.0
3352         aStatus = IntWalk_StepTooSmall;
3353       }
3354     }
3355
3356     if(aStatus == IntWalk_PasTropGrand)
3357     {
3358       pasuv[0]*=0.5; pasuv[1]*=0.5; pasuv[2]*=0.5; pasuv[3]*=0.5;
3359       return aStatus;
3360     }
3361
3362     if(aStatus == IntWalk_StepTooSmall)
3363     {
3364       pasuv[0] = Max(pasuv[0], AbsDu1);
3365       pasuv[1] = Max(pasuv[1], AbsDv1);
3366       pasuv[2] = Max(pasuv[2], AbsDu2);
3367       pasuv[3] = Max(pasuv[3], AbsDv2);
3368
3369       pasInit[0] = Max(pasInit[0], AbsDu1);
3370       pasInit[1] = Max(pasInit[1], AbsDv1);
3371       pasInit[2] = Max(pasInit[2], AbsDu2);
3372       pasInit[3] = Max(pasInit[3], AbsDv2);
3373
3374       return aStatus;
3375     }
3376   }
3377
3378   pasuv[0] = Max(myStepMin[0],Min(Min(Ratio*AbsDu1,pasuv[0]),pasInit[0]));
3379   pasuv[1] = Max(myStepMin[1],Min(Min(Ratio*AbsDv1,pasuv[1]),pasInit[1]));
3380   pasuv[2] = Max(myStepMin[2],Min(Min(Ratio*AbsDu2,pasuv[2]),pasInit[2]));
3381   pasuv[3] = Max(myStepMin[3],Min(Min(Ratio*AbsDv2,pasuv[3]),pasInit[3]));
3382
3383   if(aStatus == IntWalk_OK) STATIC_BLOCAGE_SUR_PAS_TROP_GRAND=0;
3384   return aStatus;
3385 }
3386
3387 Standard_Boolean IntWalk_PWalking::
3388 TestArret(const Standard_Boolean DejaReparti,
3389           TColStd_Array1OfReal& Param,
3390           IntImp_ConstIsoparametric&  ChoixIso)
3391
3392           //
3393           // test if the point of intersection set by these parameters remains in the 
3394           // natural domain of each square.
3395           // if the point outpasses reframe to find the best iso (border)
3396           // that intersects easiest the other square
3397           // otherwise test if closed line is present  
3398           // 
3399 {
3400   Standard_Real Uvd[4],Uvf[4],Epsuv[4],Duv[4],Uvp[4],dv,dv2,ParC[4];
3401   Standard_Real DPc,DPb;
3402   Standard_Integer i = 0, k = 0;
3403   Epsuv[0] = ResoU1;
3404   Epsuv[1] = ResoV1;
3405   Epsuv[2] = ResoU2;
3406   Epsuv[3] = ResoV2;
3407   previousPoint.Parameters(Uvp[0],Uvp[1],Uvp[2],Uvp[3]);
3408
3409   Standard_Real SolParam[4];
3410   myIntersectionOn2S.Point().Parameters(SolParam[0],SolParam[1],SolParam[2],SolParam[3]);
3411
3412   Standard_Boolean Trouve = Standard_False;
3413
3414   Uvd[0]=Um1;   Uvf[0]=UM1;   Uvd[1]=Vm1;   Uvf[1]=VM1;
3415   Uvd[2]=Um2;   Uvf[2]=UM2;   Uvd[3]=Vm2;   Uvf[3]=VM2;
3416
3417   Standard_Integer im1;
3418   for ( i = 1,im1 = 0;i<=4;i++,im1++) {
3419     switch(i) { 
3420     case 1: k=2; break;
3421     case 2: k=1; break;
3422     case 3: k=4; break;
3423     case 4: k=3; break;
3424     }
3425     if (Param(i) < (Uvd[im1]-Epsuv[im1]) ||
3426       SolParam[im1] < (Uvd[im1]-Epsuv[im1]))     //--     Current -----  Bound Inf -----  Previous
3427     {
3428       Trouve    = Standard_True;                   //-- 
3429       DPc       = Uvp[im1]-Param(i);               //--     Previous  - Current
3430       DPb       = Uvp[im1]-Uvd[im1];               //--     Previous  - Bound Inf
3431       ParC[im1] = Uvd[im1];                        //--     ParamCorrige
3432       dv        = Param(k)-Uvp[k-1];               //--     Current   - Previous (other Direction)
3433       dv2       = dv*dv;         
3434       if(dv2>RealEpsilon()) {                       //--    Progress at the other Direction ?
3435         Duv[im1]  = DPc*DPb + dv2;
3436         Duv[im1]  = Duv[im1]*Duv[im1]/(DPc*DPc+dv2)/(DPb*DPb+dv2);
3437       }
3438       else {
3439         Duv[im1]=-1.0;                              //--    If no progress, do not change  
3440       }                                             //--    the choice of iso 
3441     }   
3442     else if (Param(i) > (Uvf[im1] + Epsuv[im1]) ||
3443       SolParam[im1] > (Uvf[im1] + Epsuv[im1]))//--    Previous -----  Bound Sup -----  Current
3444     {
3445       Trouve    = Standard_True;                    //-- 
3446       DPc       = Param(i)-Uvp[im1];                //--     Current   - Previous
3447       DPb       = Uvf[im1]-Uvp[im1];                //--     Bound Sup - Previous 
3448       ParC[im1] = Uvf[im1];                         //--     Param Corrige
3449       dv        = Param(k)-Uvp[k-1];                //--     Current   - Previous (other Direction)
3450       dv2       = dv*dv;
3451       if(dv2>RealEpsilon()) {                       //--     Progress in other Direction ?
3452         Duv[im1]  =  DPc*DPb + dv2;
3453         Duv[im1]  = Duv[im1]*Duv[im1]/(DPc*DPc+dv2)/(DPb*DPb+dv2);
3454       }
3455       else {
3456         Duv[im1]=-1.0;                              //--    If no progress, do not change 
3457       }                                             //--    the choice of iso 
3458     }
3459     else { 
3460       Duv[im1]= -1.;
3461       ParC[im1]=Param(i);
3462     }
3463   }
3464
3465   if (Trouve) {
3466     //--------------------------------------------------
3467     //-- One of Parameters u1,v1,u2,v2 is outside of  --
3468     //-- the natural limits.                          -- 
3469     //-- Find the best direction of                   -- 
3470     //-- progress and reframe the parameters.        --
3471     //--------------------------------------------------
3472     Standard_Real ddv = -1.0;
3473     k=-1;
3474     for (i=0;i<=3;i++) {
3475       Param(i+1) = ParC[i];
3476       if(Duv[i]>ddv) { 
3477         ddv = Duv[i];
3478         k=i;
3479       }
3480     }
3481     if(k!=-1) { 
3482       ChoixIso   = ChoixRef[k];
3483     }
3484     else { 
3485       if((ParC[0]<=Uvd[0]+Epsuv[0]) || (ParC[0]>=Uvf[0]-Epsuv[0])) {
3486         ChoixIso = IntImp_UIsoparametricOnCaro1;
3487       }
3488       else if((ParC[1]<=Uvd[1]+Epsuv[1]) || (ParC[1]>=Uvf[1]-Epsuv[1])) {
3489         ChoixIso = IntImp_VIsoparametricOnCaro1;
3490       }
3491       else if((ParC[2]<=Uvd[2]+Epsuv[2]) || (ParC[2]>=Uvf[2]-Epsuv[2])) {
3492         ChoixIso = IntImp_UIsoparametricOnCaro2;
3493       }
3494       else if((ParC[3]<=Uvd[3]+Epsuv[3]) || (ParC[3]>=Uvf[3]-Epsuv[3])) {
3495         ChoixIso = IntImp_VIsoparametricOnCaro2;
3496       }
3497     }
3498     close = Standard_False;
3499     return Standard_True;
3500   }
3501   else 
3502   {  
3503     if (!DejaReparti) { // find if line closed
3504
3505       Standard_Real u,v;
3506       const IntSurf_PntOn2S& POn2S1=line->Value(1);
3507       //On S1
3508       POn2S1.ParametersOnS1(u,v);
3509       gp_Pnt2d P1uvS1(u,v);
3510       previousPoint.ParametersOnS1(u,v);
3511       gp_Pnt2d PrevuvS1(u,v);
3512       myIntersectionOn2S.Point().ParametersOnS1(u,v);
3513       gp_Pnt2d myIntersuvS1(u,v);
3514       Standard_Boolean close2dS1 = (P1uvS1.XY()-PrevuvS1.XY())*
3515         (P1uvS1.XY()-myIntersuvS1.XY()) < 0.0;
3516       //On S2
3517       POn2S1.ParametersOnS2(u,v);
3518       gp_Pnt2d P1uvS2(u,v);
3519       previousPoint.ParametersOnS2(u,v);
3520       gp_Pnt2d PrevuvS2(u,v);
3521       myIntersectionOn2S.Point().ParametersOnS2(u,v);
3522       gp_Pnt2d myIntersuvS2(u,v);
3523       Standard_Boolean close2dS2 = (P1uvS2.XY()-PrevuvS2.XY())*
3524         (P1uvS2.XY()-myIntersuvS2.XY()) < 0.0;
3525
3526       close = close2dS1 && close2dS2;
3527       return close;
3528     }
3529     else return Standard_False;
3530   }
3531 }
3532