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