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