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