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