OCCT Git - occt.git/blob

   1 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   2 //
   3 // This file is part of Open CASCADE Technology software library.
   4 //
   5 // This library is free software; you can redistribute it and/or modify it under
   6 // the terms of the GNU Lesser General Public License version 2.1 as published
   7 // by the Free Software Foundation, with special exception defined in the file
   8 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
   9 // distribution for complete text of the license and disclaimer of any warranty.
  10 //
  11 // Alternatively, this file may be used under the terms of Open CASCADE
  12 // commercial license or contractual agreement.
  13
  14 // 06.01.99 pdn private method SurfaceNewton PRO17015: fix against hang in Extrema
  15 // 11.01.99 pdn PRO10109 4517: protect against wrong result
  16 //%12 pdn 11.02.99 PRO9234 project degenerated
  17 // 03.03.99 rln S4135: new algorithms for IsClosed (accepts precision), Degenerated (stores
  18 // precision)
  19 //: p7 abv 10.03.99 PRO18206: adding new method IsDegenerated()
  20 //: p8 abv 11.03.99 PRO7226 #489490: improving ProjectDegenerated() for degenerated edges
  21 //: q1 pdn, abv 15.03.99 PRO7226 #525030: adding maxpreci in NextValueOfUV()
  22 //: q2 abv 16.03.99: PRO7226 #412920: applying Newton algorithm before UVFromIso()
  23 //: q6 abv 19.03.99: ie_soapbox-B.stp #390760: improving projecting point on surface
  24 // #77 rln 15.03.99: S4135: returning singularity which has minimum gap between singular point and
  25 // input 3D point
  26 //: r3 abv 30.03.99: (by smh) S4163: protect against unexpected signals
  27 //: #4 smh 07.04.99: S4163 Zero divide.
  28 // #4  szv           S4163: optimizations
  29 //: r9 abv 09.04.99: id_turbine-C.stp #3865: check degenerated 2d point by recomputing to 3d instead
  30 //: of Resolution s5 abv 22.04.99  Adding debug printouts in catch {} blocks
  31
  32 #include <Adaptor3d_Curve.hxx>
  33 #include <Adaptor3d_IsoCurve.hxx>
  34 #include <BndLib_Add3dCurve.hxx>
  35 #include <ElSLib.hxx>
  36 #include <Geom_BezierSurface.hxx>
  37 #include <Geom_BoundedSurface.hxx>
  38 #include <Geom_ConicalSurface.hxx>
  39 #include <Geom_Curve.hxx>
  40 #include <Geom_OffsetSurface.hxx>
  41 #include <Geom_RectangularTrimmedSurface.hxx>
  42 #include <Geom_SphericalSurface.hxx>
  43 #include <Geom_Surface.hxx>
  44 #include <Geom_SurfaceOfLinearExtrusion.hxx>
  45 #include <Geom_SurfaceOfRevolution.hxx>
  46 #include <Geom_ToroidalSurface.hxx>
  47 #include <GeomAdaptor_Curve.hxx>
  48 #include <GeomAdaptor_Surface.hxx>
  49 #include <gp_Pnt.hxx>
  50 #include <gp_Pnt2d.hxx>
  51 #include <Precision.hxx>
  52 #include <ShapeAnalysis.hxx>
  53 #include <ShapeAnalysis_Curve.hxx>
  54 #include <ShapeAnalysis_Surface.hxx>
  55 #include <Standard_ErrorHandler.hxx>
  56 #include <Standard_Failure.hxx>
  57 #include <Standard_Type.hxx>
  58
  59 IMPLEMENT_STANDARD_RTTIEXT(ShapeAnalysis_Surface, Standard_Transient)
  60
  61 namespace
  62 {
  63 inline void RestrictBounds(double& theFirst, double& theLast)
  64 {
  65   Standard_Boolean isFInf = Precision::IsNegativeInfinite(theFirst);
  66   Standard_Boolean isLInf = Precision::IsPositiveInfinite(theLast);
  67   if (isFInf || isLInf)
  68   {
  69     if (isFInf && isLInf)
  70     {
  71       theFirst = -1000;
  72       theLast  = 1000;
  73     }
  74     else if (isFInf)
  75     {
  76       theFirst = theLast - 2000;
  77     }
  78     else
  79     {
  80       theLast = theFirst + 2000;
  81     }
  82   }
  83 }
  84
  85 inline void RestrictBounds(double& theUf, double& theUl, double& theVf, double& theVl)
  86 {
  87   RestrictBounds(theUf, theUl);
  88   RestrictBounds(theVf, theVl);
  89 }
  90 } // namespace
  91
  92 // S4135
  93 // S4135
  94 //=================================================================================================
  95
  96 ShapeAnalysis_Surface::ShapeAnalysis_Surface(const Handle(Geom_Surface)& S)
  97     : mySurf(S),
  98       myExtOK(Standard_False), //: 30
  99       myNbDeg(-1),
 100       myIsos(Standard_False),
 101       myIsoBoxes(Standard_False),
 102       myGap(0.),
 103       myUDelt(0.01),
 104       myVDelt(0.01),
 105       myUCloseVal(-1),
 106       myVCloseVal(-1)
 107 {
 108   mySurf->Bounds(myUF, myUL, myVF, myVL);
 109   myAdSur = new GeomAdaptor_Surface(mySurf);
 110 }
 111
 112 //=================================================================================================
 113
 114 void ShapeAnalysis_Surface::Init(const Handle(Geom_Surface)& S)
 115 {
 116   if (mySurf == S)
 117     return;
 118   myExtOK     = Standard_False; //: 30
 119   mySurf      = S;
 120   myNbDeg     = -1;
 121   myUCloseVal = myVCloseVal = -1;
 122   myGap                     = 0.;
 123   mySurf->Bounds(myUF, myUL, myVF, myVL);
 124   myAdSur    = new GeomAdaptor_Surface(mySurf);
 125   myIsos     = Standard_False;
 126   myIsoBoxes = Standard_False;
 127   myIsoUF.Nullify();
 128   myIsoUL.Nullify();
 129   myIsoVF.Nullify();
 130   myIsoVL.Nullify();
 131 }
 132
 133 //=================================================================================================
 134
 135 void ShapeAnalysis_Surface::Init(const Handle(ShapeAnalysis_Surface)& other)
 136 {
 137   Init(other->Surface());
 138   myAdSur = other->TrueAdaptor3d();
 139   myNbDeg = other->myNbDeg; // rln S4135 direct transmission (to avoid computation in <other>)
 140   for (Standard_Integer i = 0; i < myNbDeg; i++)
 141   {
 142     other->Singularity(i + 1,
 143                        myPreci[i],
 144                        myP3d[i],
 145                        myFirstP2d[i],
 146                        myLastP2d[i],
 147                        myFirstPar[i],
 148                        myLastPar[i],
 149                        myUIsoDeg[i]);
 150   }
 151 }
 152
 153 //=================================================================================================
 154
 155 const Handle(GeomAdaptor_Surface)& ShapeAnalysis_Surface::Adaptor3d()
 156 {
 157   return myAdSur;
 158 }
 159
 160 //=================================================================================================
 161
 162 void ShapeAnalysis_Surface::ComputeSingularities()
 163 {
 164   // rln S4135
 165   if (myNbDeg >= 0)
 166     return;
 167   //: 51 abv 22 Dec 97: allow forcing:  if (myNbDeg >= 0) return;
 168   // CKY 27-FEV-98 : en appel direct on peut forcer. En appel interne on optimise
 169   if (mySurf.IsNull())
 170     return;
 171
 172   Standard_Real su1, sv1, su2, sv2;
 173   //  mySurf->Bounds(su1, su2, sv1, sv2);
 174   Bounds(su1, su2, sv1, sv2); // modified by rln on 12/11/97 mySurf-> is deleted
 175
 176   myNbDeg = 0; //: r3
 177
 178   if (mySurf->IsKind(STANDARD_TYPE(Geom_ConicalSurface)))
 179   {
 180     Handle(Geom_ConicalSurface) conicS = Handle(Geom_ConicalSurface)::DownCast(mySurf);
 181     Standard_Real               vApex  = -conicS->RefRadius() / Sin(conicS->SemiAngle());
 182     myPreci[0]                         = 0;
 183     myP3d[0]                           = conicS->Apex();
 184     myFirstP2d[0].SetCoord(su1, vApex);
 185     myLastP2d[0].SetCoord(su2, vApex);
 186     myFirstPar[0] = su1;
 187     myLastPar[0]  = su2;
 188     myUIsoDeg[0]  = Standard_False;
 189     myNbDeg       = 1;
 190   }
 191   else if (mySurf->IsKind(STANDARD_TYPE(Geom_ToroidalSurface)))
 192   {
 193     Handle(Geom_ToroidalSurface) toroidS = Handle(Geom_ToroidalSurface)::DownCast(mySurf);
 194     Standard_Real                minorR  = toroidS->MinorRadius();
 195     Standard_Real                majorR  = toroidS->MajorRadius();
 196     // szv#4:S4163:12Mar99 warning - possible div by zero
 197     Standard_Real Ang = ACos(Min(1., majorR / minorR));
 198     myPreci[0] = myPreci[1] = Max(0., majorR - minorR);
 199     myP3d[0]                = mySurf->Value(0., M_PI - Ang);
 200     myFirstP2d[0].SetCoord(su1, M_PI - Ang);
 201     myLastP2d[0].SetCoord(su2, M_PI - Ang);
 202     myP3d[1] = mySurf->Value(0., M_PI + Ang);
 203     myFirstP2d[1].SetCoord(su2, M_PI + Ang);
 204     myLastP2d[1].SetCoord(su1, M_PI + Ang);
 205     myFirstPar[0] = myFirstPar[1] = su1;
 206     myLastPar[0] = myLastPar[1] = su2;
 207     myUIsoDeg[0] = myUIsoDeg[1] = Standard_False;
 208     myNbDeg                     = (majorR > minorR ? 1 : 2);
 209   }
 210   else if (mySurf->IsKind(STANDARD_TYPE(Geom_SphericalSurface)))
 211   {
 212     myPreci[0] = myPreci[1] = 0;
 213     myP3d[0]                = mySurf->Value(su1, sv2); // Northern pole is first
 214     myP3d[1]                = mySurf->Value(su1, sv1);
 215     myFirstP2d[0].SetCoord(su2, sv2);
 216     myLastP2d[0].SetCoord(su1, sv2);
 217     myFirstP2d[1].SetCoord(su1, sv1);
 218     myLastP2d[1].SetCoord(su2, sv1);
 219     myFirstPar[0] = myFirstPar[1] = su1;
 220     myLastPar[0] = myLastPar[1] = su2;
 221     myUIsoDeg[0] = myUIsoDeg[1] = Standard_False;
 222     myNbDeg                     = 2;
 223   }
 224   else if ((mySurf->IsKind(STANDARD_TYPE(Geom_BoundedSurface)))
 225            || (mySurf->IsKind(STANDARD_TYPE(Geom_SurfaceOfRevolution))) || //: b2 abv 18 Feb 98
 226            (mySurf->IsKind(STANDARD_TYPE(Geom_OffsetSurface))))
 227   { // rln S4135
 228
 229     // rln S4135 //:r3
 230     myP3d[0] = myAdSur->Value(su1, 0.5 * (sv1 + sv2));
 231     myFirstP2d[0].SetCoord(su1, sv2);
 232     myLastP2d[0].SetCoord(su1, sv1);
 233
 234     myP3d[1] = myAdSur->Value(su2, 0.5 * (sv1 + sv2));
 235     myFirstP2d[1].SetCoord(su2, sv1);
 236     myLastP2d[1].SetCoord(su2, sv2);
 237
 238     myP3d[2] = myAdSur->Value(0.5 * (su1 + su2), sv1);
 239     myFirstP2d[2].SetCoord(su1, sv1);
 240     myLastP2d[2].SetCoord(su2, sv1);
 241
 242     myP3d[3] = myAdSur->Value(0.5 * (su1 + su2), sv2);
 243     myFirstP2d[3].SetCoord(su2, sv2);
 244     myLastP2d[3].SetCoord(su1, sv2);
 245
 246     myFirstPar[0] = myFirstPar[1] = sv1;
 247     myLastPar[0] = myLastPar[1] = sv2;
 248     myUIsoDeg[0] = myUIsoDeg[1] = Standard_True;
 249
 250     myFirstPar[2] = myFirstPar[3] = su1;
 251     myLastPar[2] = myLastPar[3] = su2;
 252     myUIsoDeg[2] = myUIsoDeg[3] = Standard_False;
 253
 254     gp_Pnt Corner1 = myAdSur->Value(su1, sv1);
 255     gp_Pnt Corner2 = myAdSur->Value(su1, sv2);
 256     gp_Pnt Corner3 = myAdSur->Value(su2, sv1);
 257     gp_Pnt Corner4 = myAdSur->Value(su2, sv2);
 258
 259     myPreci[0] =
 260       Max(Corner1.Distance(Corner2), Max(myP3d[0].Distance(Corner1), myP3d[0].Distance(Corner2)));
 261     myPreci[1] =
 262       Max(Corner3.Distance(Corner4), Max(myP3d[1].Distance(Corner3), myP3d[1].Distance(Corner4)));
 263     myPreci[2] =
 264       Max(Corner1.Distance(Corner3), Max(myP3d[2].Distance(Corner1), myP3d[2].Distance(Corner3)));
 265     myPreci[3] =
 266       Max(Corner2.Distance(Corner4), Max(myP3d[3].Distance(Corner2), myP3d[3].Distance(Corner4)));
 267
 268     myNbDeg = 4;
 269   }
 270   SortSingularities();
 271 }
 272
 273 //=================================================================================================
 274
 275 Standard_Boolean ShapeAnalysis_Surface::HasSingularities(const Standard_Real preci)
 276 {
 277   return NbSingularities(preci) > 0;
 278 }
 279
 280 //=================================================================================================
 281
 282 Standard_Integer ShapeAnalysis_Surface::NbSingularities(const Standard_Real preci)
 283 {
 284   if (myNbDeg < 0)
 285     ComputeSingularities();
 286   Standard_Integer Nb = 0;
 287   for (Standard_Integer i = 1; i <= myNbDeg; i++)
 288     if (myPreci[i - 1] <= preci)
 289       Nb++;
 290   return Nb;
 291 }
 292
 293 //=================================================================================================
 294
 295 Standard_Boolean ShapeAnalysis_Surface::Singularity(const Standard_Integer num,
 296                                                     Standard_Real&         preci,
 297                                                     gp_Pnt&                P3d,
 298                                                     gp_Pnt2d&              firstP2d,
 299                                                     gp_Pnt2d&              lastP2d,
 300                                                     Standard_Real&         firstpar,
 301                                                     Standard_Real&         lastpar,
 302                                                     Standard_Boolean&      uisodeg)
 303 {
 304   //  ATTENTION, les champs sont des tableaux C, n0s partent de 0. num part de 1
 305   if (myNbDeg < 0)
 306     ComputeSingularities();
 307   if (num < 1 || num > myNbDeg)
 308     return Standard_False;
 309   P3d      = myP3d[num - 1];
 310   preci    = myPreci[num - 1];
 311   firstP2d = myFirstP2d[num - 1];
 312   lastP2d  = myLastP2d[num - 1];
 313   firstpar = myFirstPar[num - 1];
 314   lastpar  = myLastPar[num - 1];
 315   uisodeg  = myUIsoDeg[num - 1];
 316   return Standard_True;
 317 }
 318
 319 //=================================================================================================
 320
 321 Standard_Boolean ShapeAnalysis_Surface::IsDegenerated(const gp_Pnt& P3d, const Standard_Real preci)
 322 {
 323   if (myNbDeg < 0)
 324     ComputeSingularities();
 325   for (Standard_Integer i = 0; i < myNbDeg && myPreci[i] <= preci; i++)
 326   {
 327     myGap = myP3d[i].Distance(P3d);
 328     // rln S4135
 329     if (myGap <= preci)
 330       return Standard_True;
 331   }
 332   return Standard_False;
 333 }
 334
 335 //=================================================================================================
 336
 337 Standard_Boolean ShapeAnalysis_Surface::DegeneratedValues(const gp_Pnt&       P3d,
 338                                                           const Standard_Real preci,
 339                                                           gp_Pnt2d&           firstP2d,
 340                                                           gp_Pnt2d&           lastP2d,
 341                                                           Standard_Real&      firstPar,
 342                                                           Standard_Real&      lastPar,
 343                                                           const Standard_Boolean /*forward*/)
 344 {
 345   if (myNbDeg < 0)
 346     ComputeSingularities();
 347   // #77 rln S4135: returning singularity which has minimum gap between singular point and input 3D
 348   // point
 349   Standard_Integer indMin = -1;
 350   Standard_Real    gapMin = RealLast();
 351   for (Standard_Integer i = 0; i < myNbDeg && myPreci[i] <= preci; i++)
 352   {
 353     myGap = myP3d[i].Distance(P3d);
 354     // rln S4135
 355     if (myGap <= preci)
 356       if (gapMin > myGap)
 357       {
 358         gapMin = myGap;
 359         indMin = i;
 360       }
 361   }
 362   if (indMin >= 0)
 363   {
 364     myGap    = gapMin;
 365     firstP2d = myFirstP2d[indMin];
 366     lastP2d  = myLastP2d[indMin];
 367     firstPar = myFirstPar[indMin];
 368     lastPar  = myLastPar[indMin];
 369     return Standard_True;
 370   }
 371   return Standard_False;
 372 }
 373
 374 //=================================================================================================
 375
 376 Standard_Boolean ShapeAnalysis_Surface::ProjectDegenerated(const gp_Pnt&       P3d,
 377                                                            const Standard_Real preci,
 378                                                            const gp_Pnt2d&     neighbour,
 379                                                            gp_Pnt2d&           result)
 380 {
 381   if (myNbDeg < 0)
 382     ComputeSingularities();
 383   // added by rln on 03/12/97
 384   //: c1 abv 23 Feb 98: preci (3d) -> Resolution (2d)
 385   // #77 rln S4135
 386   Standard_Integer indMin = -1;
 387   Standard_Real    gapMin = RealLast();
 388   for (Standard_Integer i = 0; i < myNbDeg && myPreci[i] <= preci; i++)
 389   {
 390     Standard_Real gap2 = myP3d[i].SquareDistance(P3d);
 391     if (gap2 > preci * preci)
 392       gap2 = Min(gap2, myP3d[i].SquareDistance(Value(result)));
 393     // rln S4135
 394     if (gap2 <= preci * preci && gapMin > gap2)
 395     {
 396       gapMin = gap2;
 397       indMin = i;
 398     }
 399   }
 400   if (indMin < 0)
 401     return Standard_False;
 402   myGap = Sqrt(gapMin);
 403   if (!myUIsoDeg[indMin])
 404     result.SetX(neighbour.X());
 405   else
 406     result.SetY(neighbour.Y());
 407   return Standard_True;
 408 }
 409
 410 // pdn %12 11.02.99 PRO9234 entity 15402
 411 //=================================================================================================
 412
 413 Standard_Boolean ShapeAnalysis_Surface::ProjectDegenerated(const Standard_Integer      nbrPnt,
 414                                                            const TColgp_SequenceOfPnt& points,
 415                                                            TColgp_SequenceOfPnt2d&     pnt2d,
 416                                                            const Standard_Real         preci,
 417                                                            const Standard_Boolean      direct)
 418 {
 419   if (myNbDeg < 0)
 420     ComputeSingularities();
 421
 422   Standard_Integer step = (direct ? 1 : -1);
 423   // #77 rln S4135
 424   Standard_Integer indMin = -1;
 425   Standard_Real    gapMin = RealLast(), prec2 = preci * preci;
 426   Standard_Integer j = (direct ? 1 : nbrPnt);
 427   for (Standard_Integer i = 0; i < myNbDeg && myPreci[i] <= preci; i++)
 428   {
 429     Standard_Real gap2 = myP3d[i].SquareDistance(points(j));
 430     if (gap2 > prec2)
 431       gap2 = Min(gap2, myP3d[i].SquareDistance(Value(pnt2d(j))));
 432     if (gap2 <= prec2 && gapMin > gap2)
 433     {
 434       gapMin = gap2;
 435       indMin = i;
 436     }
 437   }
 438   if (indMin < 0)
 439     return Standard_False;
 440
 441   myGap = Sqrt(gapMin);
 442   gp_Pnt2d pk;
 443
 444   Standard_Integer k; // svv Jan11 2000 : porting on DEC
 445   for (k = j + step; k <= nbrPnt && k >= 1; k += step)
 446   {
 447     pk        = pnt2d(k);
 448     gp_Pnt P1 = points(k);
 449     if (myP3d[indMin].SquareDistance(P1) > prec2 && myP3d[indMin].SquareDistance(Value(pk)) > prec2)
 450       break;
 451   }
 452
 453   //: p8 abv 11 Mar 99: PRO7226 #489490: if whole pcurve is degenerate, distribute evenly
 454   if (k < 1 || k > nbrPnt)
 455   {
 456     Standard_Real x1 = (myUIsoDeg[indMin] ? pnt2d(1).Y() : pnt2d(1).X());
 457     Standard_Real x2 = (myUIsoDeg[indMin] ? pnt2d(nbrPnt).Y() : pnt2d(nbrPnt).X());
 458     for (j = 1; j <= nbrPnt; j++)
 459     {
 460       // szv#4:S4163:12Mar99 warning - possible div by zero
 461       Standard_Real x = (x1 * (nbrPnt - j) + x2 * (j - 1)) / (nbrPnt - 1);
 462       if (!myUIsoDeg[indMin])
 463         pnt2d(j).SetX(x);
 464       else
 465         pnt2d(j).SetY(x);
 466     }
 467     return Standard_True;
 468   }
 469
 470   for (j = k - step; j <= nbrPnt && j >= 1; j -= step)
 471   {
 472     if (!myUIsoDeg[indMin])
 473       pnt2d(j).SetX(pk.X());
 474     else
 475       pnt2d(j).SetY(pk.Y());
 476   }
 477   return Standard_True;
 478 }
 479
 480 //=======================================================================
 481 // method : IsDegenerated
 482 // purpose:
 483 //=======================================================================
 484
 485 Standard_Boolean ShapeAnalysis_Surface::IsDegenerated(const gp_Pnt2d&     p2d1,
 486                                                       const gp_Pnt2d&     p2d2,
 487                                                       const Standard_Real tol,
 488                                                       const Standard_Real ratio)
 489 {
 490   gp_Pnt        p1    = Value(p2d1);
 491   gp_Pnt        p2    = Value(p2d2);
 492   gp_Pnt        pm    = Value(0.5 * (p2d1.XY() + p2d2.XY()));
 493   Standard_Real max3d = Max(p1.Distance(p2), Max(pm.Distance(p1), pm.Distance(p2)));
 494   if (max3d > tol)
 495     return Standard_False;
 496
 497   GeomAdaptor_Surface& SA = *Adaptor3d();
 498   Standard_Real        RU = SA.UResolution(1.);
 499   Standard_Real        RV = SA.VResolution(1.);
 500
 501   if (RU < Precision::PConfusion() || RV < Precision::PConfusion())
 502     return 0;
 503   Standard_Real du = Abs(p2d1.X() - p2d2.X()) / RU;
 504   Standard_Real dv = Abs(p2d1.Y() - p2d2.Y()) / RV;
 505   max3d *= ratio;
 506   return du * du + dv * dv > max3d * max3d;
 507 }
 508
 509 //=======================================================================
 510 // static : ComputeIso
 511 // purpose  :
 512 //=======================================================================
 513
 514 static Handle(Geom_Curve) ComputeIso(const Handle(Geom_Surface)& surf,
 515                                      const Standard_Boolean      utype,
 516                                      const Standard_Real         par)
 517 {
 518   Handle(Geom_Curve) iso;
 519   try
 520   {
 521     OCC_CATCH_SIGNALS
 522     if (utype)
 523       iso = surf->UIso(par);
 524     else
 525       iso = surf->VIso(par);
 526   }
 527   catch (Standard_Failure const& anException)
 528   {
 529 #ifdef OCCT_DEBUG
 530     //: s5
 531     std::cout << "\nWarning: ShapeAnalysis_Surface, ComputeIso(): Exception in UVIso(): ";
 532     anException.Print(std::cout);
 533     std::cout << std::endl;
 534 #endif
 535     (void)anException;
 536     iso.Nullify();
 537   }
 538   return iso;
 539 }
 540
 541 //=================================================================================================
 542
 543 void ShapeAnalysis_Surface::ComputeBoundIsos()
 544 {
 545   if (myIsos)
 546     return;
 547   myIsos  = Standard_True;
 548   myIsoUF = ComputeIso(mySurf, Standard_True, myUF);
 549   myIsoUL = ComputeIso(mySurf, Standard_True, myUL);
 550   myIsoVF = ComputeIso(mySurf, Standard_False, myVF);
 551   myIsoVL = ComputeIso(mySurf, Standard_False, myVL);
 552 }
 553
 554 //=================================================================================================
 555
 556 Handle(Geom_Curve) ShapeAnalysis_Surface::UIso(const Standard_Real U)
 557 {
 558   if (U == myUF)
 559   {
 560     ComputeBoundIsos();
 561     return myIsoUF;
 562   }
 563   if (U == myUL)
 564   {
 565     ComputeBoundIsos();
 566     return myIsoUL;
 567   }
 568   return ComputeIso(mySurf, Standard_True, U);
 569 }
 570
 571 //=================================================================================================
 572
 573 Handle(Geom_Curve) ShapeAnalysis_Surface::VIso(const Standard_Real V)
 574 {
 575   if (V == myVF)
 576   {
 577     ComputeBoundIsos();
 578     return myIsoVF;
 579   }
 580   if (V == myVL)
 581   {
 582     ComputeBoundIsos();
 583     return myIsoVL;
 584   }
 585   return ComputeIso(mySurf, Standard_False, V);
 586 }
 587
 588 //=================================================================================================
 589
 590 Standard_Boolean ShapeAnalysis_Surface::IsUClosed(const Standard_Real preci)
 591 {
 592   Standard_Real prec      = Max(preci, Precision::Confusion());
 593   Standard_Real anUmidVal = -1.;
 594   if (myUCloseVal < 0)
 595   {
 596     //    Faut calculer : calculs minimaux
 597     Standard_Real uf, ul, vf, vl;
 598     Bounds(uf, ul, vf, vl); // modified by rln on 12/11/97 mySurf-> is deleted
 599     // mySurf->Bounds (uf,ul,vf,vl);
 600     RestrictBounds(uf, ul, vf, vl);
 601     myUDelt = Abs(ul - uf) / 20; // modified by rln 11/11/97 instead of 10
 602                                  // because of the example when 10 was not enough
 603     if (mySurf->IsUClosed())
 604     {
 605       myUCloseVal = 0.;
 606       myUDelt     = 0.;
 607       myGap       = 0.;
 608       return Standard_True;
 609     }
 610
 611     // Calculs adaptes
 612     // #67 rln S4135
 613     GeomAdaptor_Surface& SurfAdapt = *Adaptor3d();
 614     GeomAbs_SurfaceType  surftype  = SurfAdapt.GetType();
 615     if (mySurf->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)))
 616     {
 617       surftype = GeomAbs_OtherSurface;
 618     }
 619
 620     switch (surftype)
 621     {
 622       case GeomAbs_Plane: {
 623         myUCloseVal = RealLast();
 624         break;
 625       }
 626       case GeomAbs_SurfaceOfExtrusion: { //: c8 abv 03 Mar 98: UKI60094 #753: process
 627                                          //: Geom_SurfaceOfLinearExtrusion
 628         Handle(Geom_SurfaceOfLinearExtrusion) extr =
 629           Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(mySurf);
 630         Handle(Geom_Curve) crv = extr->BasisCurve();
 631         Standard_Real      f   = crv->FirstParameter();
 632         Standard_Real      l   = crv->LastParameter();
 633         //: r3 abv (smh) 30 Mar 99: protect against unexpected signals
 634         if (!Precision::IsInfinite(f) && !Precision::IsInfinite(l))
 635         {
 636           gp_Pnt p1   = crv->Value(f);
 637           gp_Pnt p2   = crv->Value(l);
 638           myUCloseVal = p1.SquareDistance(p2);
 639           gp_Pnt pm   = crv->Value((f + l) / 2.);
 640           anUmidVal   = p1.SquareDistance(pm);
 641         }
 642         else
 643         {
 644           myUCloseVal = RealLast();
 645         }
 646         break;
 647       }
 648       case GeomAbs_BSplineSurface: {
 649         Handle(Geom_BSplineSurface) bs      = Handle(Geom_BSplineSurface)::DownCast(mySurf);
 650         Standard_Integer            nbup    = bs->NbUPoles();
 651         Standard_Real               distmin = RealLast();
 652         if (bs->IsUPeriodic())
 653         {
 654           myUCloseVal = 0;
 655           myUDelt     = 0;
 656         }
 657         else if (nbup < 3)
 658         { // modified by rln on 12/11/97
 659           myUCloseVal = RealLast();
 660         }
 661         else if (bs->IsURational() ||
 662                  // #6 rln 20/02/98 ProSTEP ug_exhaust-A.stp entity #18360 (Uclosed BSpline,
 663                  // but multiplicity of boundary knots != degree + 1)
 664                  bs->UMultiplicity(1) != bs->UDegree() + 1 || // #6 //:h4: #6 moved
 665                  bs->UMultiplicity(bs->NbUKnots()) != bs->UDegree() + 1)
 666         { // #6 //:h4
 667           Standard_Integer nbvk = bs->NbVKnots();
 668           Standard_Real    v    = bs->VKnot(1);
 669           gp_Pnt           p1   = SurfAdapt.Value(uf, v);
 670           gp_Pnt           p2   = SurfAdapt.Value(ul, v);
 671           myUCloseVal           = p1.SquareDistance(p2);
 672           gp_Pnt pm             = SurfAdapt.Value((uf + ul) / 2., v);
 673           anUmidVal             = p1.SquareDistance(pm);
 674           distmin               = myUCloseVal;
 675           for (Standard_Integer i = 2; i <= nbvk; i++)
 676           {
 677             v                   = 0.5 * (bs->VKnot(i - 1) + bs->VKnot(i));
 678             p1                  = bs->Value(uf, v);
 679             p2                  = bs->Value(ul, v);
 680             Standard_Real aDist = p1.SquareDistance(p2);
 681             if (aDist > myUCloseVal)
 682             {
 683               myUCloseVal = aDist;
 684               pm          = bs->Value((uf + ul) / 2., v);
 685               anUmidVal   = p1.SquareDistance(pm);
 686             }
 687             else
 688             {
 689               distmin = Min(distmin, aDist);
 690             }
 691           }
 692           distmin = Sqrt(distmin);
 693           myUDelt = Min(myUDelt, 0.5 * SurfAdapt.UResolution(distmin)); // #4 smh
 694         }
 695         else
 696         {
 697           Standard_Integer nbvp = bs->NbVPoles();
 698           myUCloseVal           = bs->Pole(1, 1).SquareDistance(bs->Pole(nbup, 1));
 699           anUmidVal             = bs->Pole(1, 1).SquareDistance(bs->Pole(nbup / 2 + 1, 1));
 700           distmin               = myUCloseVal;
 701           for (Standard_Integer i = 2; i <= nbvp; i++)
 702           {
 703             Standard_Real aDist = bs->Pole(1, i).SquareDistance(bs->Pole(nbup, i));
 704             if (aDist > myUCloseVal)
 705             {
 706               myUCloseVal = aDist;
 707               anUmidVal   = bs->Pole(1, i).SquareDistance(bs->Pole(nbup / 2 + 1, i));
 708             }
 709             else
 710             {
 711               distmin = Min(distmin, aDist);
 712             }
 713           }
 714           distmin = Sqrt(distmin);
 715           myUDelt = Min(myUDelt, 0.5 * SurfAdapt.UResolution(distmin)); // #4 smh
 716         }
 717         break;
 718       }
 719       case GeomAbs_BezierSurface: {
 720         Handle(Geom_BezierSurface) bz      = Handle(Geom_BezierSurface)::DownCast(mySurf);
 721         Standard_Integer           nbup    = bz->NbUPoles();
 722         Standard_Real              distmin = RealLast();
 723         if (nbup < 3)
 724         {
 725           myUCloseVal = RealLast();
 726         }
 727         else
 728         {
 729           Standard_Integer nbvp = bz->NbVPoles();
 730           myUCloseVal           = bz->Pole(1, 1).SquareDistance(bz->Pole(nbup, 1));
 731           anUmidVal             = bz->Pole(1, 1).SquareDistance(bz->Pole(nbup / 2 + 1, 1));
 732           distmin               = myUCloseVal;
 733           for (Standard_Integer i = 1; i <= nbvp; i++)
 734           {
 735             Standard_Real aDist = bz->Pole(1, i).SquareDistance(bz->Pole(nbup, i));
 736             if (aDist > myUCloseVal)
 737             {
 738               myUCloseVal = aDist;
 739               anUmidVal   = bz->Pole(1, i).SquareDistance(bz->Pole(nbup / 2 + 1, i));
 740             }
 741             else
 742             {
 743               distmin = Min(distmin, aDist);
 744             }
 745           }
 746           distmin = Sqrt(distmin);
 747           myUDelt = Min(myUDelt, 0.5 * SurfAdapt.UResolution(distmin)); // #4 smh
 748         }
 749         break;
 750       }
 751       default: { // Geom_RectangularTrimmedSurface and Geom_OffsetSurface
 752         Standard_Real    distmin  = RealLast();
 753         Standard_Integer nbpoints = 101; // can be revised
 754         gp_Pnt           p1       = SurfAdapt.Value(uf, vf);
 755         gp_Pnt           p2       = SurfAdapt.Value(ul, vf);
 756         myUCloseVal               = p1.SquareDistance(p2);
 757         gp_Pnt pm                 = SurfAdapt.Value((uf + ul) / 2, vf);
 758         anUmidVal                 = p1.SquareDistance(pm);
 759         distmin                   = myUCloseVal;
 760         for (Standard_Integer i = 1; i < nbpoints; i++)
 761         {
 762           Standard_Real vparam = vf + (vl - vf) * i / (nbpoints - 1);
 763           p1                   = SurfAdapt.Value(uf, vparam);
 764           p2                   = SurfAdapt.Value(ul, vparam);
 765           Standard_Real aDist  = p1.SquareDistance(p2);
 766           if (aDist > myUCloseVal)
 767           {
 768             myUCloseVal = aDist;
 769             pm          = SurfAdapt.Value((uf + ul) / 2, vparam);
 770             anUmidVal   = p1.SquareDistance(pm);
 771           }
 772           else
 773           {
 774             distmin = Min(distmin, aDist);
 775           }
 776         }
 777         distmin = Sqrt(distmin);
 778         myUDelt = Min(myUDelt, 0.5 * SurfAdapt.UResolution(distmin)); // #4 smh
 779         break;
 780       }
 781     } // switch
 782     myGap       = sqrt(myUCloseVal);
 783     myUCloseVal = myGap;
 784   }
 785
 786   if (anUmidVal > 0. && myUCloseVal > sqrt(anUmidVal))
 787   {
 788     myUCloseVal = RealLast();
 789     return Standard_False;
 790   }
 791
 792   return (myUCloseVal <= prec);
 793 }
 794
 795 //=================================================================================================
 796
 797 Standard_Boolean ShapeAnalysis_Surface::IsVClosed(const Standard_Real preci)
 798 {
 799   Standard_Real prec     = Max(preci, Precision::Confusion());
 800   Standard_Real aVmidVal = -1.;
 801   if (myVCloseVal < 0)
 802   {
 803     //    Faut calculer : calculs minimaux
 804     Standard_Real uf, ul, vf, vl;
 805     Bounds(uf, ul, vf, vl); // modified by rln on 12/11/97 mySurf-> is deleted
 806                             //     mySurf->Bounds (uf,ul,vf,vl);
 807     RestrictBounds(uf, ul, vf, vl);
 808     myVDelt = Abs(vl - vf) / 20; // 2; rln S4135
 809                                  // because of the example when 10 was not enough
 810     if (mySurf->IsVClosed())
 811     {
 812       myVCloseVal = 0.;
 813       myVDelt     = 0.;
 814       myGap       = 0.;
 815       return Standard_True;
 816     }
 817
 818     //    Calculs adaptes
 819     // #67 rln S4135
 820     GeomAdaptor_Surface& SurfAdapt = *Adaptor3d();
 821     GeomAbs_SurfaceType  surftype  = SurfAdapt.GetType();
 822     if (mySurf->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)))
 823     {
 824       surftype = GeomAbs_OtherSurface;
 825     }
 826
 827     switch (surftype)
 828     {
 829       case GeomAbs_Plane:
 830       case GeomAbs_Cone:
 831       case GeomAbs_Cylinder:
 832       case GeomAbs_Sphere:
 833       case GeomAbs_SurfaceOfExtrusion: {
 834         myVCloseVal = RealLast();
 835         break;
 836       }
 837       case GeomAbs_SurfaceOfRevolution: {
 838         Handle(Geom_SurfaceOfRevolution) revol = Handle(Geom_SurfaceOfRevolution)::DownCast(mySurf);
 839         Handle(Geom_Curve)               crv   = revol->BasisCurve();
 840         gp_Pnt                           p1    = crv->Value(crv->FirstParameter());
 841         gp_Pnt                           p2    = crv->Value(crv->LastParameter());
 842         myVCloseVal                            = p1.SquareDistance(p2);
 843         break;
 844       }
 845       case GeomAbs_BSplineSurface: {
 846         Handle(Geom_BSplineSurface) bs      = Handle(Geom_BSplineSurface)::DownCast(mySurf);
 847         Standard_Integer            nbvp    = bs->NbVPoles();
 848         Standard_Real               distmin = RealLast();
 849         if (bs->IsVPeriodic())
 850         {
 851           myVCloseVal = 0;
 852           myVDelt     = 0;
 853         }
 854         else if (nbvp < 3)
 855         { // modified by rln on 12/11/97
 856           myVCloseVal = RealLast();
 857         }
 858         else if (bs->IsVRational() || bs->VMultiplicity(1) != bs->VDegree() + 1 || // #6 //:h4
 859                  bs->VMultiplicity(bs->NbVKnots()) != bs->VDegree() + 1)
 860         { // #6 //:h4
 861           Standard_Integer nbuk = bs->NbUKnots();
 862           Standard_Real    u    = bs->UKnot(1);
 863           gp_Pnt           p1   = SurfAdapt.Value(u, vf);
 864           gp_Pnt           p2   = SurfAdapt.Value(u, vl);
 865           myVCloseVal           = p1.SquareDistance(p2);
 866           gp_Pnt pm             = SurfAdapt.Value(u, (vf + vl) / 2.);
 867           aVmidVal              = p1.SquareDistance(pm);
 868           distmin               = myVCloseVal;
 869           for (Standard_Integer i = 2; i <= nbuk; i++)
 870           {
 871             u                   = 0.5 * (bs->UKnot(i - 1) + bs->UKnot(i));
 872             p1                  = SurfAdapt.Value(u, vf);
 873             p2                  = SurfAdapt.Value(u, vl);
 874             Standard_Real aDist = p1.SquareDistance(p2);
 875             if (aDist > myVCloseVal)
 876             {
 877               myVCloseVal = aDist;
 878               pm          = SurfAdapt.Value(u, (vf + vl) / 2);
 879               aVmidVal    = p1.SquareDistance(pm);
 880             }
 881             else
 882             {
 883               distmin = Min(distmin, aDist);
 884             }
 885           }
 886           distmin = Sqrt(distmin);
 887           myVDelt = Min(myVDelt, 0.5 * SurfAdapt.VResolution(distmin)); // #4 smh
 888         }
 889         else
 890         {
 891           Standard_Integer nbup = bs->NbUPoles();
 892           myVCloseVal           = bs->Pole(1, 1).SquareDistance(bs->Pole(1, nbvp));
 893           aVmidVal              = bs->Pole(1, 1).SquareDistance(bs->Pole(1, nbvp / 2 + 1));
 894           distmin               = myVCloseVal;
 895           for (Standard_Integer i = 2; i <= nbup; i++)
 896           {
 897             Standard_Real aDist = bs->Pole(i, 1).SquareDistance(bs->Pole(i, nbvp));
 898             if (aDist > myVCloseVal)
 899             {
 900               myVCloseVal = aDist;
 901               aVmidVal    = bs->Pole(i, 1).SquareDistance(bs->Pole(i, nbvp / 2 + 1));
 902             }
 903             else
 904             {
 905               distmin = Min(distmin, aDist);
 906             }
 907           }
 908           distmin = Sqrt(distmin);
 909           myVDelt = Min(myVDelt, 0.5 * SurfAdapt.VResolution(distmin)); // #4 smh
 910         }
 911         break;
 912       }
 913       case GeomAbs_BezierSurface: {
 914         Handle(Geom_BezierSurface) bz      = Handle(Geom_BezierSurface)::DownCast(mySurf);
 915         Standard_Integer           nbvp    = bz->NbVPoles();
 916         Standard_Real              distmin = RealLast();
 917         if (nbvp < 3)
 918         {
 919           myVCloseVal = RealLast();
 920         }
 921         else
 922         {
 923           Standard_Integer nbup = bz->NbUPoles();
 924           myVCloseVal           = bz->Pole(1, 1).SquareDistance(bz->Pole(1, nbvp));
 925           aVmidVal              = bz->Pole(1, 1).SquareDistance(bz->Pole(1, nbvp / 2 + 1));
 926           distmin               = myVCloseVal;
 927           for (Standard_Integer i = 2; i <= nbup; i++)
 928           {
 929             Standard_Real aDist = bz->Pole(i, 1).SquareDistance(bz->Pole(i, nbvp));
 930             if (aDist > myVCloseVal)
 931             {
 932               myVCloseVal = aDist;
 933               aVmidVal    = bz->Pole(i, 1).SquareDistance(bz->Pole(i, nbvp / 2 + 1));
 934             }
 935             else
 936             {
 937               distmin = Min(distmin, aDist);
 938             }
 939           }
 940           distmin = Sqrt(distmin);
 941           myVDelt = Min(myVDelt, 0.5 * SurfAdapt.VResolution(distmin)); // #4 smh
 942         }
 943         break;
 944       }
 945       default: { // Geom_RectangularTrimmedSurface and Geom_OffsetSurface
 946         Standard_Real    distmin  = RealLast();
 947         Standard_Integer nbpoints = 101; // can be revised
 948         gp_Pnt           p1       = SurfAdapt.Value(uf, vf);
 949         gp_Pnt           p2       = SurfAdapt.Value(uf, vl);
 950         gp_Pnt           pm       = SurfAdapt.Value(uf, (vf + vl) / 2);
 951         myVCloseVal               = p1.SquareDistance(p2);
 952         aVmidVal                  = p1.SquareDistance(pm);
 953         distmin                   = myVCloseVal;
 954         for (Standard_Integer i = 1; i < nbpoints; i++)
 955         {
 956           Standard_Real uparam = uf + (ul - uf) * i / (nbpoints - 1);
 957           p1                   = SurfAdapt.Value(uparam, vf);
 958           p2                   = SurfAdapt.Value(uparam, vl);
 959           Standard_Real aDist  = p1.SquareDistance(p2);
 960           if (aDist > myVCloseVal)
 961           {
 962             myVCloseVal = aDist;
 963             pm          = SurfAdapt.Value(uparam, (vf + vl) / 2);
 964             aVmidVal    = p1.SquareDistance(pm);
 965           }
 966           else
 967           {
 968             distmin = Min(distmin, aDist);
 969           }
 970         }
 971         distmin = Sqrt(distmin);
 972         myVDelt = Min(myVDelt, 0.5 * SurfAdapt.VResolution(distmin)); // #4 smh
 973         break;
 974       }
 975     } // switch
 976     myGap       = Sqrt(myVCloseVal);
 977     myVCloseVal = myGap;
 978   }
 979
 980   if (aVmidVal > 0. && myVCloseVal > sqrt(aVmidVal))
 981   {
 982     myVCloseVal = RealLast();
 983     return Standard_False;
 984   }
 985
 986   return (myVCloseVal <= prec);
 987 }
 988
 989 //=======================================================================
 990 // function : SurfaceNewton
 991 // purpose  : Newton algo (S4030)
 992 //=======================================================================
 993 Standard_Integer ShapeAnalysis_Surface::SurfaceNewton(const gp_Pnt2d&     p2dPrev,
 994                                                       const gp_Pnt&       P3D,
 995                                                       const Standard_Real preci,
 996                                                       gp_Pnt2d&           sol)
 997 {
 998   GeomAdaptor_Surface& SurfAdapt = *Adaptor3d();
 999   Standard_Real        uf, ul, vf, vl;
1000   Bounds(uf, ul, vf, vl);
1001   Standard_Real du = SurfAdapt.UResolution(preci);
1002   Standard_Real dv = SurfAdapt.VResolution(preci);
1003   Standard_Real UF = uf - du, UL = ul + du;
1004   Standard_Real VF = vf - dv, VL = vl + dv;
1005
1006   // Standard_Integer fail = 0;
1007   constexpr Standard_Real Tol  = Precision::Confusion();
1008   constexpr Standard_Real Tol2 = Tol * Tol; //, rs2p=1e10;
1009   Standard_Real           U = p2dPrev.X(), V = p2dPrev.Y();
1010   gp_Vec                  rsfirst = P3D.XYZ() - Value(U, V).XYZ(); // pdn
1011   for (Standard_Integer i = 0; i < 25; i++)
1012   {
1013     gp_Vec ru, rv, ruu, rvv, ruv;
1014     gp_Pnt pnt;
1015     SurfAdapt.D2(U, V, pnt, ru, rv, ruu, rvv, ruv);
1016
1017     // normal
1018     Standard_Real ru2 = ru * ru, rv2 = rv * rv;
1019     gp_Vec        n    = ru ^ rv;
1020     Standard_Real nrm2 = n.SquareMagnitude();
1021     if (nrm2 < 1e-10 || Precision::IsPositiveInfinite(nrm2))
1022       break; // n == 0, use standard
1023
1024     // descriminant
1025     gp_Vec        rs   = P3D.XYZ() - Value(U, V).XYZ();
1026     Standard_Real rSuu = (rs * ruu);
1027     Standard_Real rSvv = (rs * rvv);
1028     Standard_Real rSuv = (rs * ruv);
1029     Standard_Real D =
1030       -nrm2 + rv2 * rSuu + ru2 * rSvv - 2 * rSuv * (ru * rv) + rSuv * rSuv - rSuu * rSvv;
1031     if (fabs(D) < 1e-10)
1032       break; // bad case; use standard
1033
1034     // compute step
1035     Standard_Real fract = 1. / D;
1036     du                  = (rs * ((n ^ rv) + ru * rSvv - rv * rSuv)) * fract;
1037     dv                  = (rs * ((ru ^ n) + rv * rSuu - ru * rSuv)) * fract;
1038     U += du;
1039     V += dv;
1040     if (U < UF || U > UL || V < VF || V > VL)
1041       break;
1042     // check that iterations do not diverge
1043     // pdn    Standard_Real rs2 = rs.SquareMagnitude();
1044     //    if ( rs2 > 4.*rs2p ) break;
1045     //    rs2p = rs2;
1046
1047     // test the step by uv and deviation from the solution
1048     Standard_Real aResolution = Max(1e-12, (U + V) * 10e-16);
1049     if (fabs(du) + fabs(dv) > aResolution)
1050       continue; // Precision::PConfusion()  continue;
1051
1052     // if ( U < UF || U > UL || V < VF || V > VL ) break;
1053
1054     // pdn PRO10109 4517: protect against wrong result
1055     Standard_Real rs2 = rs.SquareMagnitude();
1056     if (rs2 > rsfirst.SquareMagnitude())
1057       break;
1058
1059     Standard_Real rsn = rs * n;
1060     if (rs2 - rsn * rsn / nrm2 > Tol2)
1061       break;
1062
1063     //  if ( rs2 > 100 * preci * preci ) { fail = 6; break; }
1064
1065     // OK, return the result
1066     //  std::cout << "Newton: solution found in " << i+1 << " iterations" << std::endl;
1067     sol.SetCoord(U, V);
1068
1069     return (nrm2 < 0.01 * ru2 * rv2 ? 2 : 1); //: q6
1070   }
1071   //      std::cout << "Newton: failed after " << i+1 << " iterations (fail " << fail << " )" <<
1072   //      std::endl;
1073   return Standard_False;
1074 }
1075
1076 //=======================================================================
1077 // function : NextValueOfUV
1078 // purpose  : optimizing projection by Newton algo (S4030)
1079 //=======================================================================
1080
1081 gp_Pnt2d ShapeAnalysis_Surface::NextValueOfUV(const gp_Pnt2d&     p2dPrev,
1082                                               const gp_Pnt&       P3D,
1083                                               const Standard_Real preci,
1084                                               const Standard_Real maxpreci)
1085 {
1086   GeomAdaptor_Surface& SurfAdapt = *Adaptor3d();
1087   GeomAbs_SurfaceType  surftype  = SurfAdapt.GetType();
1088
1089   switch (surftype)
1090   {
1091     case GeomAbs_BezierSurface:
1092     case GeomAbs_BSplineSurface:
1093     case GeomAbs_SurfaceOfExtrusion:
1094     case GeomAbs_SurfaceOfRevolution:
1095     case GeomAbs_OffsetSurface:
1096
1097     {
1098       if (surftype == GeomAbs_BSplineSurface)
1099       {
1100         Handle(Geom_BSplineSurface) aBSpline = SurfAdapt.BSpline();
1101
1102         // Check near to knot position ~ near to C0 points on U isoline.
1103         if (SurfAdapt.UContinuity() == GeomAbs_C0)
1104         {
1105           Standard_Integer aMinIndex = aBSpline->FirstUKnotIndex();
1106           Standard_Integer aMaxIndex = aBSpline->LastUKnotIndex();
1107           for (Standard_Integer anIdx = aMinIndex; anIdx <= aMaxIndex; ++anIdx)
1108           {
1109             Standard_Real aKnot = aBSpline->UKnot(anIdx);
1110             if (Abs(aKnot - p2dPrev.X()) < Precision::Confusion())
1111               return ValueOfUV(P3D, preci);
1112           }
1113         }
1114
1115         // Check near to knot position ~ near to C0 points on U isoline.
1116         if (SurfAdapt.VContinuity() == GeomAbs_C0)
1117         {
1118           Standard_Integer aMinIndex = aBSpline->FirstVKnotIndex();
1119           Standard_Integer aMaxIndex = aBSpline->LastVKnotIndex();
1120           for (Standard_Integer anIdx = aMinIndex; anIdx <= aMaxIndex; ++anIdx)
1121           {
1122             Standard_Real aKnot = aBSpline->VKnot(anIdx);
1123             if (Abs(aKnot - p2dPrev.Y()) < Precision::Confusion())
1124               return ValueOfUV(P3D, preci);
1125           }
1126         }
1127       }
1128
1129       gp_Pnt2d         sol;
1130       Standard_Integer res = SurfaceNewton(p2dPrev, P3D, preci, sol);
1131       if (res != 0)
1132       {
1133         Standard_Real gap = P3D.Distance(Value(sol));
1134         if (res == 2 || //: q6 abv 19 Mar 99: protect against strange attractors
1135             (maxpreci > 0. && gap - maxpreci > Precision::Confusion()))
1136         { //: q1: check with maxpreci
1137           Standard_Real U = sol.X(), V = sol.Y();
1138           myGap = UVFromIso(P3D, preci, U, V);
1139           //      gp_Pnt2d p = ValueOfUV ( P3D, preci );
1140           if (gap >= myGap)
1141             return gp_Pnt2d(U, V);
1142         }
1143         myGap = gap;
1144         return sol;
1145       }
1146     }
1147     break;
1148     default:
1149       break;
1150   }
1151   return ValueOfUV(P3D, preci);
1152 }
1153
1154 //=================================================================================================
1155
1156 gp_Pnt2d ShapeAnalysis_Surface::ValueOfUV(const gp_Pnt& P3D, const Standard_Real preci)
1157 {
1158   GeomAdaptor_Surface& SurfAdapt = *Adaptor3d();
1159   Standard_Real        S = 0., T = 0.;
1160   myGap                     = -1.;           // devra etre calcule
1161   Standard_Boolean computed = Standard_True; // a priori
1162
1163   Standard_Real uf, ul, vf, vl;
1164   Bounds(uf, ul, vf, vl); // modified by rln on 12/11/97 mySurf-> is deleted
1165
1166   { //: c9 abv 3 Mar 98: UKI60107-1 #350: to prevent 'catch' from catching exception raising below
1167     //: it
1168     try
1169     { // ajout CKY 30-DEC-1997 (cf ProStep TR6 r_89-ug)
1170       OCC_CATCH_SIGNALS
1171       GeomAbs_SurfaceType surftype = SurfAdapt.GetType();
1172       switch (surftype)
1173       {
1174
1175         case GeomAbs_Plane: {
1176           gp_Pln Plane = SurfAdapt.Plane();
1177           ElSLib::Parameters(Plane, P3D, S, T);
1178           break;
1179         }
1180         case GeomAbs_Cylinder: {
1181           gp_Cylinder Cylinder = SurfAdapt.Cylinder();
1182           ElSLib::Parameters(Cylinder, P3D, S, T);
1183           S += ShapeAnalysis::AdjustByPeriod(S, 0.5 * (uf + ul), 2 * M_PI);
1184           break;
1185         }
1186         case GeomAbs_Cone: {
1187           gp_Cone Cone = SurfAdapt.Cone();
1188           ElSLib::Parameters(Cone, P3D, S, T);
1189           S += ShapeAnalysis::AdjustByPeriod(S, 0.5 * (uf + ul), 2 * M_PI);
1190           break;
1191         }
1192         case GeomAbs_Sphere: {
1193           gp_Sphere Sphere = SurfAdapt.Sphere();
1194           ElSLib::Parameters(Sphere, P3D, S, T);
1195           S += ShapeAnalysis::AdjustByPeriod(S, 0.5 * (uf + ul), 2 * M_PI);
1196           break;
1197         }
1198         case GeomAbs_Torus: {
1199           gp_Torus Torus = SurfAdapt.Torus();
1200           ElSLib::Parameters(Torus, P3D, S, T);
1201           S += ShapeAnalysis::AdjustByPeriod(S, 0.5 * (uf + ul), 2 * M_PI);
1202           T += ShapeAnalysis::AdjustByPeriod(T, 0.5 * (vf + vl), 2 * M_PI);
1203           break;
1204         }
1205         case GeomAbs_BezierSurface:
1206         case GeomAbs_BSplineSurface:
1207         case GeomAbs_SurfaceOfExtrusion:
1208         case GeomAbs_SurfaceOfRevolution:
1209         case GeomAbs_OffsetSurface: //: d0 abv 3 Mar 98: UKI60107-1 #350
1210         {
1211           S = (uf + ul) / 2;
1212           T = (vf + vl) / 2; // yaura aumoins qqchose
1213                              // pdn to fix hangs PRO17015
1214           if ((surftype == GeomAbs_SurfaceOfExtrusion) && Precision::IsInfinite(uf)
1215               && Precision::IsInfinite(ul))
1216           {
1217             // conic case
1218             gp_Pnt2d prev(S, T);
1219             gp_Pnt2d solution;
1220             if (SurfaceNewton(prev, P3D, preci, solution) != 0)
1221             {
1222 #ifdef OCCT_DEBUG
1223               std::cout << "Newton found point on conic extrusion" << std::endl;
1224 #endif
1225               return solution;
1226             }
1227 #ifdef OCCT_DEBUG
1228             std::cout << "Newton failed point on conic extrusion" << std::endl;
1229 #endif
1230             uf = -500;
1231             ul = 500;
1232           }
1233
1234           RestrictBounds(uf, ul, vf, vl);
1235
1236           //: 30 by abv 2.12.97: speed optimization
1237           // code is taken from GeomAPI_ProjectPointOnSurf
1238           if (!myExtOK)
1239           {
1240             //      Standard_Real du = Abs(ul-uf)/100;  Standard_Real dv = Abs(vl-vf)/100;
1241             //      if (IsUClosed()) du = 0;  if (IsVClosed()) dv = 0;
1242             //  Forcer appel a IsU-VClosed
1243             if (myUCloseVal < 0)
1244               IsUClosed();
1245             if (myVCloseVal < 0)
1246               IsVClosed();
1247             Standard_Real du = 0., dv = 0.;
1248             // extension of the surface range is limited to non-offset surfaces as the latter
1249             // can throw exception (e.g. Geom_UndefinedValue) when computing value - see id23943
1250             if (!mySurf->IsKind(STANDARD_TYPE(Geom_OffsetSurface)))
1251             {
1252               // modified by rln during fixing CSR # BUC60035 entity #D231
1253               du = Min(myUDelt, SurfAdapt.UResolution(preci));
1254               dv = Min(myVDelt, SurfAdapt.VResolution(preci));
1255             }
1256             constexpr Standard_Real Tol = Precision::PConfusion();
1257             myExtPS.SetFlag(Extrema_ExtFlag_MIN);
1258             myExtPS.Initialize(SurfAdapt, uf - du, ul + du, vf - dv, vl + dv, Tol, Tol);
1259             myExtOK = Standard_True;
1260           }
1261           myExtPS.Perform(P3D);
1262           Standard_Integer nPSurf = (myExtPS.IsDone() ? myExtPS.NbExt() : 0);
1263
1264           if (nPSurf > 0)
1265           {
1266             Standard_Real    dist2Min = myExtPS.SquareDistance(1);
1267             Standard_Integer indMin   = 1;
1268             for (Standard_Integer sol = 2; sol <= nPSurf; sol++)
1269             {
1270               Standard_Real dist2 = myExtPS.SquareDistance(sol);
1271               if (dist2Min > dist2)
1272               {
1273                 dist2Min = dist2;
1274                 indMin   = sol;
1275               }
1276             }
1277             myExtPS.Point(indMin).Parameter(S, T);
1278             // PTV 26.06.2002 WORKAROUND protect OCC486. Remove after fix bug.
1279             // file CEA_cuve-V5.igs Entityes 244, 259, 847, 925
1280             // if project point3D on SurfaceOfRevolution Extreme recompute 2d point, but
1281             // returns an old distance from 3d to solution :-(
1282             gp_Pnt aCheckPnt = SurfAdapt.Value(S, T);
1283             dist2Min         = P3D.SquareDistance(aCheckPnt);
1284             // end of WORKAROUND
1285             Standard_Real disSurf = sqrt(dist2Min); //, disCurv =1.e10;
1286
1287             // Test de projection merdeuse sur les bords :
1288             Standard_Real UU = S, VV = T,
1289                           DistMinOnIso =
1290                             RealLast(); // myGap;
1291                                         //      ForgetNewton(P3D, mySurf, preci, UU, VV, DistMinOnIso);
1292
1293             // test added by rln on 08/12/97
1294             //  DistMinOnIso = UVFromIso (P3D, preci, UU, VV);
1295             Standard_Boolean possLockal = Standard_False; //: study S4030 (optimizing)
1296             if (disSurf > preci)
1297             {
1298               gp_Pnt2d pp(UU, VV);
1299               if (SurfaceNewton(pp, P3D, preci, pp) != 0)
1300               { //: q2 abv 16 Mar 99: PRO7226 #412920
1301                 Standard_Real dist = P3D.Distance(Value(pp));
1302                 if (dist < disSurf)
1303                 {
1304                   disSurf = dist;
1305                   S = UU = pp.X();
1306                   T = VV = pp.Y();
1307                 }
1308               }
1309               if (disSurf < 10 * preci)
1310                 if (mySurf->Continuity() != GeomAbs_C0)
1311                 {
1312                   constexpr Standard_Real Tol = Precision::Confusion();
1313                   gp_Vec                  D1U, D1V;
1314                   gp_Pnt                  pnt;
1315                   SurfAdapt.D1(UU, VV, pnt, D1U, D1V);
1316                   gp_Vec        b = D1U.Crossed(D1V);
1317                   gp_Vec        a(pnt, P3D);
1318                   Standard_Real ab   = a.Dot(b);
1319                   Standard_Real nrm2 = b.SquareMagnitude();
1320                   if (nrm2 > 1e-10)
1321                   {
1322                     Standard_Real dist = a.SquareMagnitude() - (ab * ab) / nrm2;
1323                     possLockal         = (dist < Tol * Tol);
1324                   }
1325                 }
1326               if (!possLockal)
1327               {
1328                 DistMinOnIso = UVFromIso(P3D, preci, UU, VV);
1329               }
1330             }
1331
1332             if (disSurf > DistMinOnIso)
1333             {
1334               // On prend les parametres UU et VV;
1335               S     = UU;
1336               T     = VV;
1337               myGap = DistMinOnIso;
1338             }
1339             else
1340             {
1341               myGap = disSurf;
1342             }
1343
1344             // On essaie Intersection Droite Passant par P3D / Surface
1345             //  if ((myGap > preci)&&(!possLockal) ) {
1346             //    Standard_Real SS, TT;
1347             //    disCurv = FindUV(P3D, mySurf, S, T, SS, TT);
1348             //    if (disCurv < preci || disCurv < myGap) {
1349             //      S = SS;
1350             //      T = TT;
1351             //    }
1352             //  }
1353           }
1354           else
1355           {
1356 #ifdef OCCT_DEBUG
1357             std::cout << "Warning: ShapeAnalysis_Surface::ValueOfUV(): Extrema failed, doing Newton"
1358                       << std::endl;
1359 #endif
1360             // on essai sur les bords
1361             Standard_Real UU = S,
1362                           VV = T; //, DistMinOnIso;
1363                                   //    ForgetNewton(P3D, mySurf, preci, UU, VV, DistMinOnIso);
1364             myGap = UVFromIso(P3D, preci, UU, VV);
1365             //  if (DistMinOnIso > preci) {
1366             //    Standard_Real SS, TT;
1367             //    Standard_Real disCurv = FindUV(P3D, mySurf, UU, VV, SS, TT);
1368             //    if (disCurv < preci) {
1369             //      S = SS;
1370             //      T = TT;
1371             //    }
1372             //  }
1373             //  else {
1374             S = UU;
1375             T = VV;
1376             //  }
1377           }
1378         }
1379         break;
1380
1381         default:
1382           computed = Standard_False;
1383           break;
1384       }
1385
1386     } // end Try ValueOfUV (CKY 30-DEC-1997)
1387
1388     catch (Standard_Failure const& anException)
1389     {
1390 #ifdef OCCT_DEBUG
1391       //   Pas de raison mais qui sait. Mieux vaut retourner un truc faux que stopper
1392       //   L ideal serait d avoir un status ... mais qui va l interroger ?
1393       //   Avec ce status, on saurait que ce point est a sauter et voila tout
1394       //   En attendant, on met une valeur "pas idiote" mais surement fausse ...
1395       // szv#4:S4163:12Mar99 optimized
1396       //: s5
1397       std::cout << "\nWarning: ShapeAnalysis_Surface::ValueOfUV(): Exception: ";
1398       anException.Print(std::cout);
1399       std::cout << std::endl;
1400 #endif
1401       (void)anException;
1402       S = (Precision::IsInfinite(uf)) ? 0 : (uf + ul) / 2.;
1403       T = (Precision::IsInfinite(vf)) ? 0 : (vf + vl) / 2.;
1404     }
1405   } //: c9
1406     // szv#4:S4163:12Mar99 waste raise
1407     // if (!computed) throw Standard_NoSuchObject("PCurveLib_ProjectPointOnSurf::ValueOfUV untreated
1408     // surface type");
1409   if (computed)
1410   {
1411     if (myGap <= 0)
1412       myGap = P3D.Distance(SurfAdapt.Value(S, T));
1413   }
1414   else
1415   {
1416     myGap = -1.;
1417     S     = 0.;
1418     T     = 0.;
1419   }
1420   return gp_Pnt2d(S, T);
1421 }
1422
1423 //=================================================================================================
1424
1425 Standard_Real ShapeAnalysis_Surface::UVFromIso(const gp_Pnt&       P3d,
1426                                                const Standard_Real preci,
1427                                                Standard_Real&      U,
1428                                                Standard_Real&      V)
1429 {
1430   //  Projection qui considere les isos ... comme suit :
1431   //  Les 4 bords, plus les isos en U et en V
1432   //  En effet, souvent, un des deux est bon ...
1433   Standard_Real theMin = RealLast();
1434
1435   gp_Pnt        pntres;
1436   Standard_Real Cf, Cl, UU, VV;
1437
1438   //  Initialisation des recherches : point deja trouve (?)
1439   UU            = U;
1440   VV            = V;
1441   gp_Pnt depart = myAdSur->Value(U, V);
1442   theMin        = depart.Distance(P3d);
1443
1444   if (theMin < preci / 10)
1445     return theMin; // c etait deja OK
1446   ComputeBoxes();
1447   if (myIsoUF.IsNull() || myIsoUL.IsNull() || myIsoVF.IsNull() || myIsoVL.IsNull())
1448   {
1449     // no isolines
1450     // no more precise computation
1451     return theMin;
1452   }
1453   try
1454   { // RAJOUT
1455     OCC_CATCH_SIGNALS
1456     // pdn Create BndBox containing point;
1457     Bnd_Box aPBox;
1458     aPBox.Set(P3d);
1459
1460     // std::cout<<"Adaptor3d()->Surface().GetType() =
1461     // "<<Adaptor3d()->Surface().GetType()<<std::endl;
1462
1463     // modified by rln on 04/12/97 in order to use these variables later
1464     Standard_Boolean   UV  = Standard_True;
1465     Standard_Real      par = 0., other = 0., dist = 0.;
1466     Handle(Geom_Curve) iso;
1467     Adaptor3d_IsoCurve anIsoCurve(Adaptor3d());
1468     for (Standard_Integer num = 0; num < 6; num++)
1469     {
1470
1471       UV = (num < 3); // 0-1-2 : iso-U  3-4-5 : iso-V
1472       if (!(Adaptor3d()->GetType() == GeomAbs_OffsetSurface))
1473       {
1474         const Bnd_Box* anIsoBox = 0;
1475         switch (num)
1476         {
1477           case 0:
1478             par      = myUF;
1479             iso      = myIsoUF;
1480             anIsoBox = &myBndUF;
1481             break;
1482           case 1:
1483             par      = myUL;
1484             iso      = myIsoUL;
1485             anIsoBox = &myBndUL;
1486             break;
1487           case 2:
1488             par = U;
1489             iso = UIso(U);
1490             break;
1491           case 3:
1492             par      = myVF;
1493             iso      = myIsoVF;
1494             anIsoBox = &myBndVF;
1495             break;
1496           case 4:
1497             par      = myVL;
1498             iso      = myIsoVL;
1499             anIsoBox = &myBndVL;
1500             break;
1501           case 5:
1502             par = V;
1503             iso = VIso(V);
1504             break;
1505           default:
1506             break;
1507         }
1508
1509         //    On y va la-dessus
1510         if (!Precision::IsInfinite(par) && !iso.IsNull())
1511         {
1512           if (anIsoBox && anIsoBox->Distance(aPBox) > theMin)
1513             continue;
1514
1515           Cf = iso->FirstParameter();
1516           Cl = iso->LastParameter();
1517
1518           RestrictBounds(Cf, Cl);
1519           dist =
1520             ShapeAnalysis_Curve().Project(iso, P3d, preci, pntres, other, Cf, Cl, Standard_False);
1521           if (dist < theMin)
1522           {
1523             theMin = dist;
1524             //: q6      if (UV) VV = other;  else UU = other;
1525             //  Selon une isoU, on calcule le meilleur V; et lycee de Versailles
1526             UU = (UV ? par : other);
1527             VV = (UV ? other : par); //: q6: uncommented
1528           }
1529         }
1530       }
1531
1532       else
1533       {
1534         Adaptor3d_Curve*  anAdaptor = NULL;
1535         GeomAdaptor_Curve aGeomCurve;
1536
1537         const Bnd_Box* anIsoBox = 0;
1538         switch (num)
1539         {
1540           case 0:
1541             par = myUF;
1542             aGeomCurve.Load(myIsoUF);
1543             anAdaptor = &aGeomCurve;
1544             anIsoBox  = &myBndUF;
1545             break;
1546           case 1:
1547             par = myUL;
1548             aGeomCurve.Load(myIsoUL);
1549             anAdaptor = &aGeomCurve;
1550             anIsoBox  = &myBndUL;
1551             break;
1552           case 2:
1553             par = U;
1554             anIsoCurve.Load(GeomAbs_IsoU, U);
1555             anAdaptor = &anIsoCurve;
1556             break;
1557           case 3:
1558             par = myVF;
1559             aGeomCurve.Load(myIsoVF);
1560             anAdaptor = &aGeomCurve;
1561             anIsoBox  = &myBndVF;
1562             break;
1563           case 4:
1564             par = myVL;
1565             aGeomCurve.Load(myIsoVL);
1566             anAdaptor = &aGeomCurve;
1567             anIsoBox  = &myBndVL;
1568             break;
1569           case 5:
1570             par = V;
1571             anIsoCurve.Load(GeomAbs_IsoV, V);
1572             anAdaptor = &anIsoCurve;
1573             break;
1574           default:
1575             break;
1576         }
1577         if (anIsoBox && anIsoBox->Distance(aPBox) > theMin)
1578           continue;
1579         dist = ShapeAnalysis_Curve().Project(*anAdaptor, P3d, preci, pntres, other);
1580         if (dist < theMin)
1581         {
1582           theMin = dist;
1583           UU     = (UV ? par : other);
1584           VV     = (UV ? other : par); //: q6: uncommented
1585         }
1586       }
1587     }
1588
1589     // added by rln on 04/12/97 iterational process
1590     Standard_Real    PrevU = U, PrevV = V;
1591     Standard_Integer MaxIters = 5, Iters = 0;
1592     if (!(Adaptor3d()->GetType() == GeomAbs_OffsetSurface))
1593     {
1594       while (((PrevU != UU) || (PrevV != VV)) && (Iters < MaxIters) && (theMin > preci))
1595       {
1596         PrevU = UU;
1597         PrevV = VV;
1598         if (UV)
1599         {
1600           par = UU;
1601           iso = UIso(UU);
1602         }
1603         else
1604         {
1605           par = VV;
1606           iso = VIso(VV);
1607         }
1608         if (!iso.IsNull())
1609         {
1610           Cf = iso->FirstParameter();
1611           Cl = iso->LastParameter();
1612           RestrictBounds(Cf, Cl);
1613           dist =
1614             ShapeAnalysis_Curve().Project(iso, P3d, preci, pntres, other, Cf, Cl, Standard_False);
1615           if (dist < theMin)
1616           {
1617             theMin = dist;
1618             if (UV)
1619               VV = other;
1620             else
1621               UU = other;
1622           }
1623         }
1624         UV = !UV;
1625         if (UV)
1626         {
1627           par = UU;
1628           iso = UIso(UU);
1629         }
1630         else
1631         {
1632           par = VV;
1633           iso = VIso(VV);
1634         }
1635         if (!iso.IsNull())
1636         {
1637           Cf = iso->FirstParameter();
1638           Cl = iso->LastParameter();
1639           RestrictBounds(Cf, Cl);
1640           dist =
1641             ShapeAnalysis_Curve().Project(iso, P3d, preci, pntres, other, Cf, Cl, Standard_False);
1642           if (dist < theMin)
1643           {
1644             theMin = dist;
1645             if (UV)
1646               VV = other;
1647             else
1648               UU = other;
1649           }
1650         }
1651         UV = !UV;
1652         Iters++;
1653       }
1654     }
1655
1656     else
1657     {
1658       while (((PrevU != UU) || (PrevV != VV)) && (Iters < MaxIters) && (theMin > preci))
1659       {
1660         PrevU = UU;
1661         PrevV = VV;
1662         if (UV)
1663         {
1664           par = UU;
1665           anIsoCurve.Load(GeomAbs_IsoU, UU);
1666         }
1667         else
1668         {
1669           par = VV;
1670           anIsoCurve.Load(GeomAbs_IsoV, VV);
1671         }
1672         Cf = anIsoCurve.FirstParameter();
1673         Cl = anIsoCurve.LastParameter();
1674         RestrictBounds(Cf, Cl);
1675         dist = ShapeAnalysis_Curve().Project(anIsoCurve, P3d, preci, pntres, other);
1676         if (dist < theMin)
1677         {
1678           theMin = dist;
1679           if (UV)
1680             VV = other;
1681           else
1682             UU = other;
1683         }
1684         UV = !UV;
1685         if (UV)
1686         {
1687           par = UU;
1688           anIsoCurve.Load(GeomAbs_IsoU, UU);
1689         }
1690         else
1691         {
1692           par = VV;
1693           anIsoCurve.Load(GeomAbs_IsoV, VV);
1694         }
1695         Cf = anIsoCurve.FirstParameter();
1696         Cl = anIsoCurve.LastParameter();
1697         RestrictBounds(Cf, Cl);
1698         dist = ShapeAnalysis_Curve().ProjectAct(anIsoCurve, P3d, preci, pntres, other);
1699         if (dist < theMin)
1700         {
1701           theMin = dist;
1702           if (UV)
1703             VV = other;
1704           else
1705             UU = other;
1706         }
1707         UV = !UV;
1708         Iters++;
1709       }
1710     }
1711
1712     U = UU;
1713     V = VV;
1714
1715   } // fin try RAJOUT
1716   catch (Standard_Failure const& anException)
1717   {
1718 #ifdef OCCT_DEBUG
1719     //: s5
1720     std::cout << "\nWarning: ShapeAnalysis_Curve::UVFromIso(): Exception: ";
1721     anException.Print(std::cout);
1722     std::cout << std::endl;
1723 #endif
1724     (void)anException;
1725     theMin = RealLast(); // theMin de depart
1726   }
1727   return theMin;
1728 }
1729
1730 //=================================================================================================
1731
1732 void ShapeAnalysis_Surface::SortSingularities()
1733 {
1734   for (Standard_Integer i = 0; i < myNbDeg - 1; i++)
1735   {
1736     Standard_Real    minPreci = myPreci[i];
1737     Standard_Integer minIndex = i;
1738     for (Standard_Integer j = i + 1; j < myNbDeg; j++)
1739       if (minPreci > myPreci[j])
1740       {
1741         minPreci = myPreci[j];
1742         minIndex = j;
1743       }
1744     if (minIndex != i)
1745     {
1746       myPreci[minIndex]           = myPreci[i];
1747       myPreci[i]                  = minPreci;
1748       gp_Pnt tmpP3d               = myP3d[minIndex];
1749       myP3d[minIndex]             = myP3d[i];
1750       myP3d[i]                    = tmpP3d;
1751       gp_Pnt2d tmpP2d             = myFirstP2d[minIndex];
1752       myFirstP2d[minIndex]        = myFirstP2d[i];
1753       myFirstP2d[i]               = tmpP2d;
1754       tmpP2d                      = myLastP2d[minIndex];
1755       myLastP2d[minIndex]         = myLastP2d[i];
1756       myLastP2d[i]                = tmpP2d;
1757       Standard_Real tmpPar        = myFirstPar[minIndex];
1758       myFirstPar[minIndex]        = myFirstPar[i];
1759       myFirstPar[i]               = tmpPar;
1760       tmpPar                      = myLastPar[minIndex];
1761       myLastPar[minIndex]         = myLastPar[i];
1762       myLastPar[i]                = tmpPar;
1763       Standard_Boolean tmpUIsoDeg = myUIsoDeg[minIndex];
1764       myUIsoDeg[minIndex]         = myUIsoDeg[i];
1765       myUIsoDeg[i]                = tmpUIsoDeg;
1766     }
1767   }
1768 }
1769
1770 //=================================================================================================
1771
1772 void ShapeAnalysis_Surface::SetDomain(const Standard_Real U1,
1773                                       const Standard_Real U2,
1774                                       const Standard_Real V1,
1775                                       const Standard_Real V2)
1776 {
1777   myUF = U1;
1778   myUL = U2;
1779   myVF = V1;
1780   myVL = V2;
1781 }
1782
1783 void ShapeAnalysis_Surface::ComputeBoxes()
1784 {
1785   if (myIsoBoxes)
1786     return;
1787   myIsoBoxes = Standard_True;
1788   ComputeBoundIsos();
1789   if (!myIsoUF.IsNull())
1790     BndLib_Add3dCurve::Add(GeomAdaptor_Curve(myIsoUF), Precision::Confusion(), myBndUF);
1791   if (!myIsoUL.IsNull())
1792     BndLib_Add3dCurve::Add(GeomAdaptor_Curve(myIsoUL), Precision::Confusion(), myBndUL);
1793   if (!myIsoVF.IsNull())
1794     BndLib_Add3dCurve::Add(GeomAdaptor_Curve(myIsoVF), Precision::Confusion(), myBndVF);
1795   if (!myIsoVL.IsNull())
1796     BndLib_Add3dCurve::Add(GeomAdaptor_Curve(myIsoVL), Precision::Confusion(), myBndVL);
1797 }
1798
1799 const Bnd_Box& ShapeAnalysis_Surface::GetBoxUF()
1800 {
1801   ComputeBoxes();
1802   return myBndUF;
1803 }
1804
1805 const Bnd_Box& ShapeAnalysis_Surface::GetBoxUL()
1806 {
1807   ComputeBoxes();
1808   return myBndUL;
1809 }
1810
1811 const Bnd_Box& ShapeAnalysis_Surface::GetBoxVF()
1812 {
1813   ComputeBoxes();
1814   return myBndVF;
1815 }
1816
1817 const Bnd_Box& ShapeAnalysis_Surface::GetBoxVL()
1818 {
1819   ComputeBoxes();
1820   return myBndVL;
1821 }