add2c500fcbcf9e591570e4fce6f27f73dd092a9
[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(UM1 - Um1),
762                                   Epsilon(VM1 - Vm1),
763                                   Epsilon(UM2 - Um2),
764                                   Epsilon(VM2 - Vm2)};
765
766   Standard_Integer aPrevNbPoints = line->NbPoints();
767
768   Arrive = Standard_False;
769   while(!Arrive) //010
770   {
771     aPrevStatus = aStatus;
772
773     LevelOfIterWithoutAppend++;
774     if(LevelOfIterWithoutAppend>20)
775     {
776       Arrive = Standard_True; 
777       if(DejaReparti) {
778         break;
779       }
780       RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
781       LevelOfIterWithoutAppend = 0;
782     }
783     //
784     // compute f
785     f = 0.;
786     switch (ChoixIso) { 
787       case IntImp_UIsoparametricOnCaro1: f = Abs(previousd1.X()); break;
788       case IntImp_VIsoparametricOnCaro1: f = Abs(previousd1.Y()); break;
789       case IntImp_UIsoparametricOnCaro2: f = Abs(previousd2.X()); break;
790       case IntImp_VIsoparametricOnCaro2: f = Abs(previousd2.Y()); break;
791       default:break;
792     }
793     //
794     if(f<0.1) {
795       f=0.1;
796     }
797     //
798     previousPoint.Parameters(Param(1),Param(2),Param(3),Param(4));
799     //
800     //--ofv.begin
801     Standard_Real aIncKey, aEps, dP1, dP2, dP3, dP4;
802     //
803     dP1 = sensCheminement * pasuv[0] * previousd1.X() /f;
804     dP2 = sensCheminement * pasuv[1] * previousd1.Y() /f;
805     dP3 = sensCheminement * pasuv[2] * previousd2.X() /f; 
806     dP4 = sensCheminement * pasuv[3] * previousd2.Y() /f;
807     //
808     aIncKey=5.*(Standard_Real)IncKey;
809     aEps=1.e-7;
810     if(ChoixIso == IntImp_UIsoparametricOnCaro1 && Abs(dP1) < aEps)
811     {
812       dP1 *= aIncKey;
813     }
814
815     if(ChoixIso == IntImp_VIsoparametricOnCaro1 && Abs(dP2) < aEps)
816     {
817       dP2 *= aIncKey;
818     }
819
820     if(ChoixIso == IntImp_UIsoparametricOnCaro2 && Abs(dP3) < aEps)
821     {
822       dP3 *= aIncKey;
823     }
824
825     if(ChoixIso == IntImp_VIsoparametricOnCaro2 && Abs(dP4) < aEps)
826     {
827       dP4 *= aIncKey;
828     }
829     //--ofv.end
830     //
831     Param(1) += dP1;
832     Param(2) += dP2;
833     Param(3) += dP3; 
834     Param(4) += dP4;
835     //==========================
836     SvParam[0]=Param(1); 
837     SvParam[1]=Param(2);
838     SvParam[2]=Param(3);
839     SvParam[3]=Param(4);
840     //
841     Standard_Integer aTryNumber = 0;
842     Standard_Boolean isBadPoint = Standard_False;
843     IntImp_ConstIsoparametric aBestIso = ChoixIso;
844     do
845     {
846       isBadPoint = Standard_False;
847
848       ChoixIso= myIntersectionOn2S.Perform(Param, Rsnld, aBestIso);
849
850       if (myIntersectionOn2S.IsDone() && !myIntersectionOn2S.IsEmpty())
851       {
852         //If we go along any surface boundary then it is possible 
853         //to find "outboundaried" point.
854         //Nevertheless, if this deflection is quite small, we will be
855         //able to adjust this point to the boundary.
856
857         Standard_Real aNewPnt[4], anAbsParamDist[4];
858         myIntersectionOn2S.Point().Parameters(aNewPnt[0], aNewPnt[1], aNewPnt[2], aNewPnt[3]);
859         const Standard_Real aParMin[4] = {Um1, Vm1, Um2, Vm2};
860         const Standard_Real aParMax[4] = {UM1, VM1, UM2, VM2};
861
862         for(Standard_Integer i = 0; i < 4; i++)
863         {
864           if(Abs(aNewPnt[i] - aParMin[i]) < aTol[i])
865             aNewPnt[i] = aParMin[i];
866           else if(Abs(aNewPnt[i] - aParMax[i]) < aTol[i])
867             aNewPnt[i] = aParMax[i];
868         }
869
870         if (aNewPnt[0] < Um1 || aNewPnt[0] > UM1 ||
871             aNewPnt[1] < Vm1 || aNewPnt[1] > VM1 ||
872             aNewPnt[2] < Um2 || aNewPnt[2] > UM2 ||
873             aNewPnt[3] < Vm2 || aNewPnt[3] > VM2)
874         {
875           // Out of borders, process it later.
876           break; 
877         }
878
879         myIntersectionOn2S.ChangePoint().SetValue(aNewPnt[0],
880                                                   aNewPnt[1],
881                                                   aNewPnt[2],
882                                                   aNewPnt[3]);
883
884         anAbsParamDist[0] = Abs(Param(1) - dP1 - aNewPnt[0]);
885         anAbsParamDist[1] = Abs(Param(2) - dP2 - aNewPnt[1]);
886         anAbsParamDist[2] = Abs(Param(3) - dP3 - aNewPnt[2]);
887         anAbsParamDist[3] = Abs(Param(4) - dP4 - aNewPnt[3]);
888         if (anAbsParamDist[0] < ResoU1 &&
889             anAbsParamDist[1] < ResoV1 &&
890             anAbsParamDist[2] < ResoU2 &&
891             anAbsParamDist[3] < ResoV2 &&
892             aStatus != IntWalk_PasTropGrand)
893         {
894           isBadPoint = Standard_True;
895           aBestIso = IntImp_ConstIsoparametric((aBestIso + 1) % 4);
896         }
897       }
898     } while (isBadPoint && ++aTryNumber <= 4);
899     //
900     if (!myIntersectionOn2S.IsDone())
901     {
902       //end of line, division
903       Arrive = Standard_False;
904       Param(1)=SvParam[0]; 
905       Param(2)=SvParam[1]; 
906       Param(3)=SvParam[2];
907       Param(4)=SvParam[3];
908       RepartirOuDiviser(DejaReparti, ChoixIso, Arrive);
909     }
910     else  //009 
911     {
912       //== Calculation of exact point from Param(.) is possible
913       if (myIntersectionOn2S.IsEmpty())
914       {
915         Standard_Real u1,v1,u2,v2;
916         previousPoint.Parameters(u1,v1,u2,v2);
917         //
918         Arrive = Standard_False;
919         if(u1<UFirst1 || u1>ULast1)
920         {
921           Arrive=Standard_True;
922         }       
923
924         if(u2<UFirst2 || u2>ULast2)
925         {
926           Arrive=Standard_True;
927         }
928
929         if(v1<VFirst1 || v1>VLast1)
930         {
931           Arrive=Standard_True;
932         }
933
934         if(v2<VFirst2 || v2>VLast2)
935         {
936           Arrive=Standard_True;
937         }
938
939         RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
940         LevelOfEmptyInmyIntersectionOn2S++;
941         //
942         if(LevelOfEmptyInmyIntersectionOn2S>10)
943         {
944           pasuv[0]=pasSav[0]; 
945           pasuv[1]=pasSav[1]; 
946           pasuv[2]=pasSav[2]; 
947           pasuv[3]=pasSav[3];
948         }
949       }
950       else //008
951       {
952         //============================================================
953         //== A point has been found :  T E S T   D E F L E C T I O N 
954         //============================================================
955         if(NoTestDeflection)
956         {
957           NoTestDeflection = Standard_False;
958         }                 
959         else
960         {
961           if(--LevelOfEmptyInmyIntersectionOn2S<=0)
962           {
963             LevelOfEmptyInmyIntersectionOn2S=0;
964             if(LevelOfIterWithoutAppend < 10)
965             {
966               aStatus = TestDeflection(ChoixIso, aStatus);
967             }                   
968             else
969             {
970               if (aStatus != IntWalk_StepTooSmall)
971               {
972                 // Bug #0029682 (Boolean intersection with
973                 // fuzzy-option hangs). 
974                 // If aStatus == IntWalk_StepTooSmall then
975                 // the counter "LevelOfIterWithoutAppend" will
976                 // be nulified in the future if the step is smaller
977                 // (see "case IntWalk_StepTooSmall:" below).
978                 // Here, we forbid this nulification and thereby provide out from
979                 // the algorithm.
980                 pasuv[0] *= 0.5;
981                 pasuv[1] *= 0.5;
982                 pasuv[2] *= 0.5;
983                 pasuv[3] *= 0.5;
984               }
985             }
986           }
987         }
988
989         //============================================================
990         //==       T r a i t e m e n t   s u r   S t a t u s        ==
991         //============================================================
992         if(LevelOfPointConfondu > 5)
993         { 
994           aStatus = IntWalk_ArretSurPoint;
995           LevelOfPointConfondu = 0;  
996         }
997         //
998         if(aStatus==IntWalk_OK)
999         { 
1000           NbPasOKConseq++;
1001           if(NbPasOKConseq >= 5)
1002           {
1003             NbPasOKConseq=0;
1004             Standard_Boolean pastroppetit;
1005             Standard_Real t;
1006             //
1007             do
1008             {
1009               pastroppetit=Standard_True;
1010               //
1011               if(pasuv[0]<pasInit[0])
1012               {
1013                 t = (pasInit[0]-pasuv[0])*0.25;
1014                 if(t>0.1*pasInit[0])
1015                 {
1016                   t=0.1*pasuv[0];
1017                 }
1018
1019                 pasuv[0]+=t; 
1020                 pastroppetit=Standard_False;
1021               }
1022
1023               if(pasuv[1]<pasInit[1])
1024               {
1025                 t = (pasInit[1]-pasuv[1])*0.25;
1026                 if(t>0.1*pasInit[1]) {
1027                   t=0.1*pasuv[1];
1028                 }               
1029
1030                 pasuv[1]+=t; 
1031                 pastroppetit=Standard_False;
1032               }
1033
1034               if(pasuv[2]<pasInit[2])
1035               {
1036                 t = (pasInit[2]-pasuv[2])*0.25;
1037                 if(t>0.1*pasInit[2])
1038                 {
1039                   t=0.1*pasuv[2];
1040                 }
1041
1042                 pasuv[2]+=t; 
1043                 pastroppetit=Standard_False;
1044               }
1045
1046               if(pasuv[3]<pasInit[3])
1047               {
1048                 t = (pasInit[3]-pasuv[3])*0.25;
1049                 if(t>0.1*pasInit[3]) {
1050                   t=0.1*pasuv[3];
1051                 }
1052                 pasuv[3]+=t; 
1053                 pastroppetit=Standard_False;
1054               }
1055               if(pastroppetit)
1056               {
1057                 if(pasMax<0.1)
1058                 {
1059                   pasMax*=1.1;
1060                   pasInit[0]*=1.1; 
1061                   pasInit[1]*=1.1; 
1062                   pasInit[2]*=1.1; 
1063                   pasInit[3]*=1.1; 
1064                 }
1065                 else
1066                 {
1067                   pastroppetit=Standard_False;
1068                 }
1069               }
1070             }
1071             while(pastroppetit);
1072           }
1073         }//aStatus==IntWalk_OK
1074         else
1075           NbPasOKConseq=0;
1076
1077         //
1078         switch(aStatus)//007
1079         {
1080         case IntWalk_ArretSurPointPrecedent:
1081           {
1082             Arrive = Standard_False;
1083             RepartirOuDiviser(DejaReparti, ChoixIso, Arrive);
1084             break;
1085           }
1086         case IntWalk_PasTropGrand:
1087           {
1088             Param(1)=SvParam[0];
1089             Param(2)=SvParam[1]; 
1090             Param(3)=SvParam[2]; 
1091             Param(4)=SvParam[3];
1092
1093             if(LevelOfIterWithoutAppend > 5)
1094             {
1095               for (Standard_Integer i = 0; i < 4; i++)
1096               {
1097                 if (pasSav[i] > pasInit[i])
1098                   continue;
1099
1100                 const Standard_Real aDelta = (pasInit[i]-pasSav[i])*0.25;
1101
1102                 if(aDelta > Epsilon(pasInit[i]))
1103                 {
1104                   pasInit[i] -= aDelta;
1105                   if ((aPrevStatus != IntWalk_StepTooSmall) &&
1106                       (line->NbPoints() != aPrevNbPoints))
1107                   {
1108                     LevelOfIterWithoutAppend = 0;
1109                   }
1110
1111                   aPrevNbPoints = line->NbPoints();
1112                 }
1113               }
1114             }
1115
1116             break;
1117           }
1118         case IntWalk_PointConfondu:
1119           {
1120             LevelOfPointConfondu++;
1121
1122             if(LevelOfPointConfondu>5)
1123             {
1124               Standard_Boolean pastroppetit;
1125               //
1126               do
1127               {
1128                 pastroppetit=Standard_True;
1129
1130                 for(Standard_Integer i = 0; i < 4; i++)
1131                 {
1132                   if(pasuv[i]<pasInit[i])
1133                   {
1134                     pasuv[i]+=(pasInit[i]-pasuv[i])*0.25;
1135                     pastroppetit=Standard_False;
1136                   }
1137                 }
1138
1139                 if(pastroppetit)
1140                 {
1141                   if(pasMax<0.1)
1142                   {
1143                     pasMax*=1.1;
1144                     pasInit[0]*=1.1;
1145                     pasInit[1]*=1.1;
1146                     pasInit[2]*=1.1;
1147                     pasInit[3]*=1.1; 
1148                   }
1149                   else
1150                   {
1151                     pastroppetit=Standard_False;
1152                   }
1153                 }
1154               }
1155               while(pastroppetit);
1156             }
1157
1158             break;
1159           }
1160         case IntWalk_StepTooSmall:
1161           {
1162             Standard_Boolean hasStepBeenIncreased = Standard_False;
1163
1164             for(Standard_Integer i = 0; i < 4; i++)
1165             {
1166               const Standard_Real aNewStep = Min(1.5*pasuv[i], pasInit[i]);
1167               if(aNewStep > pasuv[i])
1168               {
1169                 pasuv[i] = aNewStep;
1170                 hasStepBeenIncreased = Standard_True;
1171               }
1172             }
1173
1174             if(hasStepBeenIncreased)
1175             {
1176               Param(1)=SvParam[0];
1177               Param(2)=SvParam[1];
1178               Param(3)=SvParam[2];
1179               Param(4)=SvParam[3];
1180
1181               // In order to avoid cyclic changes
1182               // (PasTropGrand --> Decrease step --> 
1183               // StepTooSmall --> Increase step --> PasTropGrand...)
1184               // nullify LevelOfIterWithoutAppend only if the condition
1185               // is satisfied:
1186               if ((aPrevStatus != IntWalk_PasTropGrand) &&
1187                   (line->NbPoints() != aPrevNbPoints))
1188               {
1189                 LevelOfIterWithoutAppend = 0;
1190               }
1191
1192               aPrevNbPoints = line->NbPoints();
1193
1194               break;
1195             }
1196           }
1197           Standard_FALLTHROUGH
1198         case IntWalk_OK:
1199         case IntWalk_ArretSurPoint://006
1200           {
1201             //=======================================================
1202             //== Stop Test t   :  Frame on Param(.)     ==
1203             //=======================================================
1204             //xft arrive here
1205             Arrive = TestArret(DejaReparti,Param,ChoixIso); 
1206             // JMB 30th December 1999. 
1207             // Some statement below should not be put in comment because they are useful.
1208             // See grid CTO 909 A1 which infinitely loops 
1209             if(Arrive==Standard_False && aStatus==IntWalk_ArretSurPoint)
1210             {
1211               Arrive=Standard_True;
1212 #ifdef OCCT_DEBUG
1213               std::cout << "IntWalk_PWalking_1.gxx: Problems with intersection"<<std::endl;
1214 #endif
1215             }
1216
1217             if(Arrive)
1218             {
1219               NbPasOKConseq = -10;
1220             }
1221
1222             if(!Arrive)//005
1223             {
1224               //=====================================================
1225               //== Param(.) is in the limits                       ==
1226               //==  and does not end a closed  line                ==
1227               //=====================================================
1228               //== Check on the current point of myInters
1229               Standard_Boolean pointisvalid = Standard_False;
1230               {
1231                 Standard_Real u1,v1,u2,v2; 
1232                 myIntersectionOn2S.Point().Parameters(u1,v1,u2,v2);
1233
1234                 //
1235                 if(u1 <= UM1  && u2 <= UM2 && v1 <= VM1 && 
1236                   v2 <= VM2  && u1 >= Um1 && u2 >= Um2 &&
1237                   v1 >= Vm1  && v2 >= Vm2)
1238                 {
1239                   pointisvalid=Standard_True;
1240                 }
1241               }
1242
1243               //
1244               if(pointisvalid)
1245               {
1246                 previousPoint = myIntersectionOn2S.Point();
1247                 previoustg = myIntersectionOn2S.IsTangent();
1248
1249                 if(!previoustg)
1250                 {
1251                   previousd  = myIntersectionOn2S.Direction();
1252                   previousd1 = myIntersectionOn2S.DirectionOnS1();
1253                   previousd2 = myIntersectionOn2S.DirectionOnS2();
1254                 }
1255                 //=====================================================
1256                 //== Check on the previous Point
1257                 {
1258                   Standard_Real u1,v1,u2,v2;
1259                   previousPoint.Parameters(u1,v1,u2,v2); 
1260                   if( u1 <= UM1  && u2 <= UM2 && v1 <= VM1 &&
1261                     v2 <= VM2  && u1 >= Um1 && u2 >= Um2 &&
1262                     v1 >= Vm1  && v2 >= Vm2)
1263                   {
1264                     pl = previousPoint.Value();
1265                     if(bTestFirstPoint)
1266                     {
1267                       if(pf.SquareDistance(pl) < aSQDistMax)
1268                       {
1269                         IncKey++;
1270                         if(IncKey == 5000)
1271                           return;
1272                         else
1273                           continue;
1274                       }
1275                       else
1276                       {
1277                         bTestFirstPoint = Standard_False;
1278                       }
1279                     }
1280                     //
1281                     AddAPoint(previousPoint);
1282                     RejectIndex++;
1283
1284                     if(RejectIndex >= RejectIndexMAX)
1285                     {
1286                       Arrive = Standard_True;
1287                       break;
1288                     }
1289
1290                     //
1291                     LevelOfIterWithoutAppend = 0;
1292                   }
1293                 }
1294               }//pointisvalid
1295               //====================================================
1296
1297               if (aStatus == IntWalk_ArretSurPoint)
1298               {
1299                 RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
1300               }
1301               else
1302               {
1303                 if (line->NbPoints() == 2)
1304                 {
1305                   pasSav[0] = pasuv[0];
1306                   pasSav[1] = pasuv[1];
1307                   pasSav[2] = pasuv[2];
1308                   pasSav[3] = pasuv[3];
1309                   
1310                   if ((aPrevStatus == IntWalk_PasTropGrand) &&
1311                       (LevelOfIterWithoutAppend > 0))
1312                   {
1313                     pasInit[0] = pasuv[0];
1314                     pasInit[1] = pasuv[1];
1315                     pasInit[2] = pasuv[2];
1316                     pasInit[3] = pasuv[3];
1317                   }
1318                 }
1319               }
1320             }//005 if(!Arrive)
1321             else  //004
1322             {
1323               if(close)
1324               {
1325                 //================= la ligne est fermee ===============
1326                 AddAPoint(line->Value(1)); //ligne fermee
1327                 LevelOfIterWithoutAppend=0;
1328               }
1329               else    //$$$
1330               {
1331                 //====================================================
1332                 //== Param was not in the limits (was reframed)
1333                 //====================================================
1334                 Standard_Boolean bPrevNotTangent = !previoustg || !myIntersectionOn2S.IsTangent();
1335
1336                 IntImp_ConstIsoparametric SauvChoixIso = ChoixIso;
1337                 ChoixIso = myIntersectionOn2S.Perform(Param,Rsnld,ChoixIso);
1338                 //
1339                 if(!myIntersectionOn2S.IsEmpty()) //002
1340                 {
1341                   // mutially outpasses in the square or intersection in corner
1342
1343                   if(TestArret(Standard_True,Param,ChoixIso))
1344                   {
1345                     NbPasOKConseq = -10;
1346                     ChoixIso = myIntersectionOn2S.Perform(Param,Rsnld,ChoixIso);
1347
1348                     if(!myIntersectionOn2S.IsEmpty())
1349                     {
1350                       previousPoint = myIntersectionOn2S.Point();
1351                       previoustg = myIntersectionOn2S.IsTangent();
1352
1353                       if (!previoustg)
1354                       {
1355                         previousd  = myIntersectionOn2S.Direction();
1356                         previousd1 = myIntersectionOn2S.DirectionOnS1();
1357                         previousd2 = myIntersectionOn2S.DirectionOnS2();
1358                       }
1359
1360                       pl = previousPoint.Value();
1361
1362                       if(bTestFirstPoint)
1363                       {
1364                         if(pf.SquareDistance(pl) < aSQDistMax)
1365                         {
1366                           IncKey++;
1367                           if(IncKey == 5000)
1368                             return;
1369                           else
1370                             continue;
1371                         }
1372                         else
1373                         {
1374                           bTestFirstPoint = Standard_False;
1375                         }
1376                       }
1377                       //
1378                       AddAPoint(previousPoint);
1379                       RejectIndex++;
1380
1381                       if(RejectIndex >= RejectIndexMAX)
1382                       {
1383                         Arrive = Standard_True;
1384                         break;
1385                       }
1386
1387                       //
1388                       LevelOfIterWithoutAppend=0;
1389                       RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
1390                     }
1391                     else
1392                     {
1393                       //fail framing divides the step
1394                       Arrive = Standard_False;
1395                       RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
1396                       NoTestDeflection = Standard_True;
1397                       ChoixIso = SauvChoixIso;
1398                     }
1399                   }//if(TestArret())
1400                   else
1401                   {
1402                     // save the last point
1403                     // to revert to it if the current point is out of bounds
1404
1405                     IntSurf_PntOn2S previousPointSave = previousPoint;
1406                     Standard_Boolean previoustgSave   = previoustg;
1407                     gp_Dir previousdSave              = previousd;
1408                     gp_Dir2d previousd1Save           = previousd1;
1409                     gp_Dir2d previousd2Save           = previousd2;
1410
1411                     previousPoint = myIntersectionOn2S.Point();
1412                     previoustg = myIntersectionOn2S.IsTangent();
1413                     Arrive = Standard_False;
1414
1415                     if(!previoustg)
1416                     {
1417                       previousd  = myIntersectionOn2S.Direction();
1418                       previousd1 = myIntersectionOn2S.DirectionOnS1();
1419                       previousd2 = myIntersectionOn2S.DirectionOnS2();
1420                     }
1421
1422                     //========================================
1423                     //== Check on PreviousPoint @@
1424
1425                     {
1426                       Standard_Real u1,v1,u2,v2;
1427                       previousPoint.Parameters(u1,v1,u2,v2);
1428
1429                       //To save initial 2d points
1430                       gp_Pnt2d ParamPntOnS1(Param(1), Param(2));
1431                       gp_Pnt2d ParamPntOnS2(Param(3), Param(4));
1432
1433                       ///////////////////////////
1434                       Param(1) = u1;
1435                       Param(2) = v1;
1436                       Param(3) = u2;
1437                       Param(4) = v2;
1438                       //
1439
1440                       //xf
1441                       Standard_Boolean bFlag1, bFlag2;
1442                       Standard_Real aTol2D=1.e-11;
1443                       //
1444                       bFlag1=u1 >= Um1-aTol2D && v1 >= Vm1-aTol2D && u1 <= UM1+aTol2D && v1 <= VM1+aTol2D;
1445                       bFlag2=u2 >= Um2-aTol2D && v2 >= Vm2-aTol2D && u2 <= UM2+aTol2D && v2 <= VM2+aTol2D;
1446                       if (bFlag1 && bFlag2)
1447                       {
1448                         if (line->NbPoints() > 1)
1449                         {
1450                           IntSurf_PntOn2S prevprevPoint = line->Value(line->NbPoints()-1);
1451                           Standard_Real ppU1, ppV1, ppU2, ppV2;
1452                           prevprevPoint.Parameters(ppU1, ppV1, ppU2, ppV2);
1453                           Standard_Real pU1, pV1, pU2, pV2;
1454                           previousPointSave.Parameters(pU1, pV1, pU2, pV2);
1455                           gp_Vec2d V1onS1(gp_Pnt2d(ppU1, ppV1), gp_Pnt2d(pU1, pV1));
1456                           gp_Vec2d V2onS1(gp_Pnt2d(pU1, pV1), gp_Pnt2d(u1, v1));
1457                           gp_Vec2d V1onS2(gp_Pnt2d(ppU2, ppV2), gp_Pnt2d(pU2, pV2));
1458                           gp_Vec2d V2onS2(gp_Pnt2d(pU2, pV2), gp_Pnt2d(u2, v2));
1459
1460                           const Standard_Real aDot1 = V1onS1 * V2onS1;
1461                           const Standard_Real aDot2 = V1onS2 * V2onS2;
1462
1463                           if ((aDot1 < 0.0) || (aDot2 < 0.0))
1464                           {
1465                             Arrive = Standard_True;
1466                             break;
1467                           }
1468                         }
1469                         /*
1470                         if(u1 <= UM1  && u2 <= UM2 && v1 <= VM1 &&
1471                         v2 <= VM2  && u1 >= Um1 && u2 >= Um2 &&
1472                         v1 >= Vm1  && v2 >= Vm2)  {
1473                         */                      
1474                         //xt
1475                         pl = previousPoint.Value();
1476
1477                         if(bTestFirstPoint)
1478                         {
1479                           if(pf.SquareDistance(pl) < aSQDistMax)
1480                           {
1481                             IncKey++;
1482
1483                             if(IncKey == 5000)
1484                               return;
1485                             else
1486                               continue;
1487                           }
1488                           else
1489                           {
1490                             bTestFirstPoint = Standard_False;
1491                           }
1492                         }
1493
1494                         //To avoid walking around the same point
1495                         //in the tangent zone near a border
1496
1497                         if (previoustg)
1498                         {
1499                           //There are three consecutive points:
1500                           //previousPointSave -> ParamPnt -> curPnt.
1501
1502                           Standard_Real prevU1, prevV1, prevU2, prevV2;
1503                           previousPointSave.Parameters(prevU1, prevV1, prevU2, prevV2);
1504                           gp_Pnt2d prevPntOnS1(prevU1, prevV1), prevPntOnS2(prevU2, prevV2);
1505                           gp_Pnt2d curPntOnS1(u1, v1), curPntOnS2(u2, v2);
1506                           gp_Vec2d PrevToParamOnS1(prevPntOnS1, ParamPntOnS1);
1507                           gp_Vec2d PrevToCurOnS1(prevPntOnS1, curPntOnS1);
1508                           gp_Vec2d PrevToParamOnS2(prevPntOnS2, ParamPntOnS2);
1509                           gp_Vec2d PrevToCurOnS2(prevPntOnS2, curPntOnS2);
1510                           Standard_Real MaxAngle = 3*M_PI/4;
1511                           Standard_Real anAngleS1 = 0.0, anAngleS2 = 0.0;
1512                           const Standard_Real aSQMParS1 = PrevToParamOnS1.SquareMagnitude();
1513                           const Standard_Real aSQMParS2 = PrevToParamOnS2.SquareMagnitude();
1514                           const Standard_Real aSQMCurS1 = PrevToCurOnS1.SquareMagnitude();
1515                           const Standard_Real aSQMCurS2 = PrevToCurOnS2.SquareMagnitude();
1516
1517                           if(aSQMCurS1 < gp::Resolution())
1518                           {
1519                             //We came back to the one of previos point.
1520                             //Therefore, we must break;
1521
1522                             anAngleS1 = M_PI;
1523                           }
1524                           else if(aSQMParS1 < gp::Resolution())
1525                           {
1526                             //We are walking along tangent zone.
1527                             //It should be continued.
1528                             anAngleS1 = 0.0;
1529                           }
1530                           else
1531                           {
1532                             anAngleS1 = Abs(PrevToParamOnS1.Angle(PrevToCurOnS1));
1533                           }
1534
1535                           if(aSQMCurS2 < gp::Resolution())
1536                           {
1537                             //We came back to the one of previos point.
1538                             //Therefore, we must break;
1539
1540                             anAngleS2 = M_PI;
1541                           }
1542                           else if(aSQMParS2 < gp::Resolution())
1543                           {
1544                             //We are walking along tangent zone.
1545                             //It should be continued;
1546                             anAngleS2 = 0.0;
1547                           }
1548                           else
1549                           {
1550                             anAngleS2 = Abs(PrevToParamOnS2.Angle(PrevToCurOnS2));
1551                           }
1552
1553                           if ((anAngleS1 > MaxAngle) && (anAngleS2 > MaxAngle))
1554                           {
1555                             Arrive = Standard_True;
1556                             break;
1557                           }
1558
1559                           {
1560                             //Check singularity.
1561                             //I.e. check if we are walking along direction, which does not
1562                             //result in comming to any point (i.e. derivative 
1563                             //3D-intersection curve along this direction is equal to 0).
1564                             //A sphere with direction {dU=1, dV=0} from point
1565                             //(U=0, V=M_PI/2) can be considered as example for
1566                             //this case (we cannot find another 3D-point if we go thus).
1567
1568                             //Direction chosen along 1st and 2nd surface correspondingly
1569                             const gp_Vec2d  aDirS1(prevPntOnS1, curPntOnS1),
1570                                             aDirS2(prevPntOnS2, curPntOnS2);
1571
1572                             gp_Pnt aPtemp;
1573                             gp_Vec aDuS1, aDvS1, aDuS2, aDvS2;
1574
1575                             myIntersectionOn2S.Function().AuxillarSurface1()->
1576                                   D1(curPntOnS1.X(), curPntOnS1.Y(), aPtemp, aDuS1, aDvS1);
1577                             myIntersectionOn2S.Function().AuxillarSurface2()->
1578                                   D1(curPntOnS2.X(), curPntOnS2.Y(), aPtemp, aDuS2, aDvS2);
1579
1580                             //Derivative WLine along (it is vector-function indeed)
1581                             //directions chosen
1582                             //(https://en.wikipedia.org/wiki/Directional_derivative#Variation_using_only_direction_of_vector).
1583                             //F1 - on the 1st surface, F2 - on the 2nd surface.
1584                             //x, y, z - coordinates of derivative vector.
1585                             const Standard_Real aF1x =  aDuS1.X()*aDirS1.X() + 
1586                                                         aDvS1.X()*aDirS1.Y();
1587                             const Standard_Real aF1y =  aDuS1.Y()*aDirS1.X() +
1588                                                         aDvS1.Y()*aDirS1.Y();
1589                             const Standard_Real aF1z =  aDuS1.Z()*aDirS1.X() +
1590                                                         aDvS1.Z()*aDirS1.Y();
1591                             const Standard_Real aF2x =  aDuS2.X()*aDirS2.X() +
1592                                                         aDvS2.X()*aDirS2.Y();
1593                             const Standard_Real aF2y =  aDuS2.Y()*aDirS2.X() +
1594                                                         aDvS2.Y()*aDirS2.Y();
1595                             const Standard_Real aF2z =  aDuS2.Z()*aDirS2.X() +
1596                                                         aDvS2.Z()*aDirS2.Y();
1597
1598                             const Standard_Real aF1 = aF1x*aF1x + aF1y*aF1y + aF1z*aF1z;
1599                             const Standard_Real aF2 = aF2x*aF2x + aF2y*aF2y + aF2z*aF2z;
1600
1601                             if((aF1 < gp::Resolution()) && (aF2 < gp::Resolution()))
1602                             {
1603                               //All derivative are equal to 0. Therefore, there is
1604                               //no point in going along direction chosen.
1605                               Arrive = Standard_True;
1606                               break;
1607                             }
1608                           }
1609                         }//if (previoustg) cond.
1610
1611                         ////////////////////////////////////////
1612                         AddAPoint(previousPoint);
1613                         RejectIndex++;
1614
1615                         if(RejectIndex >= RejectIndexMAX)
1616                         {
1617                           Arrive = Standard_True;
1618                           break;
1619                         }
1620
1621                         //
1622
1623                         LevelOfIterWithoutAppend=0;
1624                         Arrive = Standard_True;
1625                       }
1626                       else
1627                       {
1628                         // revert to the last correctly calculated point
1629                         previousPoint = previousPointSave;
1630                         previoustg    = previoustgSave;
1631                         previousd     = previousdSave;
1632                         previousd1    = previousd1Save;
1633                         previousd2    = previousd2Save;
1634                       }
1635                     }
1636
1637                     //
1638                     Standard_Boolean wasExtended = Standard_False;
1639
1640                     if(Arrive && myIntersectionOn2S.IsTangent() && bPrevNotTangent)
1641                     {
1642                       if(ExtendLineInCommonZone(SauvChoixIso, DejaReparti))
1643                       {
1644                         wasExtended = Standard_True;
1645                         Arrive = Standard_False;
1646                         ChoixIso = SauvChoixIso;
1647                       }
1648                     }
1649
1650                     RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
1651
1652                     if(Arrive && 
1653                       myIntersectionOn2S.IsDone() && !myIntersectionOn2S.IsEmpty() &&
1654                       myIntersectionOn2S.IsTangent() && bPrevNotTangent &&
1655                       !wasExtended)
1656                     {
1657                       if(ExtendLineInCommonZone(SauvChoixIso, DejaReparti))
1658                       {
1659                         wasExtended = Standard_True;
1660                         Arrive = Standard_False;
1661                         ChoixIso = SauvChoixIso;
1662                       }
1663                     }
1664                   }//else !TestArret() $
1665                 }//$$ end successful framing on border (!myIntersectionOn2S.IsEmpty())
1666                 else
1667                 {
1668                   //echec framing on border; division of step 
1669                   Arrive = Standard_False;
1670                   NoTestDeflection = Standard_True;
1671                   RepartirOuDiviser(DejaReparti,ChoixIso,Arrive);
1672                 }
1673               }//$$$ end framing on border (!close)
1674             }//004 fin TestArret return Arrive = True
1675           } // 006case IntWalk_ArretSurPoint:  end Processing Status = OK  or ArretSurPoint 
1676         } //007  switch(aStatus)
1677       } //008 end processing point  (TEST DEFLECTION)
1678     } //009 end processing line (else if myIntersectionOn2S.IsDone())
1679   }  //010 end if first departure point allows marching  while (!Arrive)
1680
1681   done = Standard_True;
1682 }
1683 // ===========================================================================================================
1684 // function: ExtendLineInCommonZone
1685 // purpose:  Extends already computed line inside tangent zone in the direction given by theChoixIso.
1686 //           Returns Standard_True if the line was extended through tangent zone and the last computed point 
1687 //           is outside the tangent zone (but it is not put into the line). Otherwise returns Standard_False.
1688 // ===========================================================================================================
1689 Standard_Boolean IntWalk_PWalking::ExtendLineInCommonZone(const IntImp_ConstIsoparametric theChoixIso,
1690                                                           const Standard_Boolean          theDirectionFlag) 
1691 {
1692   Standard_Boolean bOutOfTangentZone = Standard_False;
1693   Standard_Boolean bStop = !myIntersectionOn2S.IsTangent();
1694   Standard_Integer dIncKey = 1;
1695   TColStd_Array1OfReal Param(1,4);
1696   IntWalk_StatusDeflection aStatus = IntWalk_OK;
1697   Standard_Integer nbIterWithoutAppend = 0;
1698   Standard_Integer nbEqualPoints = 0;
1699   Standard_Integer parit = 0;
1700   Standard_Integer uvit = 0;
1701   IntSurf_SequenceOfPntOn2S aSeqOfNewPoint;
1702
1703   while (!bStop) {
1704     nbIterWithoutAppend++;
1705
1706     if((nbIterWithoutAppend > 20) || (nbEqualPoints > 20)) {
1707 #ifdef OCCT_DEBUG
1708       std::cout<<"Infinite loop detected. Stop iterations (IntWalk_PWalking_1.gxx)" << std::endl;
1709 #endif
1710       bStop = Standard_True;
1711       break;
1712     }
1713     Standard_Real f = 0.;
1714
1715     switch (theChoixIso)
1716     { 
1717     case IntImp_UIsoparametricOnCaro1: f = Abs(previousd1.X()); break;
1718     case IntImp_VIsoparametricOnCaro1: f = Abs(previousd1.Y()); break;
1719     case IntImp_UIsoparametricOnCaro2: f = Abs(previousd2.X()); break;
1720     case IntImp_VIsoparametricOnCaro2: f = Abs(previousd2.Y()); break;
1721     }
1722
1723     if(f<0.1) f=0.1;
1724
1725     previousPoint.Parameters(Param(1),Param(2),Param(3),Param(4));
1726
1727     Standard_Real dP1 = sensCheminement * pasuv[0] * previousd1.X() /f;
1728     Standard_Real dP2 = sensCheminement * pasuv[1] * previousd1.Y() /f;
1729     Standard_Real dP3 = sensCheminement * pasuv[2] * previousd2.X() /f; 
1730     Standard_Real dP4 = sensCheminement * pasuv[3] * previousd2.Y() /f;
1731
1732     if(theChoixIso == IntImp_UIsoparametricOnCaro1 && Abs(dP1) < 1.e-7) dP1 *= (5. * (Standard_Real)dIncKey);
1733     if(theChoixIso == IntImp_VIsoparametricOnCaro1 && Abs(dP2) < 1.e-7) dP2 *= (5. * (Standard_Real)dIncKey);
1734     if(theChoixIso == IntImp_UIsoparametricOnCaro2 && Abs(dP3) < 1.e-7) dP3 *= (5. * (Standard_Real)dIncKey);
1735     if(theChoixIso == IntImp_VIsoparametricOnCaro2 && Abs(dP4) < 1.e-7) dP4 *= (5. * (Standard_Real)dIncKey);
1736
1737     Param(1) += dP1;
1738     Param(2) += dP2;
1739     Param(3) += dP3; 
1740     Param(4) += dP4;
1741     Standard_Real SvParam[4];
1742     IntImp_ConstIsoparametric ChoixIso = theChoixIso;
1743
1744     for(parit = 0; parit < 4; parit++) {
1745       SvParam[parit] = Param(parit+1);
1746     }
1747     math_FunctionSetRoot  Rsnld(myIntersectionOn2S.Function());
1748     ChoixIso = myIntersectionOn2S.Perform(Param,Rsnld, theChoixIso);
1749
1750     if (!myIntersectionOn2S.IsDone()) {
1751       return bOutOfTangentZone;
1752     }
1753     else {
1754       if (myIntersectionOn2S.IsEmpty()) {
1755         return bOutOfTangentZone;
1756       }
1757
1758       aStatus = TestDeflection(ChoixIso, aStatus);
1759
1760       if(aStatus == IntWalk_OK) {
1761
1762         for(uvit = 0; uvit < 4; uvit++) {
1763           if(pasuv[uvit] < pasInit[uvit]) {
1764             pasuv[uvit] = pasInit[uvit];
1765           }
1766         }
1767       }
1768
1769       switch(aStatus) {
1770       case  IntWalk_ArretSurPointPrecedent:
1771         {
1772           bStop = Standard_True;
1773           bOutOfTangentZone = !myIntersectionOn2S.IsTangent();
1774           break;
1775         }
1776       case IntWalk_PasTropGrand:
1777         {
1778           for(parit = 0; parit < 4; parit++) {
1779             Param(parit+1) = SvParam[parit];
1780           }
1781           Standard_Boolean bDecrease = Standard_False;
1782
1783           for(uvit = 0; uvit < 4; uvit++) {
1784             if(pasSav[uvit] < pasInit[uvit]) { 
1785               pasInit[uvit] -= (pasInit[uvit] - pasSav[uvit]) * 0.1;
1786               bDecrease = Standard_True;
1787             }
1788           }
1789
1790           if(bDecrease) nbIterWithoutAppend--;
1791           break;
1792         }
1793       case IntWalk_PointConfondu:
1794         {
1795           for(uvit = 0; uvit < 4; uvit++) {
1796             if(pasuv[uvit] < pasInit[uvit]) {
1797               pasuv[uvit] += (pasInit[uvit] - pasuv[uvit]) * 0.1;
1798             }
1799           }
1800           break;
1801         }
1802       case IntWalk_OK:
1803       case IntWalk_ArretSurPoint:
1804         {
1805           //
1806           bStop = TestArret(theDirectionFlag, Param, ChoixIso);
1807           //
1808
1809           //
1810           if(!bStop) {
1811             Standard_Real u11,v11,u12,v12; 
1812             myIntersectionOn2S.Point().Parameters(u11,v11,u12,v12); 
1813             Standard_Real u21,v21,u22,v22;
1814             previousPoint.Parameters(u21,v21,u22,v22); 
1815
1816             if(((fabs(u11-u21) < ResoU1) && (fabs(v11-v21) < ResoV1)) ||
1817               ((fabs(u12-u22) < ResoU2) && (fabs(v12-v22) < ResoV2))) {
1818                 nbEqualPoints++;
1819             }
1820             else {
1821               nbEqualPoints = 0;
1822             }
1823           }
1824           //
1825
1826           bStop = bStop || !myIntersectionOn2S.IsTangent();
1827           bOutOfTangentZone = !myIntersectionOn2S.IsTangent();
1828
1829           if(!bStop) {
1830             Standard_Boolean pointisvalid = Standard_False;
1831             Standard_Real u1,v1,u2,v2; 
1832             myIntersectionOn2S.Point().Parameters(u1,v1,u2,v2); 
1833
1834             if(u1 <= UM1  && u2 <= UM2 && v1 <= VM1 && 
1835               v2 <= VM2  && u1 >= Um1 && u2 >= Um2 &&
1836               v1 >= Vm1  && v2 >= Vm2) 
1837               pointisvalid = Standard_True;
1838
1839             if(pointisvalid) {
1840               previousPoint = myIntersectionOn2S.Point();
1841               previoustg = myIntersectionOn2S.IsTangent();
1842
1843               if(!previoustg) {
1844                 previousd  = myIntersectionOn2S.Direction();
1845                 previousd1 = myIntersectionOn2S.DirectionOnS1();
1846                 previousd2 = myIntersectionOn2S.DirectionOnS2();
1847               }
1848               Standard_Boolean bAddPoint = Standard_True;
1849
1850               if(line->NbPoints() >= 1) {
1851                 gp_Pnt pf = line->Value(1).Value();
1852                 gp_Pnt pl = previousPoint.Value(); 
1853
1854                 if(pf.Distance(pl) < Precision::Confusion()) { 
1855                   dIncKey++; 
1856                   if(dIncKey == 5000) return bOutOfTangentZone; 
1857                   else bAddPoint = Standard_False;
1858                 }
1859               }
1860
1861               if(bAddPoint) {
1862                 aSeqOfNewPoint.Append(previousPoint);
1863                 nbIterWithoutAppend = 0;
1864               }
1865             }
1866
1867             if (line->NbPoints() == 2) {
1868               for(uvit = 0; uvit < 4; uvit++) {
1869                 pasSav[uvit] = pasuv[uvit]; 
1870               }
1871             }
1872
1873             if ( !pointisvalid ) {
1874               // decrease step if out of bounds
1875               // otherwise the same calculations will be 
1876               // repeated several times
1877               if ( ( u1 > UM1 ) || ( u1 < Um1 ) )
1878                 pasuv[0] *= 0.5;
1879
1880               if ( ( v1 > VM1 ) || ( v1 < Vm1 ) ) 
1881                 pasuv[1] *= 0.5;
1882
1883               if ( ( u2 > UM2 ) || ( u2 < Um2 ) )
1884                 pasuv[2] *= 0.5;
1885
1886               if ( ( v2 > VM2 ) || ( v2 < Vm2 ) )
1887                 pasuv[3] *= 0.5;
1888             }
1889           } // end if(!bStop)
1890           else { //if(bStop)
1891             if(close && (line->NbPoints() >= 1)) { 
1892
1893               if(!bOutOfTangentZone) {
1894                 aSeqOfNewPoint.Append(line->Value(1)); // line end
1895               }
1896               nbIterWithoutAppend = 0;
1897             }
1898             else {
1899               ChoixIso = myIntersectionOn2S.Perform(Param, Rsnld, theChoixIso);
1900
1901               if(myIntersectionOn2S.IsEmpty()) { 
1902                 bStop = Standard_True;// !myIntersectionOn2S.IsTangent();
1903                 bOutOfTangentZone = Standard_False; // !myIntersectionOn2S.IsTangent();
1904               }
1905               else {
1906                 Standard_Boolean bAddPoint = Standard_True;
1907                 Standard_Boolean pointisvalid = Standard_False;
1908
1909                 previousPoint = myIntersectionOn2S.Point();
1910                 Standard_Real u1,v1,u2,v2; 
1911                 previousPoint.Parameters(u1,v1,u2,v2); 
1912
1913                 if(u1 <= UM1  && u2 <= UM2 && v1 <= VM1 && 
1914                   v2 <= VM2  && u1 >= Um1 && u2 >= Um2 &&
1915                   v1 >= Vm1  && v2 >= Vm2) 
1916                   pointisvalid = Standard_True;
1917
1918                 if(pointisvalid) {
1919
1920                   if(line->NbPoints() >= 1) {
1921                     gp_Pnt pf = line->Value(1).Value();
1922                     gp_Pnt pl = previousPoint.Value(); 
1923
1924                     if(pf.Distance(pl) < Precision::Confusion()) { 
1925                       dIncKey++; 
1926                       if(dIncKey == 5000) return bOutOfTangentZone; 
1927                       else bAddPoint = Standard_False;
1928                     }
1929                   }
1930
1931                   if(bAddPoint && !bOutOfTangentZone) {
1932                     aSeqOfNewPoint.Append(previousPoint);
1933                     nbIterWithoutAppend = 0;
1934                   }
1935                 }
1936               }
1937             }
1938           }
1939           break;
1940         }
1941       default:
1942         {
1943           break;
1944         }
1945       }
1946     }
1947   }
1948   Standard_Boolean bExtendLine = Standard_False;
1949   Standard_Real u1 = 0., v1 = 0., u2 = 0., v2 = 0.; 
1950
1951   Standard_Integer pit = 0;
1952
1953   for(pit = 0; !bExtendLine && (pit < 2); pit++) {
1954     if(pit == 0)
1955       previousPoint.Parameters(u1,v1,u2,v2); 
1956     else {
1957       if(aSeqOfNewPoint.Length() > 0)
1958         aSeqOfNewPoint.Value(aSeqOfNewPoint.Length()).Parameters(u1,v1,u2,v2); 
1959       else
1960         break;
1961     }
1962
1963     if(((u1 - Um1) < ResoU1) ||
1964       ((UM1 - u1) < ResoU1) ||
1965       ((u2 - Um2) < ResoU2) ||
1966       ((UM2 - u2) < ResoU2) ||
1967       ((v1 - Vm1) < ResoV1) ||
1968       ((VM1 - v1) < ResoV1) ||
1969       ((v2 - Vm2) < ResoV2) ||
1970       ((VM2 - v2) < ResoV2))
1971       bExtendLine = Standard_True;
1972   }
1973
1974   if(!bExtendLine) {
1975     //    if(aStatus == IntWalk_OK || aStatus == IntWalk_ArretSurPoint) {
1976     if(aStatus == IntWalk_OK) {
1977       bExtendLine = Standard_True;
1978
1979       if(aSeqOfNewPoint.Length() > 1) {
1980         TColStd_Array1OfReal FirstParams(0, 3), LastParams(0, 3), Resolutions(0, 3);
1981         Resolutions(0) = ResoU1; Resolutions(1) = ResoV1; Resolutions(2) = ResoU2; Resolutions(3) = ResoV2;
1982
1983         aSeqOfNewPoint(1).Parameters(FirstParams.ChangeValue(0), FirstParams.ChangeValue(1),
1984           FirstParams.ChangeValue(2), FirstParams.ChangeValue(3));
1985         aSeqOfNewPoint(aSeqOfNewPoint.Length()).Parameters(LastParams.ChangeValue(0), 
1986           LastParams.ChangeValue(1),
1987           LastParams.ChangeValue(2), 
1988           LastParams.ChangeValue(3)); 
1989         Standard_Integer indexofiso = 0;
1990
1991         if(theChoixIso == IntImp_UIsoparametricOnCaro1) indexofiso = 0;
1992         if(theChoixIso == IntImp_VIsoparametricOnCaro1) indexofiso = 1;
1993         if(theChoixIso == IntImp_UIsoparametricOnCaro2) indexofiso = 2;
1994         if(theChoixIso == IntImp_VIsoparametricOnCaro2) indexofiso = 3;
1995
1996         Standard_Integer afirstindex = (indexofiso < 2) ? 0 : 2;
1997         gp_Vec2d aTangentZoneDir(gp_Pnt2d(FirstParams.Value(afirstindex), FirstParams.Value(afirstindex + 1)),
1998           gp_Pnt2d(LastParams.Value(afirstindex), LastParams.Value(afirstindex + 1)));
1999
2000         gp_Dir2d anIsoDir(0, 1);
2001
2002         if((indexofiso == 1) || (indexofiso == 3))
2003           anIsoDir = gp_Dir2d(1, 0);
2004
2005         if(aTangentZoneDir.SquareMagnitude() > gp::Resolution()) {
2006           Standard_Real piquota = M_PI*0.25;
2007
2008           if(fabs(aTangentZoneDir.Angle(anIsoDir)) > piquota) {
2009             Standard_Integer ii = 1, nextii = 2;
2010             gp_Vec2d d1(0, 0);
2011             Standard_Real asqresol = gp::Resolution();
2012             asqresol *= asqresol;
2013
2014             do {
2015               aSeqOfNewPoint(ii).Parameters(FirstParams.ChangeValue(0), FirstParams.ChangeValue(1),
2016                 FirstParams.ChangeValue(2), FirstParams.ChangeValue(3));
2017               aSeqOfNewPoint(ii + 1).Parameters(LastParams.ChangeValue(0), LastParams.ChangeValue(1),
2018                 LastParams.ChangeValue(2), LastParams.ChangeValue(3));
2019               d1 = gp_Vec2d(gp_Pnt2d(FirstParams.Value(afirstindex),
2020                 FirstParams.Value(afirstindex + 1)),
2021                 gp_Pnt2d(LastParams.Value(afirstindex),
2022                 LastParams.Value(afirstindex + 1)));
2023               ii++;
2024             }
2025             while((d1.SquareMagnitude() < asqresol) &&
2026               (ii < aSeqOfNewPoint.Length()));
2027
2028             nextii = ii;
2029
2030             while(nextii < aSeqOfNewPoint.Length()) {
2031
2032               gp_Vec2d nextd1(0, 0);
2033               Standard_Integer jj = nextii;
2034
2035               do {
2036                 aSeqOfNewPoint(jj).Parameters(FirstParams.ChangeValue(0), FirstParams.ChangeValue(1),
2037                   FirstParams.ChangeValue(2), FirstParams.ChangeValue(3));
2038                 aSeqOfNewPoint(jj + 1).Parameters(LastParams.ChangeValue(0), LastParams.ChangeValue(1),
2039                   LastParams.ChangeValue(2), LastParams.ChangeValue(3));
2040                 nextd1 = gp_Vec2d(gp_Pnt2d(FirstParams.Value(afirstindex),
2041                   FirstParams.Value(afirstindex + 1)),
2042                   gp_Pnt2d(LastParams.Value(afirstindex),
2043                   LastParams.Value(afirstindex + 1)));
2044                 jj++;
2045
2046               }
2047               while((nextd1.SquareMagnitude() < asqresol) &&
2048                 (jj < aSeqOfNewPoint.Length()));
2049               nextii = jj;
2050
2051               if(fabs(d1.Angle(nextd1)) > piquota) {
2052                 bExtendLine = Standard_False;
2053                 break;
2054               }
2055               d1 = nextd1;
2056             }
2057           }
2058           // end if(fabs(aTangentZoneDir.Angle(anIsoDir)
2059         }
2060       }
2061     }
2062   }
2063
2064   if(!bExtendLine) {
2065     return Standard_False;
2066   }
2067   Standard_Integer i = 0;
2068
2069   for(i = 1; i <= aSeqOfNewPoint.Length(); i++) {
2070     AddAPoint(aSeqOfNewPoint.Value(i));
2071   }
2072
2073   return bOutOfTangentZone;
2074 }
2075
2076 //=======================================================================
2077 //function : DistanceMinimizeByGradient
2078 //purpose  : 
2079 //
2080 // ATTENTION!!!
2081 //  theInit should be initialized before function calling.
2082 //=======================================================================
2083 Standard_Boolean IntWalk_PWalking::
2084 DistanceMinimizeByGradient( const Handle(Adaptor3d_HSurface)& theASurf1,
2085                            const Handle(Adaptor3d_HSurface)& theASurf2,
2086                            TColStd_Array1OfReal& theInit,
2087                            const Standard_Real* theStep0)
2088 {
2089   const Standard_Integer aNbIterMAX = 60;
2090   const Standard_Real aTol = 1.0e-14;
2091   const Standard_Real aTolNul = 1.0 / Precision::Infinite();
2092
2093   // I.e. if theU1 = 0.0 then Epsilon(theU1) = DBL_MIN (~1.0e-308).
2094   // Work with this number is impossible: there is a dangerous to 
2095   // obtain Floating-point-overflow. Therefore, we limit this value.
2096   const Standard_Real aMinAddValU1 = Max(Epsilon(theInit(1)), aTolNul);
2097   const Standard_Real aMinAddValV1 = Max(Epsilon(theInit(2)), aTolNul);
2098   const Standard_Real aMinAddValU2 = Max(Epsilon(theInit(3)), aTolNul);
2099   const Standard_Real aMinAddValV2 = Max(Epsilon(theInit(4)), aTolNul);
2100
2101   Handle(Geom_Surface) aS1, aS2;
2102
2103   if (theASurf1->GetType() != GeomAbs_BezierSurface &&
2104       theASurf1->GetType() != GeomAbs_BSplineSurface)
2105       return Standard_True;
2106   if (theASurf2->GetType() != GeomAbs_BezierSurface &&
2107       theASurf2->GetType() != GeomAbs_BSplineSurface)
2108       return Standard_True;
2109
2110   Standard_Boolean aStatus = Standard_False;
2111
2112   gp_Pnt aP1, aP2;
2113   gp_Vec aD1u, aD1v, aD2U, aD2V;
2114
2115   theASurf1->D1(theInit(1), theInit(2), aP1, aD1u, aD1v);
2116   theASurf2->D1(theInit(3), theInit(4), aP2, aD2U, aD2V);
2117
2118   Standard_Real aSQDistPrev = aP1.SquareDistance(aP2);
2119
2120   gp_Vec aP12(aP1, aP2);
2121
2122   Standard_Real aGradFu(-aP12.Dot(aD1u));
2123   Standard_Real aGradFv(-aP12.Dot(aD1v));
2124   Standard_Real aGradFU( aP12.Dot(aD2U));
2125   Standard_Real aGradFV( aP12.Dot(aD2V));
2126
2127   Standard_Real aStepU1 = 1.0e-6, aStepV1 = 1.0e-6,
2128                 aStepU2 = 1.0e-6, aStepV2 = 1.0e-6;
2129     
2130   if (theStep0)
2131   {
2132     aStepU1 = theStep0[0];
2133     aStepV1 = theStep0[1];
2134     aStepU2 = theStep0[2];
2135     aStepV2 = theStep0[3];
2136   }
2137
2138   Standard_Boolean flRepeat = Standard_True;
2139   Standard_Integer aNbIter = aNbIterMAX;
2140
2141   while(flRepeat)
2142   {
2143     Standard_Real anAdd = aGradFu*aStepU1;
2144     const Standard_Real aPARu = theInit(1) - Sign(Max(Abs(anAdd), aMinAddValU1), anAdd);
2145     
2146     anAdd = aGradFv*aStepV1;
2147     const Standard_Real aPARv = theInit(2) - Sign(Max(Abs(anAdd), aMinAddValV1), anAdd);
2148
2149     anAdd = aGradFU*aStepU2;
2150     const Standard_Real aParU = theInit(3) - Sign(Max(Abs(anAdd), aMinAddValU2), anAdd);
2151
2152     anAdd = aGradFV*aStepV2;
2153     const Standard_Real aParV = theInit(4) - Sign(Max(Abs(anAdd), aMinAddValV2), anAdd);
2154
2155     gp_Pnt aPt1, aPt2;
2156
2157     theASurf1->D1(aPARu, aPARv, aPt1, aD1u, aD1v);
2158     theASurf2->D1(aParU, aParV, aPt2, aD2U, aD2V);
2159
2160     Standard_Real aSQDist = aPt1.SquareDistance(aPt2);
2161
2162     if(aSQDist < aSQDistPrev)
2163     {
2164       aSQDistPrev = aSQDist;
2165       theInit(1) = aPARu;
2166       theInit(2) = aPARv;
2167       theInit(3) = aParU;
2168       theInit(4) = aParV;
2169
2170       aStatus = aSQDistPrev < aTol;
2171       aStepU1 *= 1.2;
2172       aStepV1 *= 1.2;
2173       aStepU2 *= 1.2;
2174       aStepV2 *= 1.2;
2175     }
2176     else
2177     {
2178       if(--aNbIter < 0)
2179       {
2180         flRepeat = Standard_False;
2181       }
2182       else
2183       {
2184         theASurf1->D1(theInit(1), theInit(2), aPt1, aD1u, aD1v);
2185         theASurf2->D1(theInit(3), theInit(4), aPt2, aD2U, aD2V);
2186
2187         gp_Vec aPt12(aPt1, aPt2);
2188         aGradFu = -aPt12.Dot(aD1u);
2189         aGradFv = -aPt12.Dot(aD1v);
2190         aGradFU = aPt12.Dot(aD2U);
2191         aGradFV = aPt12.Dot(aD2V);
2192
2193         if (theStep0)
2194         {
2195           aStepU1 = theStep0[0];
2196           aStepV1 = theStep0[1];
2197           aStepU2 = theStep0[2];
2198           aStepV2 = theStep0[3];
2199         }
2200         else
2201         {
2202           aStepU1 = aStepV1 = aStepU2 = aStepV2 = 1.0e-6;
2203         }
2204       }
2205     }
2206   }
2207
2208   return aStatus;
2209 }
2210
2211 //=======================================================================
2212 //function : DistanceMinimizeByExtrema
2213 //purpose  : 
2214 //
2215 // ATTENTION!!!
2216 //  theP0, theU0 and theV0 parameters should be initialized
2217 //    before the function calling.
2218 //=======================================================================
2219 Standard_Boolean IntWalk_PWalking::
2220 DistanceMinimizeByExtrema(const Handle(Adaptor3d_HSurface)& theASurf, 
2221                           const gp_Pnt& theP0,
2222                           Standard_Real& theU0,
2223                           Standard_Real& theV0,
2224                           const Standard_Real* theStep0)
2225 {
2226   const Standard_Real aTol = 1.0e-14;
2227   gp_Pnt aPS;
2228   gp_Vec aD1Su, aD1Sv, aD2Su, aD2Sv, aD2SuvTemp;
2229   Standard_Real aSQDistPrev = RealLast();
2230   Standard_Real aU = theU0, aV = theV0;
2231
2232   const Standard_Real aStep0[2] = { theStep0 ? theStep0[0] : 1.0,
2233                                     theStep0 ? theStep0[1] : 1.0 };
2234
2235   Standard_Integer aNbIter = 10;
2236   do
2237   {
2238     theASurf->D2(aU, aV, aPS, aD1Su, aD1Sv, aD2Su, aD2Sv, aD2SuvTemp);
2239
2240     gp_Vec aVec(theP0, aPS);
2241
2242     Standard_Real aSQDist = aVec.SquareMagnitude();
2243
2244     if(aSQDist >= aSQDistPrev)
2245       break;
2246
2247     aSQDistPrev = aSQDist;
2248     theU0 = aU;
2249     theV0 = aV;
2250     aNbIter--;
2251
2252     if(aSQDistPrev < aTol)
2253       break;
2254
2255     //Functions
2256     const Standard_Real aF1 = aD1Su.Dot(aVec), aF2 = aD1Sv.Dot(aVec);
2257
2258     //Derivatives
2259     const Standard_Real aDf1u = aD2Su.Dot(aVec) + aD1Su.Dot(aD1Su),
2260       aDf1v = aD2Su.Dot(aD1Sv),
2261       aDf2u = aDf1v,
2262       aDf2v = aD2Sv.Dot(aVec) + aD1Sv.Dot(aD1Sv);
2263
2264     const Standard_Real aDet = aDf1u*aDf2v - aDf1v*aDf2u;
2265     aU -= aStep0[0]*(aDf2v*aF1 - aDf1v*aF2) / aDet;
2266     aV += aStep0[1]*(aDf2u*aF1 - aDf1u*aF2) / aDet;
2267   }
2268   while(aNbIter > 0);
2269
2270   return (aSQDistPrev < aTol);
2271 }
2272
2273 //=======================================================================
2274 //function : HandleSingleSingularPoint
2275 //purpose  : 
2276 //=======================================================================
2277 Standard_Boolean IntWalk_PWalking::HandleSingleSingularPoint(const Handle(Adaptor3d_HSurface)& theASurf1,
2278                                                              const Handle(Adaptor3d_HSurface)& theASurf2,
2279                                                              const Standard_Real the3DTol,
2280                                                              TColStd_Array1OfReal &thePnt)
2281 {
2282   // u1, v1, u2, v2 order is used.
2283   Standard_Real aLowBorder[4] = {theASurf1->FirstUParameter(),
2284                                  theASurf1->FirstVParameter(),
2285                                  theASurf2->FirstUParameter(),
2286                                  theASurf2->FirstVParameter()};
2287   Standard_Real aUppBorder[4] = {theASurf1->LastUParameter(),
2288                                  theASurf1->LastVParameter(),
2289                                  theASurf2->LastUParameter(),
2290                                  theASurf2->LastVParameter()};
2291   IntImp_ConstIsoparametric aLockedDir[4] = {IntImp_UIsoparametricOnCaro1,
2292                                              IntImp_VIsoparametricOnCaro1,
2293                                              IntImp_UIsoparametricOnCaro2,
2294                                              IntImp_VIsoparametricOnCaro2};
2295
2296   // Create new intersector with new tolerance.
2297   IntWalk_TheInt2S anInt(theASurf1, theASurf2, the3DTol);
2298   math_FunctionSetRoot aRsnld(anInt.Function());
2299
2300   for (Standard_Integer i = 1; i <= 4; ++i)
2301   {
2302     if ( Abs(thePnt(i) - aLowBorder[i - 1]) < Precision::PConfusion() ||
2303          Abs(thePnt(i) - aUppBorder[i - 1]) < Precision::PConfusion())
2304     {
2305
2306       anInt.Perform(thePnt,aRsnld, aLockedDir[i - 1]);
2307
2308       if (!anInt.IsDone())
2309         continue;
2310
2311       if (anInt.IsEmpty())
2312         continue;
2313
2314       anInt.Point().Parameters(thePnt(1), thePnt(2), thePnt(3), thePnt(4));
2315
2316       Standard_Boolean isInDomain = Standard_True;
2317       for (Standard_Integer j = 1; isInDomain && (j <= 4); ++j)
2318       {
2319         if ((thePnt(j) - aLowBorder[j - 1] + Precision::PConfusion())*
2320             (thePnt(j) - aUppBorder[j - 1] - Precision::PConfusion()) > 0.0)
2321         {
2322           isInDomain = Standard_False;
2323         }
2324       }
2325
2326       if (isInDomain)
2327         return Standard_True;
2328     }
2329   }
2330
2331   return Standard_False;
2332 }
2333
2334 //=======================================================================
2335 //function : SeekPointOnBoundary
2336 //purpose  : 
2337 //=======================================================================
2338 Standard_Boolean IntWalk_PWalking::
2339         SeekPointOnBoundary(const Handle(Adaptor3d_HSurface)& theASurf1,
2340                             const Handle(Adaptor3d_HSurface)& theASurf2,
2341                             const Standard_Real theU1,
2342                             const Standard_Real theV1,
2343                             const Standard_Real theU2,
2344                             const Standard_Real theV2,
2345                             const Standard_Boolean isTheFirst)
2346 {
2347   Standard_Boolean isOK = Standard_False;
2348
2349   // u1, v1, u2, v2 order is used.
2350   const Standard_Real aLowBorder[4] = {theASurf1->FirstUParameter(),
2351                                        theASurf1->FirstVParameter(),
2352                                        theASurf2->FirstUParameter(),
2353                                        theASurf2->FirstVParameter()};
2354   const Standard_Real aUppBorder[4] = {theASurf1->LastUParameter(),
2355                                        theASurf1->LastVParameter(),
2356                                        theASurf2->LastUParameter(),
2357                                        theASurf2->LastVParameter()};
2358
2359   // Tune solution tolerance according with object size.
2360   const Standard_Real aRes1 = Max(Precision::PConfusion() / theASurf1->UResolution(1.0),
2361                                   Precision::PConfusion() / theASurf1->VResolution(1.0));
2362   const Standard_Real aRes2 = Max(Precision::PConfusion() / theASurf2->UResolution(1.0),
2363                                   Precision::PConfusion() / theASurf2->VResolution(1.0));
2364   const Standard_Real a3DTol = Max(aRes1, aRes2);
2365   const Standard_Real aTol = Max(Precision::Confusion(), a3DTol);
2366
2367   // u1, v1, u2, v2 order is used.
2368   TColStd_Array1OfReal aPnt(1,4);
2369   aPnt(1) = theU1; aPnt(2) = theV1; aPnt(3) = theU2; aPnt(4) = theV2;
2370   TColStd_Array1OfReal aSingularPnt(aPnt);
2371
2372   Standard_Integer aNbIter = 20;
2373   Standard_Boolean aStatus = Standard_False;
2374   do
2375   {
2376     aNbIter--;
2377     aStatus = DistanceMinimizeByGradient(theASurf1, theASurf2, aPnt);
2378     if (aStatus && !AdjustToDomain(4, &aPnt(1), &aLowBorder[0], &aUppBorder[0]))
2379       break;
2380
2381     aStatus = DistanceMinimizeByExtrema(theASurf1, theASurf2->Value(aPnt(3), aPnt(4)),
2382                                         aPnt(1), aPnt(2));
2383     if (aStatus && !AdjustToDomain(2, &aPnt(1), &aLowBorder[0], &aUppBorder[0]))
2384       break;
2385
2386     aStatus = DistanceMinimizeByExtrema(theASurf2, theASurf1->Value(aPnt(1), aPnt(2)),
2387                                         aPnt(3), aPnt(4));
2388     if (aStatus && !AdjustToDomain(2, &aPnt(3), &aLowBorder[2], &aUppBorder[2]))
2389       break;
2390   }
2391   while(!aStatus && (aNbIter > 0));
2392
2393   // Handle singular points.
2394   Standard_Boolean aSingularStatus = HandleSingleSingularPoint(theASurf1,
2395                                                                theASurf2,
2396                                                                aTol,
2397                                                                aSingularPnt);
2398   if (aSingularStatus)
2399     aPnt = aSingularPnt;
2400
2401   if (!aStatus && !aSingularStatus)
2402   {
2403     return isOK;
2404   }
2405
2406   gp_Pnt aP1 = theASurf1->Value(aPnt(1), aPnt(2));
2407   gp_Pnt aP2 = theASurf2->Value(aPnt(3), aPnt(4));
2408   const gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));
2409
2410   const Standard_Real aSQDist = aPInt.SquareDistance(aP1);
2411   if (aSQDist > aTol * aTol)
2412   {
2413     return isOK;
2414   }
2415
2416   //Found point is true intersection point
2417   IntSurf_PntOn2S anIP;
2418   anIP.SetValue(aPInt, aPnt(1), aPnt(2), aPnt(3), aPnt(4));
2419
2420   //The main idea of checks below is to define if insertion of
2421   //addition point (on the boundary) does not lead to invalid
2422   //intersection curve (e.g. having a loop).
2423   //
2424   //Loops are detected with rotation angle of the Walking-line (WL).
2425   //If there is hairpin bend then insertion is forbidden.
2426
2427   //There are at least two possible problems:
2428   //  1. There are some cases when two neighbor points of the WL
2429   //      are almost coincident (the distance between them is less
2430   //      than Precision::Confusion). It is impossible to define
2431   //      rotation angle in these cases. Therefore, points with
2432   //      "good" distances should be selected.
2433
2434   //  2. Intersection point on the surface boundary has highest
2435   //      priority in compare with other "middle" points. Therefore,
2436   //      if insertion of new point will result in a bend then some
2437   //      "middle" points should be deleted in order to provide
2438   //      correct insertion.
2439
2440   //Problem test cases:
2441   //  test bugs modalg_5 bug24585_1
2442   //  test boolean bcut_complex G7
2443   //  test bugs moddata_2 bug469
2444
2445   if (isTheFirst)
2446   {
2447     while (line->NbPoints() > 1)
2448     {
2449       const Standard_Integer aNbPnts = line->NbPoints();
2450
2451       Standard_Integer aPInd = 1;
2452       for (; aPInd <= aNbPnts; aPInd++)
2453       {
2454         aP1.SetXYZ(line->Value(aPInd).Value().XYZ());
2455         if (aP1.SquareDistance(aPInt) > Precision::SquareConfusion())
2456         {
2457           break;
2458         }
2459         else if (aPInd == 1)
2460         {
2461           // After insertion, we will obtain
2462           // two coincident points in the line.
2463           // Therefore, insertion is forbidden.
2464           return isOK;
2465         }
2466       }
2467
2468       for (++aPInd; aPInd <= aNbPnts; aPInd++)
2469       {
2470         aP2.SetXYZ(line->Value(aPInd).Value().XYZ());
2471         if (aP1.SquareDistance(aP2) > Precision::SquareConfusion())
2472           break;
2473       }
2474
2475       if (aPInd > aNbPnts)
2476       {
2477         return isOK;
2478       }
2479
2480       const gp_XYZ aDir01(aP1.XYZ() - aPInt.XYZ());
2481       const gp_XYZ aDir12(aP2.XYZ() - aP1.XYZ());
2482
2483       if (aDir01.Dot(aDir12) > 0.0)
2484       {
2485         break;
2486       }
2487
2488       RemoveAPoint(1);
2489     }
2490
2491     line->InsertBefore(1, anIP);
2492     isOK = Standard_True;
2493   }
2494   else
2495   {
2496     while (line->NbPoints() > 1)
2497     {
2498       const Standard_Integer aNbPnts = line->NbPoints();
2499
2500       gp_Pnt aPPrev, aPCurr;
2501       Standard_Integer aPInd = aNbPnts;
2502       for (; aPInd > 0; aPInd--)
2503       {
2504         aPCurr.SetXYZ(line->Value(aPInd).Value().XYZ());
2505         if (aPCurr.SquareDistance(aPInt) > Precision::SquareConfusion())
2506         {
2507           break;
2508         }
2509         else if (aPInd == aNbPnts)
2510         {
2511           // After insertion, we will obtain
2512           // two coincident points in the line.
2513           // Therefore, insertion is forbidden.
2514           return isOK;
2515         }
2516       }
2517
2518       for (--aPInd; aPInd > 0; aPInd--)
2519       {
2520         aPPrev.SetXYZ(line->Value(aPInd).Value().XYZ());
2521         if (aPCurr.SquareDistance(aPPrev) > Precision::SquareConfusion())
2522           break;
2523       }
2524
2525       if (aPInd < 1)
2526       {
2527         return isOK;
2528       }
2529
2530       const gp_XYZ aDirPC(aPCurr.XYZ() - aPPrev.XYZ());
2531       const gp_XYZ aDirCN(aPInt.XYZ() - aPCurr.XYZ());
2532
2533       if (aDirPC.Dot(aDirCN) > 0.0)
2534       {
2535         break;
2536       }
2537
2538       RemoveAPoint(aNbPnts);
2539     }
2540
2541     line->Add(anIP);
2542     isOK = Standard_True;
2543   }
2544
2545   return isOK;
2546 }
2547
2548 //=======================================================================
2549 //function : PutToBoundary
2550 //purpose  : 
2551 //=======================================================================
2552 Standard_Boolean IntWalk_PWalking::
2553 PutToBoundary(const Handle(Adaptor3d_HSurface)& theASurf1,
2554               const Handle(Adaptor3d_HSurface)& theASurf2)
2555 {
2556   const Standard_Real aTolMin = Precision::Confusion();
2557
2558   Standard_Boolean hasBeenAdded = Standard_False;
2559
2560   const Standard_Real aU1bFirst = theASurf1->FirstUParameter();
2561   const Standard_Real aU1bLast = theASurf1->LastUParameter();
2562   const Standard_Real aU2bFirst = theASurf2->FirstUParameter();
2563   const Standard_Real aU2bLast = theASurf2->LastUParameter();
2564   const Standard_Real aV1bFirst = theASurf1->FirstVParameter();
2565   const Standard_Real aV1bLast = theASurf1->LastVParameter();
2566   const Standard_Real aV2bFirst = theASurf2->FirstVParameter();
2567   const Standard_Real aV2bLast = theASurf2->LastVParameter();
2568
2569   Standard_Real aTol = 1.0;
2570   aTol = Min(aTol, aU1bLast - aU1bFirst);
2571   aTol = Min(aTol, aU2bLast - aU2bFirst);
2572   aTol = Min(aTol, aV1bLast - aV1bFirst);
2573   aTol = Min(aTol, aV2bLast - aV2bFirst)*1.0e-3;
2574
2575   if(aTol <= 2.0*aTolMin)
2576     return hasBeenAdded;
2577
2578   Standard_Boolean isNeedAdding = Standard_False;
2579   Standard_Boolean isU1parallel = Standard_False, isV1parallel = Standard_False;
2580   Standard_Boolean isU2parallel = Standard_False, isV2parallel = Standard_False;
2581   IsParallel(line, Standard_True, aTol, isU1parallel, isV1parallel);
2582   IsParallel(line, Standard_False, aTol, isU2parallel, isV2parallel);
2583
2584   Standard_Real u1, v1, u2, v2;
2585   line->Value(1).Parameters(u1, v1, u2, v2);
2586   Standard_Real aDelta = 0.0;
2587
2588   if(!isV1parallel)
2589   {
2590     aDelta = u1 - aU1bFirst;
2591     if((aTolMin < aDelta) && (aDelta < aTol))
2592     {
2593       u1 = aU1bFirst;
2594       isNeedAdding = Standard_True;
2595     }
2596     else
2597     {
2598       aDelta = aU1bLast - u1;
2599       if((aTolMin < aDelta) && (aDelta < aTol))
2600       {
2601         u1 = aU1bLast;
2602         isNeedAdding = Standard_True;
2603       }
2604     }
2605   }
2606
2607   if(!isV2parallel)
2608   {
2609     aDelta = u2 - aU2bFirst;
2610     if((aTolMin < aDelta) && (aDelta < aTol))
2611     {
2612       u2 = aU2bFirst;
2613       isNeedAdding = Standard_True;
2614     }
2615     else
2616     {
2617       aDelta = aU2bLast - u2;
2618       if((aTolMin < aDelta) && (aDelta < aTol))
2619       {
2620         u2 = aU2bLast;
2621         isNeedAdding = Standard_True;
2622       }
2623     }
2624   }
2625
2626   if(!isU1parallel)
2627   {
2628     aDelta = v1 - aV1bFirst;
2629     if((aTolMin < aDelta) && (aDelta < aTol))
2630     {
2631       v1 = aV1bFirst;
2632       isNeedAdding = Standard_True;
2633     }
2634     else
2635     {
2636       aDelta = aV1bLast - v1;
2637       if((aTolMin < aDelta) && (aDelta < aTol))
2638       {
2639         v1 = aV1bLast;
2640         isNeedAdding = Standard_True;
2641       }
2642     }
2643   }
2644
2645   if(!isU2parallel)
2646   {
2647     aDelta = v2 - aV2bFirst;
2648     if((aTolMin < aDelta) && (aDelta < aTol))
2649     {
2650       v2 = aV2bFirst;
2651       isNeedAdding = Standard_True;
2652     }
2653     else
2654     {
2655       aDelta = aV2bLast - v2;
2656       if((aTolMin < aDelta) && (aDelta < aTol))
2657       {
2658         v2 = aV2bLast;
2659         isNeedAdding = Standard_True;
2660       }
2661     }
2662   }
2663
2664   if(isNeedAdding)
2665   {
2666     hasBeenAdded = 
2667       SeekPointOnBoundary(theASurf1, theASurf2, u1, 
2668       v1, u2, v2, Standard_True);
2669   }
2670
2671   const Standard_Integer aNbPnts = line->NbPoints();
2672   isNeedAdding = Standard_False;
2673   line->Value(aNbPnts).Parameters(u1, v1, u2, v2);
2674
2675   if(!isV1parallel)
2676   {
2677     aDelta = u1 - aU1bFirst;
2678     if((aTolMin < aDelta) && (aDelta < aTol))
2679     {
2680       u1 = aU1bFirst;
2681       isNeedAdding = Standard_True;
2682     }
2683     else
2684     {
2685       aDelta = aU1bLast - u1;
2686       if((aTolMin < aDelta) && (aDelta < aTol))
2687       {
2688         u1 = aU1bLast;
2689         isNeedAdding = Standard_True;
2690       }
2691     }
2692   }
2693
2694   if(!isV2parallel)
2695   {
2696     aDelta = u2 - aU2bFirst;
2697     if((aTolMin < aDelta) && (aDelta < aTol))
2698     {
2699       u2 = aU2bFirst;
2700       isNeedAdding = Standard_True;
2701     }
2702     else
2703     {
2704       aDelta = aU2bLast - u2;
2705       if((aTolMin < aDelta) && (aDelta < aTol))
2706       {
2707         u2 = aU2bLast;
2708         isNeedAdding = Standard_True;
2709       }
2710     }
2711   }
2712
2713   if(!isU1parallel)
2714   {
2715     aDelta = v1 - aV1bFirst;
2716     if((aTolMin < aDelta) && (aDelta < aTol))
2717     {
2718       v1 = aV1bFirst;
2719       isNeedAdding = Standard_True;
2720     }
2721     else
2722     {
2723       aDelta = aV1bLast - v1;
2724       if((aTolMin < aDelta) && (aDelta < aTol))
2725       {
2726         v1 = aV1bLast;
2727         isNeedAdding = Standard_True;
2728       }
2729     }
2730   }
2731
2732   if(!isU2parallel)
2733   {
2734     aDelta = v2 - aV2bFirst;
2735     if((aTolMin < aDelta) && (aDelta < aTol))
2736     {
2737       v2 = aV2bFirst;
2738       isNeedAdding = Standard_True;
2739     }
2740     else
2741     {
2742       aDelta = aV2bLast - v2;
2743       if((aTolMin < aDelta) && (aDelta < aTol))
2744       {
2745         v2 = aV2bLast;
2746         isNeedAdding = Standard_True;
2747       }
2748     }
2749   }
2750
2751   if(isNeedAdding)
2752   {
2753     hasBeenAdded = 
2754       SeekPointOnBoundary(theASurf1, theASurf2, u1, 
2755       v1, u2, v2, Standard_False);
2756   }
2757
2758   return hasBeenAdded;
2759 }
2760
2761 //=======================================================================
2762 //function : SeekAdditionalPoints
2763 //purpose  : 
2764 //=======================================================================
2765 Standard_Boolean IntWalk_PWalking::
2766 SeekAdditionalPoints( const Handle(Adaptor3d_HSurface)& theASurf1,
2767                      const Handle(Adaptor3d_HSurface)& theASurf2,
2768                      const Standard_Integer theMinNbPoints)
2769 {
2770   const Standard_Real aTol = 1.0e-14;
2771   Standard_Integer aNbPoints = line->NbPoints();
2772   if(aNbPoints > theMinNbPoints)
2773     return Standard_True;
2774
2775   const Standard_Real aU1bFirst = theASurf1->FirstUParameter();
2776   const Standard_Real aU1bLast = theASurf1->LastUParameter();
2777   const Standard_Real aU2bFirst = theASurf2->FirstUParameter();
2778   const Standard_Real aU2bLast = theASurf2->LastUParameter();
2779   const Standard_Real aV1bFirst = theASurf1->FirstVParameter();
2780   const Standard_Real aV1bLast = theASurf1->LastVParameter();
2781   const Standard_Real aV2bFirst = theASurf2->FirstVParameter();
2782   const Standard_Real aV2bLast = theASurf2->LastVParameter();
2783
2784
2785   Standard_Boolean isPrecise = Standard_False;
2786
2787   TColStd_Array1OfReal aPnt(1, 4);
2788   aPnt.Init(0.0);
2789
2790   Standard_Integer aNbPointsPrev = 0;
2791   while(aNbPoints < theMinNbPoints && (aNbPoints != aNbPointsPrev))
2792   {
2793     aNbPointsPrev = aNbPoints;
2794     for(Standard_Integer fp = 1, lp = 2; fp < aNbPoints; fp = lp + 1)
2795     {
2796       Standard_Real U1f, V1f, U2f, V2f; //first point in 1st and 2nd surafaces
2797       Standard_Real U1l, V1l, U2l, V2l; //last  point in 1st and 2nd surafaces
2798
2799       lp = fp+1;
2800       line->Value(fp).Parameters(U1f, V1f, U2f, V2f);
2801       line->Value(lp).Parameters(U1l, V1l, U2l, V2l);
2802
2803       aPnt(1) = 0.5*(U1f + U1l);
2804       if(aPnt(1) < aU1bFirst)
2805         aPnt(1) = aU1bFirst;
2806       if(aPnt(1) > aU1bLast)
2807         aPnt(1) = aU1bLast;
2808
2809       aPnt(2) = 0.5*(V1f+V1l);
2810       if(aPnt(2) < aV1bFirst)
2811         aPnt(2) = aV1bFirst;
2812       if(aPnt(2) > aV1bLast)
2813         aPnt(2) = aV1bLast;
2814
2815       aPnt(3) = 0.5*(U2f+U2l);
2816       if(aPnt(3) < aU2bFirst)
2817         aPnt(3) = aU2bFirst;
2818       if(aPnt(3) > aU2bLast)
2819         aPnt(3) = aU2bLast;
2820
2821       aPnt(4) = 0.5*(V2f+V2l);
2822       if(aPnt(4) < aV2bFirst)
2823         aPnt(4) = aV2bFirst;
2824       if(aPnt(4) > aV2bLast)
2825         aPnt(4) = aV2bLast;
2826
2827       Standard_Boolean aStatus = Standard_False;
2828       Standard_Integer aNbIter = 5;
2829       do
2830       {
2831         aStatus = DistanceMinimizeByGradient(theASurf1, theASurf2, aPnt);
2832         if(aStatus)
2833         {
2834           break;
2835         }
2836
2837         aStatus = DistanceMinimizeByExtrema(theASurf1, theASurf2->Value(aPnt(3), aPnt(4)), aPnt(1), aPnt(2));
2838         if(aStatus)
2839         {
2840           break;
2841         }
2842
2843         aStatus = DistanceMinimizeByExtrema(theASurf2, theASurf1->Value(aPnt(1), aPnt(2)), aPnt(3), aPnt(4));
2844         if(aStatus)
2845         {
2846           break;
2847         }
2848       }
2849       while(!aStatus && (--aNbIter > 0));
2850
2851       if(aStatus)
2852       {
2853         gp_Pnt  aP1 = theASurf1->Value(aPnt(1), aPnt(2)),
2854           aP2 = theASurf2->Value(aPnt(3), aPnt(4));
2855         gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));
2856
2857         const Standard_Real aSQDist1 = aPInt.SquareDistance(aP1),
2858           aSQDist2 = aPInt.SquareDistance(aP2);
2859
2860         if((aSQDist1 < aTol) && (aSQDist2 < aTol))
2861         {
2862           IntSurf_PntOn2S anIP;
2863           anIP.SetValue(aPInt, aPnt(1), aPnt(2), aPnt(3), aPnt(4));
2864           line->InsertBefore(lp, anIP);
2865
2866           isPrecise = Standard_True;
2867
2868           if(++aNbPoints >= theMinNbPoints)
2869             break;
2870         }
2871         else
2872         {
2873           lp--;
2874         }
2875       }
2876     }
2877   }
2878
2879   return isPrecise;
2880 }
2881
2882 void IntWalk_PWalking::
2883 RepartirOuDiviser(Standard_Boolean& DejaReparti,
2884                   IntImp_ConstIsoparametric& ChoixIso,
2885                   Standard_Boolean& Arrive) 
2886
2887                   // at the neighborhood of a point, there is a fail of marching 
2888                   // it is required to divide the steps to try to continue
2889                   // if the step is too small if we are on border
2890                   // restart in another direction if it was not done, otherwise stop
2891
2892 {
2893   //  Standard_Integer i;
2894   if (Arrive) {    //restart in the other direction
2895     if (!DejaReparti ) {
2896       Arrive = Standard_False;
2897       DejaReparti = Standard_True;
2898       previousPoint = line->Value(1);
2899       previoustg = Standard_False;
2900       previousd1 = firstd1;
2901       previousd2 = firstd2;
2902       previousd = tgdir;
2903       myTangentIdx = line->NbPoints();
2904       tgdir.Reverse();
2905       line->Reverse();
2906
2907       //-- printf("\nIntWalk_PWalking_2.gxx Reverse %3d\n",indextg);
2908       sensCheminement = -1;
2909       tgfirst      = tglast;
2910       tglast       = Standard_False;
2911       ChoixIso     = choixIsoSav;
2912 #if 0
2913       pasuv[0]=pasSav[0];
2914       pasuv[1]=pasSav[1];
2915       pasuv[2]=pasSav[2];
2916       pasuv[3]=pasSav[3];
2917 #else 
2918       Standard_Real u1,v1,u2,v2;
2919       Standard_Real U1,V1,U2,V2;
2920       Standard_Integer nn=line->NbPoints();
2921       if(nn>2) { 
2922         line->Value(nn).Parameters(u1,v1,u2,v2);
2923         line->Value(nn-1).Parameters(U1,V1,U2,V2);
2924         pasuv[0]=Abs(u1-U1);
2925         pasuv[1]=Abs(v1-V1);
2926         pasuv[2]=Abs(u2-U2);
2927         pasuv[3]=Abs(v2-V2);
2928       }
2929 #endif
2930
2931     }
2932   }  
2933   else  {
2934     if (    pasuv[0]*0.5 < ResoU1
2935       &&  pasuv[1]*0.5 < ResoV1
2936       &&  pasuv[2]*0.5 < ResoU2
2937       &&  pasuv[3]*0.5 < ResoV2
2938       ) {
2939         if (!previoustg) {
2940           tglast = Standard_True;      // IS IT ENOUGH ????
2941         }
2942
2943         if (!DejaReparti)
2944         {  //restart in the other direction
2945           DejaReparti = Standard_True;
2946           previousPoint = line->Value(1);
2947           previoustg = Standard_False;
2948           previousd1 = firstd1;
2949           previousd2 = firstd2;
2950           previousd = tgdir;
2951           myTangentIdx = line->NbPoints();
2952           tgdir.Reverse();
2953           line->Reverse();
2954
2955           //-- printf("\nIntWalk_PWalking_2.gxx Reverse %3d\n",indextg);
2956
2957           sensCheminement   = -1;
2958           tgfirst           = tglast;
2959           tglast            = Standard_False;
2960           ChoixIso          = choixIsoSav;
2961
2962 #if 0 
2963           pasuv[0]=pasSav[0];
2964           pasuv[1]=pasSav[1];
2965           pasuv[2]=pasSav[2];
2966           pasuv[3]=pasSav[3];
2967 #else 
2968           Standard_Real u1,v1,u2,v2;
2969           Standard_Real U1,V1,U2,V2;
2970           Standard_Integer nn=line->NbPoints();
2971           if(nn>2) { 
2972             line->Value(nn).Parameters(u1,v1,u2,v2);
2973             line->Value(nn-1).Parameters(U1,V1,U2,V2);
2974             pasuv[0]=Abs(u1-U1);
2975             pasuv[1]=Abs(v1-V1);
2976             pasuv[2]=Abs(u2-U2);
2977             pasuv[3]=Abs(v2-V2);
2978           }
2979 #endif
2980         }
2981         else Arrive = Standard_True;
2982     }
2983     else {
2984       pasuv[0]*=0.5;
2985       pasuv[1]*=0.5;
2986       pasuv[2]*=0.5;
2987       pasuv[3]*=0.5; 
2988     }
2989   }
2990 }
2991
2992 namespace {
2993   //OCC431(apo): modified ->
2994   static const Standard_Real CosRef2D =  Cos(M_PI/9.0),  AngRef2D = M_PI/2.0; 
2995
2996   static const Standard_Real d = 7.0;
2997 }
2998
2999 IntWalk_StatusDeflection  IntWalk_PWalking::TestDeflection(const IntImp_ConstIsoparametric choixIso,
3000                                                            const IntWalk_StatusDeflection  theStatus)
3001
3002 // test if vector is observed by calculating an increase of vector 
3003 //     or the previous point and its tangent, the new calculated point and its  
3004 //     tangent; it is possible to find a cube passing by the 2 points and having as a 
3005 //     derivative the tangents of the intersection
3006 //     calculate the point with parameter 0.5 on cube=p1 
3007 //     calculate the medium point of 2 points of intersection=p2
3008 //   if arrow/2<=||p1p2||<= arrow consider that the vector is observed
3009 //   otherwise adjust the step depending on the ratio ||p1p2||/vector
3010 //   and the previous step 
3011 // test if in  2 tangent planes of surfaces there is no too great angle2d 
3012 // grand : if yes divide the step
3013 // test if there is no change of side
3014 //  
3015 {
3016   if(line->NbPoints() ==1 ) { 
3017     STATIC_BLOCAGE_SUR_PAS_TROP_GRAND=STATIC_PRECEDENT_INFLEXION=0;
3018   }
3019
3020   IntWalk_StatusDeflection aStatus = IntWalk_OK;
3021   Standard_Real FlecheCourante , Ratio = 1.0;
3022
3023   // Caro1 and Caro2
3024   const Handle(Adaptor3d_HSurface)& Caro1 = myIntersectionOn2S.Function().AuxillarSurface1();
3025   const Handle(Adaptor3d_HSurface)& Caro2 = myIntersectionOn2S.Function().AuxillarSurface2();
3026
3027   const IntSurf_PntOn2S& CurrentPoint = myIntersectionOn2S.Point(); 
3028   //==================================================================================
3029   //=========               S t o p   o n   p o i n t                 ============
3030   //================================================================================== 
3031   if (myIntersectionOn2S.IsTangent())  { 
3032     return IntWalk_ArretSurPoint;  
3033   }
3034
3035   const gp_Dir& TgCourante = myIntersectionOn2S.Direction();
3036
3037   const Standard_Real aCosBetweenTangent = TgCourante.Dot(previousd);
3038
3039   //==================================================================================
3040   //=========   R i s k   o f    i n f l e x i o n   p o i n t  ============
3041   //==================================================================================  
3042   if (aCosBetweenTangent < 0) {
3043     //------------------------------------------------------------
3044     //-- Risk of inflexion point : Divide the step by 2
3045     //-- Initialize STATIC_PRECEDENT_INFLEXION so that 
3046     //-- at the next call to return Pas_OK if there is no 
3047     //-- more risk of the point of inflexion
3048     //------------------------------------------------------------
3049
3050     pasuv[0]*=0.5;
3051     pasuv[1]*=0.5;
3052     pasuv[2]*=0.5;
3053     pasuv[3]*=0.5;
3054     STATIC_PRECEDENT_INFLEXION+=3; 
3055     if (pasuv[0] < ResoU1 && pasuv[1] <ResoV1 && pasuv[2] <ResoU2 && pasuv[3] < ResoV2)
3056       return IntWalk_ArretSurPointPrecedent;
3057     else 
3058       return IntWalk_PasTropGrand;
3059   }
3060   else {
3061     if(STATIC_PRECEDENT_INFLEXION  > 0) { 
3062       STATIC_PRECEDENT_INFLEXION -- ;
3063       return IntWalk_OK;
3064     }
3065   }
3066
3067   //==================================================================================
3068   //=========  D e t e c t    c o n f u s e d    P o in t s       ===========
3069   //==================================================================================
3070
3071   const Standard_Real aSqDist = previousPoint.Value().
3072     SquareDistance(CurrentPoint.Value());
3073
3074
3075   if (aSqDist < Precision::SquareConfusion()) { 
3076     pasInit[0] = Max(pasInit[0], 5.0*ResoU1);
3077     pasInit[1] = Max(pasInit[1], 5.0*ResoV1);
3078     pasInit[2] = Max(pasInit[2], 5.0*ResoU2);
3079     pasInit[3] = Max(pasInit[3], 5.0*ResoV2);
3080
3081     for(Standard_Integer i = 0; i < 4; i++)
3082     {
3083       pasuv[i] = Max(pasuv[i], Min(1.5*pasuv[i], pasInit[i]));
3084     }
3085     //Compute local resolution: for OCC26717
3086     if (Abs(pasuv[choixIso] - pasInit[choixIso]) <= Precision::Confusion())
3087     {
3088       Standard_Real CurU, CurV;
3089       if (choixIso == IntImp_UIsoparametricOnCaro1 ||
3090           choixIso == IntImp_VIsoparametricOnCaro1)
3091         previousPoint.ParametersOnS1(CurU, CurV);
3092       else
3093         previousPoint.ParametersOnS2(CurU, CurV);
3094       gp_Pnt CurPnt = (choixIso == IntImp_UIsoparametricOnCaro1 ||
3095                        choixIso == IntImp_VIsoparametricOnCaro1)?
3096         Adaptor3d_HSurfaceTool::Value(Caro1, CurU, CurV) :
3097         Adaptor3d_HSurfaceTool::Value(Caro2, CurU, CurV);
3098       gp_Pnt OffsetPnt;
3099       switch(choixIso)
3100       {
3101       case IntImp_UIsoparametricOnCaro1:
3102         OffsetPnt =
3103           Adaptor3d_HSurfaceTool::Value(Caro1,
3104                                         CurU + sensCheminement*pasuv[0],
3105                                         CurV);
3106         break;
3107       case IntImp_VIsoparametricOnCaro1:
3108         OffsetPnt =
3109           Adaptor3d_HSurfaceTool::Value(Caro1,
3110                                         CurU,
3111                                         CurV + sensCheminement*pasuv[1]);
3112         break;
3113       case IntImp_UIsoparametricOnCaro2:
3114         OffsetPnt =
3115           Adaptor3d_HSurfaceTool::Value(Caro2,
3116                                         CurU + sensCheminement*pasuv[2],
3117                                         CurV);
3118         break;
3119       case IntImp_VIsoparametricOnCaro2:
3120         OffsetPnt =
3121           Adaptor3d_HSurfaceTool::Value(Caro2,
3122                                         CurU,
3123                                         CurV + sensCheminement*pasuv[3]);
3124         break;
3125       default:break;
3126       }
3127       Standard_Real RefDist = CurPnt.Distance(OffsetPnt);
3128       Standard_Real LocalResol = 0.;
3129       if (RefDist > gp::Resolution())
3130         LocalResol = pasuv[choixIso] * tolconf / RefDist;
3131       if (pasuv[choixIso] < 2*LocalResol)
3132         pasuv[choixIso] = pasInit[choixIso] = 2*LocalResol;
3133     }
3134     ////////////////////////////////////////
3135     aStatus = IntWalk_PointConfondu;
3136   }
3137
3138   //==================================================================================
3139   Standard_Real Up1,Vp1,Uc1,Vc1,Du1,Dv1,AbsDu1,AbsDu2,AbsDv1,AbsDv2;
3140   Standard_Real Up2,Vp2,Uc2,Vc2,Du2,Dv2;
3141
3142   previousPoint.Parameters(Up1,Vp1,Up2,Vp2);
3143   CurrentPoint.Parameters(Uc1,Vc1,Uc2,Vc2);               
3144
3145   Du1 = Uc1 - Up1;   Dv1 = Vc1 - Vp1;
3146   Du2 = Uc2 - Up2;   Dv2 = Vc2 - Vp2;
3147
3148   AbsDu1 = Abs(Du1);
3149   AbsDu2 = Abs(Du2);
3150   AbsDv1 = Abs(Dv1);
3151   AbsDv2 = Abs(Dv2);
3152   //=================================================================================
3153   //====   S t e p   o f   p  r o g r e s s i o n (between previous and Current)   =======
3154   //=================================================================================
3155   if (   AbsDu1 < ResoU1 && AbsDv1 < ResoV1 
3156     && AbsDu2 < ResoU2 && AbsDv2 < ResoV2) {
3157       pasuv[0] = ResoU1; pasuv[1] = ResoV1; pasuv[2] = ResoU2; pasuv[3] = ResoV2;
3158       return(IntWalk_ArretSurPointPrecedent);
3159   }
3160   //==================================================================================
3161
3162   Standard_Real tolArea = 100.0;
3163   if (ResoU1 < Precision::PConfusion() ||
3164     ResoV1 < Precision::PConfusion() ||
3165     ResoU2 < Precision::PConfusion() ||
3166     ResoV2 < Precision::PConfusion() )
3167     tolArea =  tolArea*2.0;
3168
3169   Standard_Real Cosi1, CosRef1, Ang1, AngRef1, ResoUV1, Duv1, d1, tolCoeff1;   
3170   Standard_Real Cosi2, CosRef2, Ang2, AngRef2, ResoUV2, Duv2, d2, tolCoeff2;   
3171   Cosi1 = Du1*previousd1.X() + Dv1*previousd1.Y();
3172   Cosi2 = Du2*previousd2.X() + Dv2*previousd2.Y();
3173   Duv1 = Du1*Du1 + Dv1*Dv1;
3174   Duv2 = Du2*Du2 + Dv2*Dv2;
3175   ResoUV1 = ResoU1*ResoU1 + ResoV1*ResoV1;
3176   ResoUV2 = ResoU2*ResoU2 + ResoV2*ResoV2;
3177   //
3178   //modified by NIZNHY-PKV Wed Nov 13 12:25:44 2002 f
3179   //
3180   Standard_Real aMinDiv2=Precision::Confusion();
3181   aMinDiv2=aMinDiv2*aMinDiv2;
3182   //
3183   d1=d;
3184   if (Duv1>aMinDiv2)  {
3185     d1 = Abs(ResoUV1/Duv1);
3186     d1 = Min(Sqrt(d1)*tolArea, d);  
3187   } 
3188   //d1 = Abs(ResoUV1/Duv1); 
3189   //d1 = Min(Sqrt(d1)*tolArea,d);  
3190   //modified by NIZNHY-PKV Wed Nov 13 12:34:30 2002 t
3191   tolCoeff1 = Exp(d1);
3192   //
3193   //modified by NIZNHY-PKV Wed Nov 13 12:34:43 2002 f
3194   d2=d;
3195   if (Duv2>aMinDiv2) {
3196     d2 = Abs(ResoUV2/Duv2); 
3197     d2 = Min(Sqrt(d2)*tolArea,d); 
3198   }
3199   //d2 = Abs(ResoUV2/Duv2); 
3200   //d2 = Min(Sqrt(d2)*tolArea,d);  
3201   //modified by NIZNHY-PKV Wed Nov 13 12:34:53 2002 t
3202   tolCoeff2 = Exp(d2);
3203   CosRef1 = CosRef2D/tolCoeff1;
3204   CosRef2 = CosRef2D/tolCoeff2;
3205   //
3206   //==================================================================================
3207   //== The points are not confused :                                           ==
3208   //== 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 ==
3209   //==                           N o t    T o o    G r e a t (angle in space UV)    ==
3210   //==                           C h a n g e    o f    s i d e                ==
3211   //==================================================================================
3212   if (aStatus != IntWalk_PointConfondu) {
3213     if(Cosi1*Cosi1 < CosRef1*Duv1 || Cosi2*Cosi2 < CosRef2*Duv2) {
3214       pasuv[0]*=0.5;  pasuv[1]*=0.5;  pasuv[2]*=0.5;  pasuv[3]*=0.5;
3215       if (pasuv[0]<ResoU1 && pasuv[1]<ResoV1 && pasuv[2]<ResoU2 && pasuv[3]<ResoV2) { 
3216         return(IntWalk_ArretSurPointPrecedent);
3217       }
3218       else {
3219         pasuv[0]*=0.5; pasuv[1]*=0.5; pasuv[2]*=0.5; pasuv[3]*=0.5;
3220         return(IntWalk_PasTropGrand);
3221       }
3222     }
3223     const gp_Dir2d& Tg2dcourante1 = myIntersectionOn2S.DirectionOnS1();
3224     const gp_Dir2d& Tg2dcourante2 = myIntersectionOn2S.DirectionOnS2();
3225     Cosi1 = Du1*Tg2dcourante1.X() + Dv1*Tg2dcourante1.Y();
3226     Cosi2 = Du2*Tg2dcourante2.X() + Dv2*Tg2dcourante2.Y();
3227     Ang1 = Abs(previousd1.Angle(Tg2dcourante1));  
3228     Ang2 = Abs(previousd2.Angle(Tg2dcourante2));  
3229     AngRef1 = AngRef2D*tolCoeff1;
3230     AngRef2 = AngRef2D*tolCoeff2;
3231     //-------------------------------------------------------
3232     //-- Test : Angle too great in space UV       -----
3233     //--        Change of  side                      -----
3234     //-------------------------------------------------------
3235     if(Cosi1*Cosi1 < CosRef1*Duv1 || Cosi2*Cosi2 < CosRef2*Duv2 || Ang1 > AngRef1 || Ang2 > AngRef2) {
3236       pasuv[0]*=0.5;  pasuv[1]*=0.5;  pasuv[2]*=0.5;  pasuv[3]*=0.5;
3237       if (pasuv[0]<ResoU1 && pasuv[1]<ResoV1 && pasuv[2]<ResoU2 && pasuv[3]<ResoV2) 
3238         return(IntWalk_ArretSurPoint);
3239       else 
3240         return(IntWalk_PasTropGrand);
3241     }
3242   }
3243   //<-OCC431(apo)
3244   //==================================================================================
3245   //== D e t e c t i o n   o f    :  Step Too Small 
3246   //==                               STEP TOO Great 
3247   //==================================================================================
3248
3249   //---------------------------------------
3250   //-- Estimate of the vector           --
3251   //---------------------------------------
3252   FlecheCourante =
3253     Sqrt(Abs((previousd.XYZ()-TgCourante.XYZ()).SquareModulus()*aSqDist))/8.;
3254
3255   if ( FlecheCourante<= fleche*0.5) {     //-- Current step too small
3256     if(FlecheCourante>1e-16) { 
3257       Ratio = 0.5*(fleche/FlecheCourante);
3258     }
3259     else { 
3260       Ratio = 10.0;
3261     }
3262     Standard_Real pasSu1 = pasuv[0];
3263     Standard_Real pasSv1 = pasuv[1];
3264     Standard_Real pasSu2 = pasuv[2];
3265     Standard_Real pasSv2 = pasuv[3];
3266
3267     //-- In  case if 
3268     //-- a point at U+DeltaU is required, ....
3269     //-- return a point at U + Epsilon
3270     //-- Epsilon << DeltaU.
3271
3272     if(pasuv[0]< AbsDu1) pasuv[0] = AbsDu1;
3273     if(pasuv[1]< AbsDv1) pasuv[1] = AbsDv1;
3274     if(pasuv[2]< AbsDu2) pasuv[2] = AbsDu2;
3275     if(pasuv[3]< AbsDv2) pasuv[3] = AbsDv2;
3276
3277     if(pasuv[0]<ResoU1) pasuv[0]=ResoU1;
3278     if(pasuv[1]<ResoV1) pasuv[1]=ResoV1;
3279     if(pasuv[2]<ResoU2) pasuv[2]=ResoU2;
3280     if(pasuv[3]<ResoV2) pasuv[3]=ResoV2;
3281     //-- if(Ratio>10.0 ) { Ratio=10.0; } 
3282     Standard_Real R1,R = pasInit[0]/pasuv[0];
3283     R1= pasInit[1]/pasuv[1];     if(R1<R) R=R1;
3284     R1= pasInit[2]/pasuv[2];     if(R1<R) R=R1;
3285     R1= pasInit[3]/pasuv[3];     if(R1<R) R=R1;
3286     if(Ratio > R) Ratio=R;
3287     pasuv[0] = Min(Ratio*pasuv[0],pasInit[0]);
3288     pasuv[1] = Min(Ratio*pasuv[1],pasInit[1]);
3289     pasuv[2] = Min(Ratio*pasuv[2],pasInit[2]);
3290     pasuv[3] = Min(Ratio*pasuv[3],pasInit[3]);
3291     if (pasuv[0] != pasSu1 || pasuv[2] != pasSu2|| 
3292       pasuv[1] != pasSv1 || pasuv[3] != pasSv2) {
3293         if(++STATIC_BLOCAGE_SUR_PAS_TROP_GRAND > 5) {
3294           STATIC_BLOCAGE_SUR_PAS_TROP_GRAND = 0;
3295           return IntWalk_PasTropGrand; 
3296         }
3297     }
3298     if(aStatus == IntWalk_OK) {
3299       STATIC_BLOCAGE_SUR_PAS_TROP_GRAND=0;
3300       //-- Try to increase the step
3301     }
3302     return aStatus;
3303   }
3304   else {                                //-- CurrentVector > vector*0.5 
3305     if (FlecheCourante > fleche) {      //-- Current step too Great
3306       Ratio = fleche/FlecheCourante; 
3307       pasuv[0] = Ratio*pasuv[0];
3308       pasuv[1] = Ratio*pasuv[1];
3309       pasuv[2] = Ratio*pasuv[2];
3310       pasuv[3] = Ratio*pasuv[3];
3311       //if(++STATIC_BLOCAGE_SUR_PAS_TROP_GRAND > 5) {
3312       //        STATIC_BLOCAGE_SUR_PAS_TROP_GRAND = 0;
3313       return IntWalk_PasTropGrand; 
3314       //}
3315     }
3316     else {                             //-- vector/2  <  CurrentVector <= vector   
3317       Ratio = 0.75 * (fleche / FlecheCourante);
3318     }
3319   }
3320
3321   if(aStatus != IntWalk_PointConfondu)
3322   {
3323     //Here, aCosBetweenTangent >= 0.0 definitely.
3324
3325     /*
3326     Brief algorithm description.
3327     We have two (not-coincindent) intersection point (P1 and P2). In every point,
3328     vector of tangent (to the intersection curve) is known (vectors T1 and T2).
3329     Basing on these data, we create osculating circle.
3330
3331                                     * - arc of osculating circle
3332                                 *      * 
3333                            P1 x----------x P2
3334                              /            \
3335                             /              \
3336                           Vec(T1)         Vec(T2)
3337
3338     Let me draw your attention to the following facts:
3339       1. Vectors T1 and T2 direct FROM (not TO) points P1 and P2. Therefore,
3340         one of previously computed vector should be reversed.
3341
3342     In this case, the absolute (!) value of the deflection between the arc of
3343     the osculating circle and the P1P2 segment can be computed as follows:
3344           e = d*(1-sin(B/2))/(2*cos(B/2)),      (1)
3345     where d is the length of P1P2 segment, B is the angle between vectors T1 and T2.
3346     At that,
3347               pi/2 <= B <= pi,
3348               cos(B/2) >= 0,
3349               sin(B/2) > 0,
3350               sin(B) > 0,
3351               cos(B) < 0.
3352
3353     Later, the normal state of algorithm work is (as we apply) 
3354               tolconf/2 <= e <= tolconf.
3355     In this case, we keep previous step.
3356
3357     If e < tolconf/2 then the local curvature of the intersection curve is small.
3358     As result, the step should be increased.
3359
3360     If e > tolconf then the step is too big. Therefore, we should decrease one.
3361
3362     Condition (1) is equivalent to
3363               sin(B/2) = 1 - 2/(1+(d/(2*e))^2) = Fs(e),
3364               cos(B) = 1 - 2*Fs(e)^2 = Fd(e),
3365     where Fs(e)and Fd(e) are some function with parameter "deflection".
3366
3367     Let mean that Fs(e) is decreasing function. Fd(e) is increasing function,