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