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