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