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