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