src/IntWalk/IntWalk_PWalking.cxx

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