src/ShapeAnalysis/ShapeAnalysis_Surface.cxx

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