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