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