0026938: Boolean operations fail between two ellipsoids
[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 Standard_Boolean IntWalk_PWalking::
2008 DistanceMinimizeByGradient( const Handle(Adaptor3d_HSurface)& theASurf1,
2009                            const Handle(Adaptor3d_HSurface)& theASurf2,
2010                            Standard_Real& theU1,
2011                            Standard_Real& theV1,
2012                            Standard_Real& theU2,
2013                            Standard_Real& theV2,
2014                            const Standard_Real theStep0U1V1,
2015                            const Standard_Real theStep0U2V2)
2016 {
2017   const Standard_Integer aNbIterMAX = 60;
2018   const Standard_Real aTol = 1.0e-14;
2019   Handle(Geom_Surface) aS1, aS2;
2020
2021   if (theASurf1->GetType() != GeomAbs_BezierSurface &&
2022       theASurf1->GetType() != GeomAbs_BSplineSurface)
2023       return Standard_True;
2024   if (theASurf2->GetType() != GeomAbs_BezierSurface &&
2025       theASurf2->GetType() != GeomAbs_BSplineSurface)
2026       return Standard_True;
2027
2028   Standard_Boolean aStatus = Standard_False;
2029
2030   gp_Pnt aP1, aP2;
2031   gp_Vec aD1u, aD1v, aD2U, aD2V;
2032
2033   theASurf1->D1(theU1, theV1, aP1, aD1u, aD1v);
2034   theASurf2->D1(theU2, theV2, aP2, aD2U, aD2V);
2035
2036   Standard_Real aSQDistPrev = aP1.SquareDistance(aP2);
2037
2038   gp_Vec aP12(aP1, aP2);
2039
2040   Standard_Real aGradFu(-aP12.Dot(aD1u));
2041   Standard_Real aGradFv(-aP12.Dot(aD1v));
2042   Standard_Real aGradFU( aP12.Dot(aD2U));
2043   Standard_Real aGradFV( aP12.Dot(aD2V));
2044
2045   Standard_Real aSTEPuv = theStep0U1V1, aStepUV = theStep0U2V2;
2046
2047   Standard_Boolean flRepeat = Standard_True;
2048   Standard_Integer aNbIter = aNbIterMAX;
2049
2050   while(flRepeat)
2051   {
2052     Standard_Real anAdd = aGradFu*aSTEPuv;
2053     Standard_Real aPARu = (anAdd >= 0.0)?
2054       (theU1 - Max(anAdd, Epsilon(theU1))) :
2055     (theU1 + Max(-anAdd, Epsilon(theU1)));
2056     anAdd = aGradFv*aSTEPuv;
2057     Standard_Real aPARv = (anAdd >= 0.0)?
2058       (theV1 - Max(anAdd, Epsilon(theV1))) :
2059     (theV1 + Max(-anAdd, Epsilon(theV1)));
2060     anAdd = aGradFU*aStepUV;
2061     Standard_Real aParU = (anAdd >= 0.0)?
2062       (theU2 - Max(anAdd, Epsilon(theU2))) :
2063     (theU2 + Max(-anAdd, Epsilon(theU2)));
2064     anAdd = aGradFV*aStepUV;
2065     Standard_Real aParV = (anAdd >= 0.0)?
2066       (theV2 - Max(anAdd, Epsilon(theV2))) :
2067     (theV2 + Max(-anAdd, Epsilon(theV2)));
2068
2069     gp_Pnt aPt1, aPt2;
2070
2071     theASurf1->D1(aPARu, aPARv, aPt1, aD1u, aD1v);
2072     theASurf2->D1(aParU, aParV, aPt2, aD2U, aD2V);
2073
2074     Standard_Real aSQDist = aPt1.SquareDistance(aPt2);
2075
2076     if(aSQDist < aSQDistPrev)
2077     {
2078       aSQDistPrev = aSQDist;
2079       theU1 = aPARu;
2080       theV1 = aPARv;
2081       theU2 = aParU;
2082       theV2 = aParV;
2083
2084       aStatus = aSQDistPrev < aTol;
2085       aSTEPuv *= 1.2;
2086       aStepUV *= 1.2;
2087     }
2088     else
2089     {
2090       if(--aNbIter < 0)
2091       {
2092         flRepeat = Standard_False;
2093       }
2094       else
2095       {
2096         theASurf1->D1(theU1, theV1, aPt1, aD1u, aD1v);
2097         theASurf2->D1(theU2, theV2, aPt2, aD2U, aD2V);
2098
2099         gp_Vec aPt12(aPt1, aPt2);
2100         aGradFu = -aPt12.Dot(aD1u);
2101         aGradFv = -aPt12.Dot(aD1v);
2102         aGradFU = aPt12.Dot(aD2U);
2103         aGradFV = aPt12.Dot(aD2V);
2104         aSTEPuv = theStep0U1V1;
2105         aStepUV = theStep0U2V2;
2106       }
2107     }
2108   }
2109
2110   return aStatus;
2111 }
2112
2113 //=======================================================================
2114 //function : DistanceMinimizeByExtrema
2115 //purpose  : 
2116 //=======================================================================
2117 Standard_Boolean IntWalk_PWalking::
2118 DistanceMinimizeByExtrema(const Handle(Adaptor3d_HSurface)& theASurf, 
2119                           const gp_Pnt& theP0,
2120                           Standard_Real& theU0,
2121                           Standard_Real& theV0,
2122                           const Standard_Real theStep0U,
2123                           const Standard_Real theStep0V)
2124 {
2125   const Standard_Real aTol = 1.0e-14;
2126   gp_Pnt aPS;
2127   gp_Vec aD1Su, aD1Sv, aD2Su, aD2Sv, aD2SuvTemp;
2128   Standard_Real aSQDistPrev = RealLast();
2129   Standard_Real aU = theU0, aV = theV0;
2130
2131   Standard_Integer aNbIter = 10;
2132   do
2133   {
2134     theASurf->D2(aU, aV, aPS, aD1Su, aD1Sv, aD2Su, aD2Sv, aD2SuvTemp);
2135
2136     gp_Vec aVec(theP0, aPS);
2137
2138     Standard_Real aSQDist = aVec.SquareMagnitude();
2139
2140     if(aSQDist >= aSQDistPrev)
2141       break;
2142
2143     aSQDistPrev = aSQDist;
2144     theU0 = aU;
2145     theV0 = aV;
2146     aNbIter--;
2147
2148     if(aSQDistPrev < aTol)
2149       break;
2150
2151     //Functions
2152     const Standard_Real aF1 = aD1Su.Dot(aVec), aF2 = aD1Sv.Dot(aVec);
2153
2154     //Derivatives
2155     const Standard_Real aDf1u = aD2Su.Dot(aVec) + aD1Su.Dot(aD1Su),
2156       aDf1v = aD2Su.Dot(aD1Sv),
2157       aDf2u = aDf1v,
2158       aDf2v = aD2Sv.Dot(aVec) + aD1Sv.Dot(aD1Sv);
2159
2160     const Standard_Real aDet = aDf1u*aDf2v - aDf1v*aDf2u;
2161     aU -= theStep0U*(aDf2v*aF1 - aDf1v*aF2)/aDet;
2162     aV += theStep0V*(aDf2u*aF1 - aDf1u*aF2)/aDet;
2163   }
2164   while(aNbIter > 0);
2165
2166   return (aSQDistPrev < aTol);
2167 }
2168
2169 //=======================================================================
2170 //function : HandleSingleSingularPoint
2171 //purpose  : 
2172 //=======================================================================
2173 Standard_Boolean IntWalk_PWalking::HandleSingleSingularPoint(const Handle(Adaptor3d_HSurface)& theASurf1,
2174                                                              const Handle(Adaptor3d_HSurface)& theASurf2,
2175                                                              const Standard_Real the3DTol,
2176                                                              TColStd_Array1OfReal &thePnt)
2177 {
2178   // u1, v1, u2, v2 order is used.
2179   Standard_Real aLowBorder[4] = {theASurf1->FirstUParameter(),
2180                                  theASurf1->FirstVParameter(),
2181                                  theASurf2->FirstUParameter(),
2182                                  theASurf2->FirstVParameter()};
2183   Standard_Real aUppBorder[4] = {theASurf1->LastUParameter(),
2184                                  theASurf1->LastVParameter(),
2185                                  theASurf2->LastUParameter(),
2186                                  theASurf2->LastVParameter()};
2187   IntImp_ConstIsoparametric aLockedDir[4] = {IntImp_UIsoparametricOnCaro1,
2188                                              IntImp_VIsoparametricOnCaro1,
2189                                              IntImp_UIsoparametricOnCaro2,
2190                                              IntImp_VIsoparametricOnCaro2};
2191
2192   // Create new intersector with new tolerance.
2193   IntWalk_TheInt2S anInt(theASurf1, theASurf2, the3DTol);
2194   math_FunctionSetRoot aRsnld(anInt.Function());
2195
2196   for (Standard_Integer i = 1; i <= 4; ++i)
2197   {
2198     if ( Abs(thePnt(i) - aLowBorder[i - 1]) < Precision::PConfusion() ||
2199          Abs(thePnt(i) - aUppBorder[i - 1]) < Precision::PConfusion())
2200     {
2201
2202       anInt.Perform(thePnt,aRsnld, aLockedDir[i - 1]);
2203
2204       if (!anInt.IsDone())
2205         continue;
2206
2207       if (anInt.IsEmpty())
2208         continue;
2209
2210       anInt.Point().Parameters(thePnt(1), thePnt(2), thePnt(3), thePnt(4));
2211       return Standard_True;
2212     }
2213   }
2214
2215   return Standard_False;
2216 }
2217
2218 //=======================================================================
2219 //function : SeekPointOnBoundary
2220 //purpose  : 
2221 //=======================================================================
2222 Standard_Boolean IntWalk_PWalking::
2223 SeekPointOnBoundary(const Handle(Adaptor3d_HSurface)& theASurf1,
2224                     const Handle(Adaptor3d_HSurface)& theASurf2,
2225                     const Standard_Real theU1,
2226                     const Standard_Real theV1,
2227                     const Standard_Real theU2,
2228                     const Standard_Real theV2,
2229                     const Standard_Boolean isTheFirst)
2230 {
2231   Standard_Boolean isOK = Standard_False;
2232
2233   // Tune solution tolerance according with object size.
2234   const Standard_Real aRes1 = Max(Precision::PConfusion() / theASurf1->UResolution(1.0),
2235                                   Precision::PConfusion() / theASurf1->VResolution(1.0));
2236   const Standard_Real aRes2 = Max(Precision::PConfusion() / theASurf2->UResolution(1.0),
2237                                   Precision::PConfusion() / theASurf2->VResolution(1.0));
2238   const Standard_Real a3DTol = Max(aRes1, aRes2);
2239   const Standard_Real aTol = Max(Precision::Confusion(), a3DTol);
2240
2241   // u1, v1, u2, v2 order is used.
2242   TColStd_Array1OfReal aPnt(1,4);
2243   aPnt(1) = theU1; aPnt(2) = theV1; aPnt(3) = theU2; aPnt(4) = theV2;
2244   TColStd_Array1OfReal aSingularPnt(aPnt);
2245
2246   Standard_Integer aNbIter = 20;
2247   Standard_Boolean aStatus = Standard_False;
2248   do
2249   {
2250     aNbIter--;
2251     aStatus = DistanceMinimizeByGradient(theASurf1, theASurf2, aPnt(1), aPnt(2), aPnt(3), aPnt(4));
2252     if(aStatus)
2253       break;
2254
2255     aStatus = DistanceMinimizeByExtrema(theASurf1, theASurf2->Value(aPnt(3), aPnt(4)), aPnt(1), aPnt(2));
2256     if(aStatus)
2257       break;
2258
2259     aStatus = DistanceMinimizeByExtrema(theASurf2, theASurf1->Value(aPnt(1), aPnt(2)), aPnt(3), aPnt(4));
2260     if(aStatus)
2261       break;
2262   }
2263   while(!aStatus && (aNbIter > 0));
2264
2265   // Handle singular points.
2266   Standard_Boolean aSingularStatus = HandleSingleSingularPoint(theASurf1, theASurf2, aTol, aSingularPnt);
2267   if (aSingularStatus)
2268     aPnt = aSingularPnt;
2269
2270   if(aStatus || aSingularStatus)
2271   {
2272     gp_Pnt  aP1 = theASurf1->Value(aPnt(1), aPnt(2)),
2273             aP2 = theASurf2->Value(aPnt(3), aPnt(4));
2274     gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));
2275
2276     const Standard_Real aSQDist = aPInt.SquareDistance(aP1);
2277     if (aSQDist < aTol * aTol)
2278     {
2279       IntSurf_PntOn2S anIP;
2280       anIP.SetValue(aPInt, aPnt(1), aPnt(2), aPnt(3), aPnt(4));
2281
2282       if(isTheFirst)
2283         line->InsertBefore(1,anIP);
2284       else
2285         line->Add(anIP);
2286
2287       isOK = Standard_True;
2288     }
2289   }
2290
2291   return isOK;
2292 }
2293
2294 //=======================================================================
2295 //function : PutToBoundary
2296 //purpose  : 
2297 //=======================================================================
2298 Standard_Boolean IntWalk_PWalking::
2299 PutToBoundary(const Handle(Adaptor3d_HSurface)& theASurf1,
2300               const Handle(Adaptor3d_HSurface)& theASurf2)
2301 {
2302   const Standard_Real aTolMin = Precision::Confusion();
2303
2304   Standard_Boolean hasBeenAdded = Standard_False;
2305
2306   const Standard_Real aU1bFirst = theASurf1->FirstUParameter();
2307   const Standard_Real aU1bLast = theASurf1->LastUParameter();
2308   const Standard_Real aU2bFirst = theASurf2->FirstUParameter();
2309   const Standard_Real aU2bLast = theASurf2->LastUParameter();
2310   const Standard_Real aV1bFirst = theASurf1->FirstVParameter();
2311   const Standard_Real aV1bLast = theASurf1->LastVParameter();
2312   const Standard_Real aV2bFirst = theASurf2->FirstVParameter();
2313   const Standard_Real aV2bLast = theASurf2->LastVParameter();
2314
2315   Standard_Real aTol = 1.0;
2316   aTol = Min(aTol, aU1bLast - aU1bFirst);
2317   aTol = Min(aTol, aU2bLast - aU2bFirst);
2318   aTol = Min(aTol, aV1bLast - aV1bFirst);
2319   aTol = Min(aTol, aV2bLast - aV2bFirst)*1.0e-3;
2320
2321   if(aTol <= 2.0*aTolMin)
2322     return hasBeenAdded;
2323
2324   Standard_Boolean isNeedAdding = Standard_False;
2325   Standard_Boolean isU1parallel = Standard_False, isV1parallel = Standard_False;
2326   Standard_Boolean isU2parallel = Standard_False, isV2parallel = Standard_False;
2327   IsParallel(line, Standard_True, aTol, isU1parallel, isV1parallel);
2328   IsParallel(line, Standard_False, aTol, isU2parallel, isV2parallel);
2329
2330   Standard_Real u1, v1, u2, v2;
2331   line->Value(1).Parameters(u1, v1, u2, v2);
2332   Standard_Real aDelta = 0.0;
2333
2334   if(!isV1parallel)
2335   {
2336     aDelta = u1 - aU1bFirst;
2337     if((aTolMin < aDelta) && (aDelta < aTol))
2338     {
2339       u1 = aU1bFirst;
2340       isNeedAdding = Standard_True;
2341     }
2342     else
2343     {
2344       aDelta = aU1bLast - u1;
2345       if((aTolMin < aDelta) && (aDelta < aTol))
2346       {
2347         u1 = aU1bLast;
2348         isNeedAdding = Standard_True;
2349       }
2350     }
2351   }
2352
2353   if(!isV2parallel)
2354   {
2355     aDelta = u2 - aU2bFirst;
2356     if((aTolMin < aDelta) && (aDelta < aTol))
2357     {
2358       u2 = aU2bFirst;
2359       isNeedAdding = Standard_True;
2360     }
2361     else
2362     {
2363       aDelta = aU2bLast - u2;
2364       if((aTolMin < aDelta) && (aDelta < aTol))
2365       {
2366         u2 = aU2bLast;
2367         isNeedAdding = Standard_True;
2368       }
2369     }
2370   }
2371
2372   if(!isU1parallel)
2373   {
2374     aDelta = v1 - aV1bFirst;
2375     if((aTolMin < aDelta) && (aDelta < aTol))
2376     {
2377       v1 = aV1bFirst;
2378       isNeedAdding = Standard_True;
2379     }
2380     else
2381     {
2382       aDelta = aV1bLast - v1;
2383       if((aTolMin < aDelta) && (aDelta < aTol))
2384       {
2385         v1 = aV1bLast;
2386         isNeedAdding = Standard_True;
2387       }
2388     }
2389   }
2390
2391   if(!isU2parallel)
2392   {
2393     aDelta = v2 - aV2bFirst;
2394     if((aTolMin < aDelta) && (aDelta < aTol))
2395     {
2396       v2 = aV2bFirst;
2397       isNeedAdding = Standard_True;
2398     }
2399     else
2400     {
2401       aDelta = aV2bLast - v2;
2402       if((aTolMin < aDelta) && (aDelta < aTol))
2403       {
2404         v2 = aV2bLast;
2405         isNeedAdding = Standard_True;
2406       }
2407     }
2408   }
2409
2410   if(isNeedAdding)
2411   {
2412     hasBeenAdded = 
2413       SeekPointOnBoundary(theASurf1, theASurf2, u1, 
2414       v1, u2, v2, Standard_True);
2415   }
2416
2417   const Standard_Integer aNbPnts = line->NbPoints();
2418   isNeedAdding = Standard_False;
2419   line->Value(aNbPnts).Parameters(u1, v1, u2, v2);
2420
2421   if(!isV1parallel)
2422   {
2423     aDelta = u1 - aU1bFirst;
2424     if((aTolMin < aDelta) && (aDelta < aTol))
2425     {
2426       u1 = aU1bFirst;
2427       isNeedAdding = Standard_True;
2428     }
2429     else
2430     {
2431       aDelta = aU1bLast - u1;
2432       if((aTolMin < aDelta) && (aDelta < aTol))
2433       {
2434         u1 = aU1bLast;
2435         isNeedAdding = Standard_True;
2436       }
2437     }
2438   }
2439
2440   if(!isV2parallel)
2441   {
2442     aDelta = u2 - aU2bFirst;
2443     if((aTolMin < aDelta) && (aDelta < aTol))
2444     {
2445       u2 = aU2bFirst;
2446       isNeedAdding = Standard_True;
2447     }
2448     else
2449     {
2450       aDelta = aU2bLast - u2;
2451       if((aTolMin < aDelta) && (aDelta < aTol))
2452       {
2453         u2 = aU2bLast;
2454         isNeedAdding = Standard_True;
2455       }
2456     }
2457   }
2458
2459   if(!isU1parallel)
2460   {
2461     aDelta = v1 - aV1bFirst;
2462     if((aTolMin < aDelta) && (aDelta < aTol))
2463     {
2464       v1 = aV1bFirst;
2465       isNeedAdding = Standard_True;
2466     }
2467     else
2468     {
2469       aDelta = aV1bLast - v1;
2470       if((aTolMin < aDelta) && (aDelta < aTol))
2471       {
2472         v1 = aV1bLast;
2473         isNeedAdding = Standard_True;
2474       }
2475     }
2476   }
2477
2478   if(!isU2parallel)
2479   {
2480     aDelta = v2 - aV2bFirst;
2481     if((aTolMin < aDelta) && (aDelta < aTol))
2482     {
2483       v2 = aV2bFirst;
2484       isNeedAdding = Standard_True;
2485     }
2486     else
2487     {
2488       aDelta = aV2bLast - v2;
2489       if((aTolMin < aDelta) && (aDelta < aTol))
2490       {
2491         v2 = aV2bLast;
2492         isNeedAdding = Standard_True;
2493       }
2494     }
2495   }
2496
2497   if(isNeedAdding)
2498   {
2499     hasBeenAdded = 
2500       SeekPointOnBoundary(theASurf1, theASurf2, u1, 
2501       v1, u2, v2, Standard_False);
2502   }
2503
2504   return hasBeenAdded;
2505 }
2506
2507 //=======================================================================
2508 //function : SeekAdditionalPoints
2509 //purpose  : 
2510 //=======================================================================
2511 Standard_Boolean IntWalk_PWalking::
2512 SeekAdditionalPoints( const Handle(Adaptor3d_HSurface)& theASurf1,
2513                      const Handle(Adaptor3d_HSurface)& theASurf2,
2514                      const Standard_Integer theMinNbPoints)
2515 {
2516   const Standard_Real aTol = 1.0e-14;
2517   Standard_Integer aNbPoints = line->NbPoints();
2518   if(aNbPoints > theMinNbPoints)
2519     return Standard_True;
2520
2521   const Standard_Real aU1bFirst = theASurf1->FirstUParameter();
2522   const Standard_Real aU1bLast = theASurf1->LastUParameter();
2523   const Standard_Real aU2bFirst = theASurf2->FirstUParameter();
2524   const Standard_Real aU2bLast = theASurf2->LastUParameter();
2525   const Standard_Real aV1bFirst = theASurf1->FirstVParameter();
2526   const Standard_Real aV1bLast = theASurf1->LastVParameter();
2527   const Standard_Real aV2bFirst = theASurf2->FirstVParameter();
2528   const Standard_Real aV2bLast = theASurf2->LastVParameter();
2529
2530
2531   Standard_Boolean isPrecise = Standard_False;
2532
2533   Standard_Real U1prec = 0.0, V1prec = 0.0, U2prec = 0.0, V2prec = 0.0;
2534
2535   Standard_Integer aNbPointsPrev = 0;
2536   while(aNbPoints < theMinNbPoints && (aNbPoints != aNbPointsPrev))
2537   {
2538     aNbPointsPrev = aNbPoints;
2539     for(Standard_Integer fp = 1, lp = 2; fp < aNbPoints; fp = lp + 1)
2540     {
2541       Standard_Real U1f, V1f, U2f, V2f; //first point in 1st and 2nd surafaces
2542       Standard_Real U1l, V1l, U2l, V2l; //last  point in 1st and 2nd surafaces
2543
2544       lp = fp+1;
2545       line->Value(fp).Parameters(U1f, V1f, U2f, V2f);
2546       line->Value(lp).Parameters(U1l, V1l, U2l, V2l);
2547
2548       U1prec = 0.5*(U1f+U1l);
2549       if(U1prec < aU1bFirst)
2550         U1prec = aU1bFirst;
2551       if(U1prec > aU1bLast)
2552         U1prec = aU1bLast;
2553
2554       V1prec = 0.5*(V1f+V1l);
2555       if(V1prec < aV1bFirst)
2556         V1prec = aV1bFirst;
2557       if(V1prec > aV1bLast)
2558         V1prec = aV1bLast;
2559
2560       U2prec = 0.5*(U2f+U2l);
2561       if(U2prec < aU2bFirst)
2562         U2prec = aU2bFirst;
2563       if(U2prec > aU2bLast)
2564         U2prec = aU2bLast;
2565
2566       V2prec = 0.5*(V2f+V2l);
2567       if(V2prec < aV2bFirst)
2568         V2prec = aV2bFirst;
2569       if(V2prec > aV2bLast)
2570         V2prec = aV2bLast;
2571
2572       Standard_Boolean aStatus = Standard_False;
2573       Standard_Integer aNbIter = 5;
2574       do
2575       {
2576         aStatus = DistanceMinimizeByGradient(theASurf1, theASurf2, U1prec, V1prec, U2prec, V2prec);
2577         if(aStatus)
2578         {
2579           break;
2580         }
2581
2582         aStatus = DistanceMinimizeByExtrema(theASurf1, theASurf2->Value(U2prec, V2prec), U1prec, V1prec);
2583         if(aStatus)
2584         {
2585           break;
2586         }
2587
2588         aStatus = DistanceMinimizeByExtrema(theASurf2, theASurf1->Value(U1prec, V1prec), U2prec, V2prec);
2589         if(aStatus)
2590         {
2591           break;
2592         }
2593       }
2594       while(!aStatus && (--aNbIter > 0));
2595
2596       if(aStatus)
2597       {
2598         gp_Pnt  aP1 = theASurf1->Value(U1prec, V1prec),
2599           aP2 = theASurf2->Value(U2prec, V2prec);
2600         gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));
2601
2602         const Standard_Real aSQDist1 = aPInt.SquareDistance(aP1),
2603           aSQDist2 = aPInt.SquareDistance(aP2);
2604
2605         if((aSQDist1 < aTol) && (aSQDist2 < aTol))
2606         {
2607           IntSurf_PntOn2S anIP;
2608           anIP.SetValue(aPInt, U1prec, V1prec, U2prec, V2prec);
2609           line->InsertBefore(lp, anIP);
2610
2611           isPrecise = Standard_True;
2612
2613           if(++aNbPoints >= theMinNbPoints)
2614             break;
2615         }
2616         else
2617         {
2618           lp--;
2619         }
2620       }
2621     }
2622   }
2623
2624   return isPrecise;
2625 }
2626
2627 void IntWalk_PWalking::
2628 RepartirOuDiviser(Standard_Boolean& DejaReparti,
2629                   IntImp_ConstIsoparametric& ChoixIso,
2630                   Standard_Boolean& Arrive) 
2631
2632                   // at the neighborhood of a point, there is a fail of marching 
2633                   // it is required to divide the steps to try to continue
2634                   // if the step is too small if we are on border
2635                   // restart in another direction if it was not done, otherwise stop
2636
2637 {
2638   //  Standard_Integer i;
2639   if (Arrive) {    //restart in the other direction
2640     if (!DejaReparti ) {
2641       Arrive        = Standard_False; 
2642       DejaReparti   = Standard_True;
2643       previousPoint = line->Value(1);
2644       previoustg    = Standard_False;
2645       previousd1    = firstd1;
2646       previousd2    = firstd2;
2647       previousd     = tgdir;
2648       indextg       = line->NbPoints();
2649       tgdir.Reverse();
2650       line->Reverse();
2651
2652       //-- printf("\nIntWalk_PWalking_2.gxx Reverse %3d\n",indextg);
2653       sensCheminement = -1;
2654       tgfirst      = tglast;
2655       tglast       = Standard_False;
2656       ChoixIso     = choixIsoSav;
2657 #if 0
2658       pasuv[0]=pasSav[0];
2659       pasuv[1]=pasSav[1];
2660       pasuv[2]=pasSav[2];
2661       pasuv[3]=pasSav[3];
2662 #else 
2663       Standard_Real u1,v1,u2,v2;
2664       Standard_Real U1,V1,U2,V2;
2665       Standard_Integer nn=line->NbPoints();
2666       if(nn>2) { 
2667         line->Value(nn).Parameters(u1,v1,u2,v2);
2668         line->Value(nn-1).Parameters(U1,V1,U2,V2);
2669         pasuv[0]=Abs(u1-U1);
2670         pasuv[1]=Abs(v1-V1);
2671         pasuv[2]=Abs(u2-U2);
2672         pasuv[3]=Abs(v2-V2);
2673       }
2674 #endif
2675
2676     }
2677   }  
2678   else  {
2679     if (    pasuv[0]*0.5 < ResoU1
2680       &&  pasuv[1]*0.5 < ResoV1
2681       &&  pasuv[2]*0.5 < ResoU2
2682       &&  pasuv[3]*0.5 < ResoV2
2683       ) {
2684         if (!previoustg) {
2685           tglast = Standard_True;      // IS IT ENOUGH ????
2686         }
2687
2688         if (!DejaReparti) {  //restart in the other direction
2689           DejaReparti       = Standard_True;
2690           previousPoint     = line->Value(1);
2691           previoustg        = Standard_False;
2692           previousd1        = firstd1;
2693           previousd2        = firstd2;
2694           previousd         = tgdir;
2695           indextg           = line->NbPoints();
2696           tgdir.Reverse();
2697           line->Reverse();
2698
2699           //-- printf("\nIntWalk_PWalking_2.gxx Reverse %3d\n",indextg);
2700
2701           sensCheminement   = -1;
2702           tgfirst           = tglast;
2703           tglast            = Standard_False;
2704           ChoixIso          = choixIsoSav;
2705
2706 #if 0 
2707           pasuv[0]=pasSav[0];
2708           pasuv[1]=pasSav[1];
2709           pasuv[2]=pasSav[2];
2710           pasuv[3]=pasSav[3];
2711 #else 
2712           Standard_Real u1,v1,u2,v2;
2713           Standard_Real U1,V1,U2,V2;
2714           Standard_Integer nn=line->NbPoints();
2715           if(nn>2) { 
2716             line->Value(nn).Parameters(u1,v1,u2,v2);
2717             line->Value(nn-1).Parameters(U1,V1,U2,V2);
2718             pasuv[0]=Abs(u1-U1);
2719             pasuv[1]=Abs(v1-V1);
2720             pasuv[2]=Abs(u2-U2);
2721             pasuv[3]=Abs(v2-V2);
2722           }
2723 #endif
2724         }
2725         else Arrive = Standard_True;
2726     }
2727     else {
2728       pasuv[0]*=0.5;
2729       pasuv[1]*=0.5;
2730       pasuv[2]*=0.5;
2731       pasuv[3]*=0.5; 
2732     }
2733   }
2734 }
2735
2736 namespace {
2737   //OCC431(apo): modified ->
2738   static const Standard_Real CosRef2D =  Cos(M_PI/9.0),  AngRef2D = M_PI/2.0; 
2739
2740   static const Standard_Real d = 7.0;
2741 }
2742
2743 IntWalk_StatusDeflection  IntWalk_PWalking::TestDeflection(const IntImp_ConstIsoparametric choixIso)
2744
2745 // test if vector is observed by calculating an increase of vector 
2746 //     or the previous point and its tangent, the new calculated point and its  
2747 //     tangent; it is possible to find a cube passing by the 2 points and having as a 
2748 //     derivative the tangents of the intersection
2749 //     calculate the point with parameter 0.5 on cube=p1 
2750 //     calculate the medium point of 2 points of intersection=p2
2751 //   if arrow/2<=||p1p2||<= arrow consider that the vector is observed
2752 //   otherwise adjust the step depending on the ratio ||p1p2||/vector
2753 //   and the previous step 
2754 // test if in  2 tangent planes of surfaces there is no too great angle2d 
2755 // grand : if yes divide the step
2756 // test if there is no change of side
2757 //  
2758 {
2759   if(line->NbPoints() ==1 ) { 
2760     STATIC_BLOCAGE_SUR_PAS_TROP_GRAND=STATIC_PRECEDENT_INFLEXION=0;
2761   }
2762
2763   IntWalk_StatusDeflection Status = IntWalk_OK;
2764   Standard_Real FlecheCourante , Ratio = 1.0;
2765
2766   // Caro1 and Caro2
2767   const Handle(Adaptor3d_HSurface)& Caro1 = myIntersectionOn2S.Function().AuxillarSurface1();
2768   const Handle(Adaptor3d_HSurface)& Caro2 = myIntersectionOn2S.Function().AuxillarSurface2();
2769
2770   const IntSurf_PntOn2S& CurrentPoint = myIntersectionOn2S.Point(); 
2771   //==================================================================================
2772   //=========               S t o p   o n   p o i n t                 ============
2773   //================================================================================== 
2774   if (myIntersectionOn2S.IsTangent())  { 
2775     return IntWalk_ArretSurPoint;  
2776   }
2777
2778   const gp_Dir& TgCourante = myIntersectionOn2S.Direction();
2779
2780   const Standard_Real aCosBetweenTangent = TgCourante.Dot(previousd);
2781
2782   //==================================================================================
2783   //=========   R i s k   o f    i n f l e x i o n   p o i n t  ============
2784   //==================================================================================  
2785   if (aCosBetweenTangent < 0) {
2786     //------------------------------------------------------------
2787     //-- Risk of inflexion point : Divide the step by 2
2788     //-- Initialize STATIC_PRECEDENT_INFLEXION so that 
2789     //-- at the next call to return Pas_OK if there is no 
2790     //-- more risk of the point of inflexion
2791     //------------------------------------------------------------
2792
2793     pasuv[0]*=0.5;
2794     pasuv[1]*=0.5;
2795     pasuv[2]*=0.5;
2796     pasuv[3]*=0.5;
2797     STATIC_PRECEDENT_INFLEXION+=3; 
2798     if (pasuv[0] < ResoU1 && pasuv[1] <ResoV1 && pasuv[2] <ResoU2 && pasuv[3] < ResoV2)
2799       return IntWalk_ArretSurPointPrecedent;
2800     else 
2801       return IntWalk_PasTropGrand;
2802   }
2803   else {
2804     if(STATIC_PRECEDENT_INFLEXION  > 0) { 
2805       STATIC_PRECEDENT_INFLEXION -- ;
2806       return IntWalk_OK;
2807     }
2808   }
2809
2810   //==================================================================================
2811   //=========  D e t e c t    c o n f u s e d    P o in t s       ===========
2812   //==================================================================================
2813
2814   const Standard_Real aSqDist = previousPoint.Value().
2815     SquareDistance(CurrentPoint.Value());
2816
2817
2818   if (aSqDist < tolconf*tolconf) { 
2819     pasInit[0] = Max(pasInit[0], 5.0*ResoU1);
2820     pasInit[1] = Max(pasInit[1], 5.0*ResoV1);
2821     pasInit[2] = Max(pasInit[2], 5.0*ResoU2);
2822     pasInit[3] = Max(pasInit[3], 5.0*ResoV2);
2823
2824     for(Standard_Integer i = 0; i < 4; i++)
2825     {
2826       pasuv[i] = Max(pasuv[i], Min(1.5*pasuv[i], pasInit[i]));
2827     }
2828     //Compute local resolution: for OCC26717
2829     if (Abs(pasuv[choixIso] - pasInit[choixIso]) <= Precision::Confusion())
2830     {
2831       Standard_Real CurU, CurV;
2832       if (choixIso == IntImp_UIsoparametricOnCaro1 ||
2833           choixIso == IntImp_VIsoparametricOnCaro1)
2834         previousPoint.ParametersOnS1(CurU, CurV);
2835       else
2836         previousPoint.ParametersOnS2(CurU, CurV);
2837       gp_Pnt CurPnt = (choixIso == IntImp_UIsoparametricOnCaro1 ||
2838                        choixIso == IntImp_VIsoparametricOnCaro1)?
2839         Adaptor3d_HSurfaceTool::Value(Caro1, CurU, CurV) :
2840         Adaptor3d_HSurfaceTool::Value(Caro2, CurU, CurV);
2841       gp_Pnt OffsetPnt;
2842       switch(choixIso)
2843       {
2844       case IntImp_UIsoparametricOnCaro1:
2845         OffsetPnt =
2846           Adaptor3d_HSurfaceTool::Value(Caro1,
2847                                         CurU + sensCheminement*pasuv[0],
2848                                         CurV);
2849         break;
2850       case IntImp_VIsoparametricOnCaro1:
2851         OffsetPnt =
2852           Adaptor3d_HSurfaceTool::Value(Caro1,
2853                                         CurU,
2854                                         CurV + sensCheminement*pasuv[1]);
2855         break;
2856       case IntImp_UIsoparametricOnCaro2:
2857         OffsetPnt =
2858           Adaptor3d_HSurfaceTool::Value(Caro2,
2859                                         CurU + sensCheminement*pasuv[2],
2860                                         CurV);
2861         break;
2862       case IntImp_VIsoparametricOnCaro2:
2863         OffsetPnt =
2864           Adaptor3d_HSurfaceTool::Value(Caro2,
2865                                         CurU,
2866                                         CurV + sensCheminement*pasuv[3]);
2867         break;
2868       default:break;
2869       }
2870       Standard_Real RefDist = CurPnt.Distance(OffsetPnt);
2871       Standard_Real LocalResol = 0.;
2872       if (RefDist > gp::Resolution())
2873         LocalResol = pasuv[choixIso] * tolconf / RefDist;
2874       if (pasuv[choixIso] < 2*LocalResol)
2875         pasuv[choixIso] = pasInit[choixIso] = 2*LocalResol;
2876     }
2877     ////////////////////////////////////////
2878     Status = IntWalk_PointConfondu;
2879   }
2880
2881   //==================================================================================
2882   Standard_Real Up1,Vp1,Uc1,Vc1,Du1,Dv1,AbsDu1,AbsDu2,AbsDv1,AbsDv2;
2883   Standard_Real Up2,Vp2,Uc2,Vc2,Du2,Dv2;
2884
2885   previousPoint.Parameters(Up1,Vp1,Up2,Vp2);
2886   CurrentPoint.Parameters(Uc1,Vc1,Uc2,Vc2);               
2887
2888   Du1 = Uc1 - Up1;   Dv1 = Vc1 - Vp1;
2889   Du2 = Uc2 - Up2;   Dv2 = Vc2 - Vp2;
2890
2891   AbsDu1 = Abs(Du1);
2892   AbsDu2 = Abs(Du2);
2893   AbsDv1 = Abs(Dv1);
2894   AbsDv2 = Abs(Dv2);
2895   //=================================================================================
2896   //====   S t e p   o f   p  r o g r e s s i o n (between previous and Current)   =======
2897   //=================================================================================
2898   if (   AbsDu1 < ResoU1 && AbsDv1 < ResoV1 
2899     && AbsDu2 < ResoU2 && AbsDv2 < ResoV2) {
2900       pasuv[0] = ResoU1; pasuv[1] = ResoV1; pasuv[2] = ResoU2; pasuv[3] = ResoV2;
2901       return(IntWalk_ArretSurPointPrecedent);
2902   }
2903   //==================================================================================
2904
2905   Standard_Real tolArea = 100.0;
2906   if (ResoU1 < Precision::PConfusion() ||
2907     ResoV1 < Precision::PConfusion() ||
2908     ResoU2 < Precision::PConfusion() ||
2909     ResoV2 < Precision::PConfusion() )
2910     tolArea =  tolArea*2.0;
2911
2912   Standard_Real Cosi1, CosRef1, Ang1, AngRef1, ResoUV1, Duv1, d1, tolCoeff1;   
2913   Standard_Real Cosi2, CosRef2, Ang2, AngRef2, ResoUV2, Duv2, d2, tolCoeff2;   
2914   Cosi1 = Du1*previousd1.X() + Dv1*previousd1.Y();
2915   Cosi2 = Du2*previousd2.X() + Dv2*previousd2.Y();
2916   Duv1 = Du1*Du1 + Dv1*Dv1;
2917   Duv2 = Du2*Du2 + Dv2*Dv2;
2918   ResoUV1 = ResoU1*ResoU1 + ResoV1*ResoV1;
2919   ResoUV2 = ResoU2*ResoU2 + ResoV2*ResoV2;
2920   //
2921   //modified by NIZNHY-PKV Wed Nov 13 12:25:44 2002 f
2922   //
2923   Standard_Real aMinDiv2=Precision::Confusion();
2924   aMinDiv2=aMinDiv2*aMinDiv2;
2925   //
2926   d1=d;
2927   if (Duv1>aMinDiv2)  {
2928     d1 = Abs(ResoUV1/Duv1);
2929     d1 = Min(Sqrt(d1)*tolArea, d);  
2930   } 
2931   //d1 = Abs(ResoUV1/Duv1); 
2932   //d1 = Min(Sqrt(d1)*tolArea,d);  
2933   //modified by NIZNHY-PKV Wed Nov 13 12:34:30 2002 t
2934   tolCoeff1 = Exp(d1);
2935   //
2936   //modified by NIZNHY-PKV Wed Nov 13 12:34:43 2002 f
2937   d2=d;
2938   if (Duv2>aMinDiv2) {
2939     d2 = Abs(ResoUV2/Duv2); 
2940     d2 = Min(Sqrt(d2)*tolArea,d); 
2941   }
2942   //d2 = Abs(ResoUV2/Duv2); 
2943   //d2 = Min(Sqrt(d2)*tolArea,d);  
2944   //modified by NIZNHY-PKV Wed Nov 13 12:34:53 2002 t
2945   tolCoeff2 = Exp(d2);
2946   CosRef1 = CosRef2D/tolCoeff1;
2947   CosRef2 = CosRef2D/tolCoeff2;
2948   //
2949   //==================================================================================
2950   //== The points are not confused :                                           ==
2951   //== 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 ==
2952   //==                           N o t    T o o    G r e a t (angle in space UV)    ==
2953   //==                           C h a n g e    o f    s i d e                ==
2954   //==================================================================================
2955   if (Status != IntWalk_PointConfondu) { 
2956     if(Cosi1*Cosi1 < CosRef1*Duv1 || Cosi2*Cosi2 < CosRef2*Duv2) {
2957       pasuv[0]*=0.5;  pasuv[1]*=0.5;  pasuv[2]*=0.5;  pasuv[3]*=0.5;
2958       if (pasuv[0]<ResoU1 && pasuv[1]<ResoV1 && pasuv[2]<ResoU2 && pasuv[3]<ResoV2) { 
2959         return(IntWalk_ArretSurPointPrecedent);
2960       }
2961       else {
2962         pasuv[0]*=0.5; pasuv[1]*=0.5; pasuv[2]*=0.5; pasuv[3]*=0.5;
2963         return(IntWalk_PasTropGrand);
2964       }
2965     }
2966     const gp_Dir2d& Tg2dcourante1 = myIntersectionOn2S.DirectionOnS1();
2967     const gp_Dir2d& Tg2dcourante2 = myIntersectionOn2S.DirectionOnS2();
2968     Cosi1 = Du1*Tg2dcourante1.X() + Dv1*Tg2dcourante1.Y();
2969     Cosi2 = Du2*Tg2dcourante2.X() + Dv2*Tg2dcourante2.Y();
2970     Ang1 = Abs(previousd1.Angle(Tg2dcourante1));  
2971     Ang2 = Abs(previousd2.Angle(Tg2dcourante2));  
2972     AngRef1 = AngRef2D*tolCoeff1;
2973     AngRef2 = AngRef2D*tolCoeff2;
2974     //-------------------------------------------------------
2975     //-- Test : Angle too great in space UV       -----
2976     //--        Change of  side                      -----
2977     //-------------------------------------------------------
2978     if(Cosi1*Cosi1 < CosRef1*Duv1 || Cosi2*Cosi2 < CosRef2*Duv2 || Ang1 > AngRef1 || Ang2 > AngRef2) {
2979       pasuv[0]*=0.5;  pasuv[1]*=0.5;  pasuv[2]*=0.5;  pasuv[3]*=0.5;
2980       if (pasuv[0]<ResoU1 && pasuv[1]<ResoV1 && pasuv[2]<ResoU2 && pasuv[3]<ResoV2) 
2981         return(IntWalk_ArretSurPoint);
2982       else 
2983         return(IntWalk_PasTropGrand);
2984     }
2985   }
2986   //<-OCC431(apo)
2987   //==================================================================================
2988   //== D e t e c t i o n   o f    :  Step Too Small 
2989   //==                               STEP TOO Great 
2990   //==================================================================================
2991
2992   //---------------------------------------
2993   //-- Estimate of the vector           --
2994   //---------------------------------------
2995   FlecheCourante =
2996     Sqrt(Abs((previousd.XYZ()-TgCourante.XYZ()).SquareModulus()*aSqDist))/8.;
2997
2998   if ( FlecheCourante<= fleche*0.5) {     //-- Current step too small
2999     if(FlecheCourante>1e-16) { 
3000       Ratio = 0.5*(fleche/FlecheCourante);
3001     }
3002     else { 
3003       Ratio = 10.0;
3004     }
3005     Standard_Real pasSu1 = pasuv[0];
3006     Standard_Real pasSv1 = pasuv[1];
3007     Standard_Real pasSu2 = pasuv[2];
3008     Standard_Real pasSv2 = pasuv[3];
3009
3010     //-- In  case if 
3011     //-- a point at U+DeltaU is required, ....
3012     //-- return a point at U + Epsilon
3013     //-- Epsilon << DeltaU.
3014
3015     if(pasuv[0]< AbsDu1) pasuv[0] = AbsDu1;
3016     if(pasuv[1]< AbsDv1) pasuv[1] = AbsDv1;
3017     if(pasuv[2]< AbsDu2) pasuv[2] = AbsDu2;
3018     if(pasuv[3]< AbsDv2) pasuv[3] = AbsDv2;
3019
3020     if(pasuv[0]<ResoU1) pasuv[0]=ResoU1;
3021     if(pasuv[1]<ResoV1) pasuv[1]=ResoV1;
3022     if(pasuv[2]<ResoU2) pasuv[2]=ResoU2;
3023     if(pasuv[3]<ResoV2) pasuv[3]=ResoV2;
3024     //-- if(Ratio>10.0 ) { Ratio=10.0; } 
3025     Standard_Real R1,R = pasInit[0]/pasuv[0];
3026     R1= pasInit[1]/pasuv[1];     if(R1<R) R=R1;
3027     R1= pasInit[2]/pasuv[2];     if(R1<R) R=R1;
3028     R1= pasInit[3]/pasuv[3];     if(R1<R) R=R1;
3029     if(Ratio > R) Ratio=R;
3030     pasuv[0] = Min(Ratio*pasuv[0],pasInit[0]);
3031     pasuv[1] = Min(Ratio*pasuv[1],pasInit[1]);
3032     pasuv[2] = Min(Ratio*pasuv[2],pasInit[2]);
3033     pasuv[3] = Min(Ratio*pasuv[3],pasInit[3]);
3034     if (pasuv[0] != pasSu1 || pasuv[2] != pasSu2|| 
3035       pasuv[1] != pasSv1 || pasuv[3] != pasSv2) {
3036         if(++STATIC_BLOCAGE_SUR_PAS_TROP_GRAND > 5) {
3037           STATIC_BLOCAGE_SUR_PAS_TROP_GRAND = 0;
3038           return IntWalk_PasTropGrand; 
3039         }
3040     }
3041     if(Status == IntWalk_OK) { 
3042       STATIC_BLOCAGE_SUR_PAS_TROP_GRAND=0;
3043       //-- Try to increase the step
3044     }
3045     return Status;
3046   }
3047   else {                                //-- CurrentVector > vector*0.5 
3048     if (FlecheCourante > fleche) {      //-- Current step too Great
3049       Ratio = fleche/FlecheCourante; 
3050       pasuv[0] = Ratio*pasuv[0];
3051       pasuv[1] = Ratio*pasuv[1];
3052       pasuv[2] = Ratio*pasuv[2];
3053       pasuv[3] = Ratio*pasuv[3];
3054       //if(++STATIC_BLOCAGE_SUR_PAS_TROP_GRAND > 5) {
3055       //        STATIC_BLOCAGE_SUR_PAS_TROP_GRAND = 0;
3056       return IntWalk_PasTropGrand; 
3057       //}
3058     }
3059     else {                             //-- vector/2  <  CurrentVector <= vector   
3060       Ratio = 0.75 * (fleche / FlecheCourante);
3061     }
3062   }
3063
3064   if(Status != IntWalk_PointConfondu)
3065   {
3066     //Here, aCosBetweenTangent >= 0.0 definitely.
3067
3068     /*
3069     Brief algorithm description.
3070     We have two (not-coincindent) intersection point (P1 and P2). In every point,
3071     vector of tangent (to the intersection curve) is known (vectors T1 and T2).
3072     Basing on these data, we create osculating circle.
3073
3074                                     * - arc of osculating circle
3075                                 *      * 
3076                            P1 x----------x P2
3077                              /            \
3078                             /              \
3079                           Vec(T1)         Vec(T2)
3080
3081     Let me draw your attention to the following facts:
3082       1. Vectors T1 and T2 direct FROM (not TO) points P1 and P2. Therefore,
3083         one of previously computed vector should be reversed.
3084
3085     In this case, the absolute (!) value of the deflection between the arc of
3086     the osculating circle and the P1P2 segment can be computed as follows:
3087           e = d*(1-sin(B/2))/(2*cos(B/2)),      (1)
3088     where d is the length of P1P2 segment, B is the angle between vectors T1 and T2.
3089     At that,
3090               pi/2 <= B <= pi,
3091               cos(B/2) >= 0,
3092               sin(B/2) > 0,
3093               sin(B) > 0,
3094               cos(B) < 0.
3095
3096     Later, the normal state of algorithm work is (as we apply) 
3097               tolconf/2 <= e <= tolconf.
3098     In this case, we keep previous step.
3099
3100     If e < tolconf/2 then the local curvature of the intersection curve is small.
3101     As result, the step should be increased.
3102
3103     If e > tolconf then the step is too big. Therefore, we should decrease one.
3104
3105     Condition (1) is equivalent to
3106               sin(B/2) = 1 - 2/(1+(d/(2*e))^2) = Fs(e),
3107               cos(B) = 1 - 2*Fs(e)^2 = Fd(e),
3108     where Fs(e)and Fd(e) are some function with parameter "deflection".
3109
3110     Let mean that Fs(e) is decreasing function. Fd(e) is increasing function,
3111     in the range, where Fs(e) > 0.0 (i.e. when e < d/2).
3112
3113     Now, let substitute required deflection (tolconf or tolconf/2) to z. Then
3114     it is necessary to check if e < z or if e > z.
3115
3116     In this case, it is enough to comapare Fs(e) and Fs(z).
3117     At that Fs(e) > 0 because sin(B/2) > 0 always.
3118
3119     Therefore, if Fs(z) < 0.0 then Fs(e) > Fs(z) ==> e < z definitely.
3120     If Fs(z) > 0.0 then we can compare Fs(z)^2 and Fs(e)^2 or, in substance,
3121     values Fd(e) and Fd(z). If Fd(e) > Fd(z) then e > z and vice versa.
3122     */
3123
3124     //Fd(e) is already known (Fd(e) == -aCosBetweenTangent)
3125
3126     const Standard_Real anInvSqAbsArcDeflMax = 0.25*aSqDist/(tolconf*tolconf);
3127     const Standard_Real aSinB2Max = 1.0 - 2.0/(1.0 + anInvSqAbsArcDeflMax);
3128
3129     if(aSinB2Max >= 0.0 && (aCosBetweenTangent <= 2.0 * aSinB2Max * aSinB2Max - 1.0))
3130     {//Real deflection is greater or equal than tolconf
3131       Status = IntWalk_PasTropGrand;
3132     }
3133     else
3134     {//Real deflection is less than tolconf
3135       const Standard_Real anInvSqAbsArcDeflMin = 4.0*anInvSqAbsArcDeflMax;
3136       const Standard_Real aSinB2Min = 1.0 - 2.0/(1.0 + anInvSqAbsArcDeflMin);
3137
3138       if((aSinB2Min < 0.0) || (aCosBetweenTangent >= 2.0 * aSinB2Min * aSinB2Min - 1.0))
3139       {//Real deflection is less than tolconf/2.0
3140         Status = IntWalk_StepTooSmall;
3141       }
3142     }
3143
3144     if(Status == IntWalk_PasTropGrand)
3145     {
3146       pasuv[0]*=0.5; pasuv[1]*=0.5; pasuv[2]*=0.5; pasuv[3]*=0.5;
3147       return Status;
3148     }
3149
3150     if(Status == IntWalk_StepTooSmall)
3151     {
3152       pasuv[0] = Max(pasuv[0], AbsDu1);
3153       pasuv[1] = Max(pasuv[1], AbsDv1);
3154       pasuv[2] = Max(pasuv[2], AbsDu2);
3155       pasuv[3] = Max(pasuv[3], AbsDv2);
3156
3157       pasInit[0] = Max(pasInit[0], AbsDu1);
3158       pasInit[1] = Max(pasInit[1], AbsDv1);
3159       pasInit[2] = Max(pasInit[2], AbsDu2);
3160       pasInit[3] = Max(pasInit[3], AbsDv2);
3161
3162       return Status;
3163     }
3164   }
3165
3166   pasuv[0] = Max(myStepMin[0],Min(Min(Ratio*AbsDu1,pasuv[0]),pasInit[0]));
3167   pasuv[1] = Max(myStepMin[1],Min(Min(Ratio*AbsDv1,pasuv[1]),pasInit[1]));
3168   pasuv[2] = Max(myStepMin[2],Min(Min(Ratio*AbsDu2,pasuv[2]),pasInit[2]));
3169   pasuv[3] = Max(myStepMin[3],Min(Min(Ratio*AbsDv2,pasuv[3]),pasInit[3]));
3170
3171   if(Status == IntWalk_OK) STATIC_BLOCAGE_SUR_PAS_TROP_GRAND=0;
3172   return Status;
3173 }
3174
3175 Standard_Boolean IntWalk_PWalking::
3176 TestArret(const Standard_Boolean DejaReparti,
3177           TColStd_Array1OfReal& Param,
3178           IntImp_ConstIsoparametric&  ChoixIso)
3179
3180           //
3181           // test if the point of intersection set by these parameters remains in the 
3182           // natural domain of each square.
3183           // if the point outpasses reframe to find the best iso (border)
3184           // that intersects easiest the other square
3185           // otherwise test if closed line is present  
3186           // 
3187 {
3188   Standard_Real Uvd[4],Uvf[4],Epsuv[4],Duv[4],Uvp[4],dv,dv2,ParC[4];
3189   Standard_Real DPc,DPb;
3190   Standard_Integer i = 0, k = 0;
3191   Epsuv[0] = ResoU1;
3192   Epsuv[1] = ResoV1;
3193   Epsuv[2] = ResoU2;
3194   Epsuv[3] = ResoV2;
3195   previousPoint.Parameters(Uvp[0],Uvp[1],Uvp[2],Uvp[3]);
3196
3197   Standard_Real SolParam[4];
3198   myIntersectionOn2S.Point().Parameters(SolParam[0],SolParam[1],SolParam[2],SolParam[3]);
3199
3200   Standard_Boolean Trouve = Standard_False;
3201
3202   Uvd[0]=Um1;   Uvf[0]=UM1;   Uvd[1]=Vm1;   Uvf[1]=VM1;
3203   Uvd[2]=Um2;   Uvf[2]=UM2;   Uvd[3]=Vm2;   Uvf[3]=VM2;
3204
3205   Standard_Integer im1;
3206   for ( i = 1,im1 = 0;i<=4;i++,im1++) {
3207     switch(i) { 
3208     case 1: k=2; break;
3209     case 2: k=1; break;
3210     case 3: k=4; break;
3211     case 4: k=3; break;
3212     }
3213     if (Param(i) < (Uvd[im1]-Epsuv[im1]) ||
3214       SolParam[im1] < (Uvd[im1]-Epsuv[im1]))     //--     Current -----  Bound Inf -----  Previous
3215     {
3216       Trouve    = Standard_True;                   //-- 
3217       DPc       = Uvp[im1]-Param(i);               //--     Previous  - Current
3218       DPb       = Uvp[im1]-Uvd[im1];               //--     Previous  - Bound Inf
3219       ParC[im1] = Uvd[im1];                        //--     ParamCorrige
3220       dv        = Param(k)-Uvp[k-1];               //--     Current   - Previous (other Direction)
3221       dv2       = dv*dv;         
3222       if(dv2>RealEpsilon()) {                       //--    Progress at the other Direction ?
3223         Duv[im1]  = DPc*DPb + dv2;
3224         Duv[im1]  = Duv[im1]*Duv[im1]/(DPc*DPc+dv2)/(DPb*DPb+dv2);
3225       }
3226       else {
3227         Duv[im1]=-1.0;                              //--    If no progress, do not change  
3228       }                                             //--    the choice of iso 
3229     }   
3230     else if (Param(i) > (Uvf[im1] + Epsuv[im1]) ||
3231       SolParam[im1] > (Uvf[im1] + Epsuv[im1]))//--    Previous -----  Bound Sup -----  Current
3232     {
3233       Trouve    = Standard_True;                    //-- 
3234       DPc       = Param(i)-Uvp[im1];                //--     Current   - Previous
3235       DPb       = Uvf[im1]-Uvp[im1];                //--     Bound Sup - Previous 
3236       ParC[im1] = Uvf[im1];                         //--     Param Corrige
3237       dv        = Param(k)-Uvp[k-1];                //--     Current   - Previous (other Direction)
3238       dv2       = dv*dv;
3239       if(dv2>RealEpsilon()) {                       //--     Progress in other Direction ?
3240         Duv[im1]  =  DPc*DPb + dv2;
3241         Duv[im1]  = Duv[im1]*Duv[im1]/(DPc*DPc+dv2)/(DPb*DPb+dv2);
3242       }
3243       else {
3244         Duv[im1]=-1.0;                              //--    If no progress, do not change 
3245       }                                             //--    the choice of iso 
3246     }
3247     else { 
3248       Duv[im1]= -1.;
3249       ParC[im1]=Param(i);
3250     }
3251   }
3252
3253   if (Trouve) {
3254     //--------------------------------------------------
3255     //-- One of Parameters u1,v1,u2,v2 is outside of  --
3256     //-- the natural limits.                          -- 
3257     //-- Find the best direction of                   -- 
3258     //-- progress and reframe the parameters.        --
3259     //--------------------------------------------------
3260     Standard_Real ddv = -1.0;
3261     k=-1;
3262     for (i=0;i<=3;i++) {
3263       Param(i+1) = ParC[i];
3264       if(Duv[i]>ddv) { 
3265         ddv = Duv[i];
3266         k=i;
3267       }
3268     }
3269     if(k!=-1) { 
3270       ChoixIso   = ChoixRef[k];
3271     }
3272     else { 
3273       if((ParC[0]<=Uvd[0]+Epsuv[0]) || (ParC[0]>=Uvf[0]-Epsuv[0])) {
3274         ChoixIso = IntImp_UIsoparametricOnCaro1;
3275       }
3276       else if((ParC[1]<=Uvd[1]+Epsuv[1]) || (ParC[1]>=Uvf[1]-Epsuv[1])) {
3277         ChoixIso = IntImp_VIsoparametricOnCaro1;
3278       }
3279       else if((ParC[2]<=Uvd[2]+Epsuv[2]) || (ParC[2]>=Uvf[2]-Epsuv[2])) {
3280         ChoixIso = IntImp_UIsoparametricOnCaro2;
3281       }
3282       else if((ParC[3]<=Uvd[3]+Epsuv[3]) || (ParC[3]>=Uvf[3]-Epsuv[3])) {
3283         ChoixIso = IntImp_VIsoparametricOnCaro2;
3284       }
3285     }
3286     close = Standard_False;
3287     return Standard_True;
3288   }
3289   else 
3290   {  
3291     if (!DejaReparti) { // find if line closed
3292
3293       Standard_Real u,v;
3294       const IntSurf_PntOn2S& POn2S1=line->Value(1);
3295       //On S1
3296       POn2S1.ParametersOnS1(u,v);
3297       gp_Pnt2d P1uvS1(u,v);
3298       previousPoint.ParametersOnS1(u,v);
3299       gp_Pnt2d PrevuvS1(u,v);
3300       myIntersectionOn2S.Point().ParametersOnS1(u,v);
3301       gp_Pnt2d myIntersuvS1(u,v);
3302       Standard_Boolean close2dS1 = (P1uvS1.XY()-PrevuvS1.XY())*
3303         (P1uvS1.XY()-myIntersuvS1.XY()) < 0.0;
3304       //On S2
3305       POn2S1.ParametersOnS2(u,v);
3306       gp_Pnt2d P1uvS2(u,v);
3307       previousPoint.ParametersOnS2(u,v);
3308       gp_Pnt2d PrevuvS2(u,v);
3309       myIntersectionOn2S.Point().ParametersOnS2(u,v);
3310       gp_Pnt2d myIntersuvS2(u,v);
3311       Standard_Boolean close2dS2 = (P1uvS2.XY()-PrevuvS2.XY())*
3312         (P1uvS2.XY()-myIntersuvS2.XY()) < 0.0;
3313
3314       close = close2dS1 && close2dS2;
3315       return close;
3316     }
3317     else return Standard_False;
3318   }
3319 }
3320