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