src/IntWalk/IntWalk_PWalking.cxx

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