0024008: ShapeAnalysis_Surface causes exception in Geom_OffsetSurface
[occt.git] / src / ShapeAnalysis / ShapeAnalysis_Surface.cxx
CommitLineData
b311480e 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
7fd59977 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
70ShapeAnalysis_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
86void 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
104void 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
119const 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
131void 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);
c6541a0c
D
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);
7fd59977 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
239Standard_Boolean ShapeAnalysis_Surface::HasSingularities (const Standard_Real preci)
240{
241 return NbSingularities(preci) > 0;
242}
243
244//=======================================================================
245//function : NbSingularities
246//purpose :
247//=======================================================================
248
249Standard_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
291Standard_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
308Standard_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
345Standard_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
379Standard_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
440Standard_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
468static 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
493void 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
508Handle(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
520Handle(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
532Standard_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
662Standard_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//=======================================================================
787Standard_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;
2e9fd4bc 831 if ( U < UF || U > UL || V < VF || V > VL ) break;
7fd59977 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
2e9fd4bc 841 //if ( U < UF || U > UL || V < VF || V > VL ) break;
7fd59977 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
867gp_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
911gp_Pnt2d ShapeAnalysis_Surface::ValueOfUV(const gp_Pnt& P3D,const Standard_Real preci)
912{
913 GeomAdaptor_Surface& SurfAdapt = Adaptor3d()->ChangeSurface();
1d47d8d0 914 Standard_Real S = 0., T = 0.;
7fd59977 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);
c6541a0c 937 S += ShapeAnalysis::AdjustByPeriod(S,0.5*(uf+ul),2*M_PI);
7fd59977 938 break;
939 }
940 case GeomAbs_Cone :
941 {
942 gp_Cone Cone = SurfAdapt.Cone();
943 ElSLib::Parameters( Cone, P3D, S, T);
c6541a0c 944 S += ShapeAnalysis::AdjustByPeriod(S,0.5*(uf+ul),2*M_PI);
7fd59977 945 break;
946 }
947 case GeomAbs_Sphere :
948 {
949 gp_Sphere Sphere = SurfAdapt.Sphere();
950 ElSLib::Parameters( Sphere, P3D, S, T);
c6541a0c 951 S += ShapeAnalysis::AdjustByPeriod(S,0.5*(uf+ul),2*M_PI);
7fd59977 952 break;
953 }
954 case GeomAbs_Torus :
955 {
956 gp_Torus Torus = SurfAdapt.Torus();
957 ElSLib::Parameters( Torus, P3D, S, T);
c6541a0c
D
958 S += ShapeAnalysis::AdjustByPeriod(S,0.5*(uf+ul),2*M_PI);
959 T += ShapeAnalysis::AdjustByPeriod(T,0.5*(vf+vl),2*M_PI);
7fd59977 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();
c15398ab
RL
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 }
7fd59977 1008 myExtSrf = mySurf;
1009 Standard_Real Tol = Precision::PConfusion();
6062bf3e 1010 myExtPS.SetFlag (Extrema_ExtFlag_MIN);
7fd59977 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
1152Standard_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
1346void 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
1375void 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
1387void 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
1402const Bnd_Box& ShapeAnalysis_Surface::GetBoxUF()
1403{
1404 ComputeBoxes();
1405 return myBndUF;
1406}
1407
1408const Bnd_Box& ShapeAnalysis_Surface::GetBoxUL()
1409{
1410 ComputeBoxes();
1411 return myBndUL;
1412}
1413
1414const Bnd_Box& ShapeAnalysis_Surface::GetBoxVF()
1415{
1416 ComputeBoxes();
1417 return myBndVF;
1418}
1419
1420const Bnd_Box& ShapeAnalysis_Surface::GetBoxVL()
1421{
1422 ComputeBoxes();
1423 return myBndVL;
1424}