src/ShapeAnalysis/ShapeAnalysis_Surface.cxx

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