b311480e |
1 | // Copyright (c) 1995-1999 Matra Datavision |
973c2be1 |
2 | // Copyright (c) 1999-2014 OPEN CASCADE SAS |
b311480e |
3 | // |
973c2be1 |
4 | // This file is part of Open CASCADE Technology software library. |
b311480e |
5 | // |
d5f74e42 |
6 | // This library is free software; you can redistribute it and/or modify it under |
7 | // the terms of the GNU Lesser General Public License version 2.1 as published |
973c2be1 |
8 | // by the Free Software Foundation, with special exception defined in the file |
9 | // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT |
10 | // distribution for complete text of the license and disclaimer of any warranty. |
b311480e |
11 | // |
973c2be1 |
12 | // Alternatively, this file may be used under the terms of Open CASCADE |
13 | // commercial license or contractual agreement. |
b311480e |
14 | |
7fd59977 |
15 | #include <LProp_Status.hxx> |
16 | #include <LProp_NotDefined.hxx> |
17 | #include <Standard_OutOfRange.hxx> |
18 | #include <Standard_DomainError.hxx> |
19 | #include <CSLib.hxx> |
20 | #include <CSLib_DerivativeStatus.hxx> |
21 | #include <CSLib_NormalStatus.hxx> |
22 | #include <TColgp_Array2OfVec.hxx> |
23 | #include <math_DirectPolynomialRoots.hxx> |
24 | |
32ca7a51 |
25 | static const Standard_Real MinStep = 1.0e-7; |
7fd59977 |
26 | |
27 | static Standard_Boolean IsTangentDefined (LProp_SLProps& SProp, |
32ca7a51 |
28 | const Standard_Integer cn, |
29 | const Standard_Real linTol, |
30 | const Standard_Integer Derivative, |
31 | Standard_Integer& Order, |
32 | LProp_Status& Status) |
7fd59977 |
33 | { |
34 | Standard_Real Tol = linTol * linTol; |
35 | gp_Vec V[2]; |
36 | Order = 0; |
32ca7a51 |
37 | |
38 | while (Order < 3) |
39 | { |
7fd59977 |
40 | Order++; |
32ca7a51 |
41 | if(cn >= Order) |
42 | { |
43 | switch(Order) |
44 | { |
45 | case 1: |
46 | V[0] = SProp.D1U(); |
47 | V[1] = SProp.D1V(); |
48 | break; |
49 | case 2: |
50 | V[0] = SProp.D2U(); |
51 | V[1] = SProp.D2V(); |
52 | break; |
53 | }//switch(Order) |
54 | |
55 | if(V[Derivative].SquareMagnitude() > Tol) |
56 | { |
57 | Status = LProp_Defined; |
58 | return Standard_True; |
7fd59977 |
59 | } |
32ca7a51 |
60 | }//if(cn >= Order) |
61 | else |
62 | { |
7fd59977 |
63 | Status = LProp_Undefined; |
64 | return Standard_False; |
65 | } |
66 | } |
32ca7a51 |
67 | |
7fd59977 |
68 | return Standard_False; |
69 | } |
70 | |
71 | LProp_SLProps::LProp_SLProps (const Surface& S, |
32ca7a51 |
72 | const Standard_Real U, |
73 | const Standard_Real V, |
74 | const Standard_Integer N, |
75 | const Standard_Real Resolution) |
76 | : mySurf(S),myDerOrder(N), myCN(4), // (Tool::Continuity(S)), |
77 | myLinTol(Resolution) |
7fd59977 |
78 | { |
7fd59977 |
79 | Standard_OutOfRange_Raise_if(N < 0 || N > 2, |
32ca7a51 |
80 | "LProp_SLProps::LProp_SLProps()"); |
7fd59977 |
81 | |
82 | SetParameters(U, V); |
83 | } |
84 | |
85 | LProp_SLProps::LProp_SLProps (const Surface& S, |
86 | const Standard_Integer N, |
87 | const Standard_Real Resolution) |
32ca7a51 |
88 | : mySurf(S), myU(RealLast()), myV(RealLast()), myDerOrder(N), |
89 | myCN(4), // (Tool::Continuity(S)) |
90 | myLinTol(Resolution), |
91 | myUTangentStatus (LProp_Undecided), |
92 | myVTangentStatus (LProp_Undecided), |
93 | myNormalStatus (LProp_Undecided), |
94 | myCurvatureStatus(LProp_Undecided) |
7fd59977 |
95 | { |
96 | Standard_OutOfRange_Raise_if(N < 0 || N > 2, |
32ca7a51 |
97 | "LProp_SLProps::LProp_SLProps()"); |
7fd59977 |
98 | } |
99 | |
100 | LProp_SLProps::LProp_SLProps (const Standard_Integer N, |
32ca7a51 |
101 | const Standard_Real Resolution) |
102 | : myU(RealLast()), myV(RealLast()), myDerOrder(N), myCN(0), |
103 | myLinTol(Resolution), |
104 | myUTangentStatus (LProp_Undecided), |
105 | myVTangentStatus (LProp_Undecided), |
106 | myNormalStatus (LProp_Undecided), |
107 | myCurvatureStatus(LProp_Undecided) |
7fd59977 |
108 | { |
109 | Standard_OutOfRange_Raise_if(N < 0 || N > 2, |
32ca7a51 |
110 | "LProp_SLProps::LProp_SLProps() bad level"); |
7fd59977 |
111 | } |
112 | |
32ca7a51 |
113 | void LProp_SLProps::SetSurface (const Surface& S ) |
114 | { |
115 | mySurf = S; |
116 | myCN = 4; // =Tool::Continuity(S); |
7fd59977 |
117 | } |
118 | |
32ca7a51 |
119 | void LProp_SLProps::SetParameters (const Standard_Real U, |
120 | const Standard_Real V) |
7fd59977 |
121 | { |
32ca7a51 |
122 | myU = U; |
123 | myV = V; |
124 | switch (myDerOrder) |
125 | { |
7fd59977 |
126 | case 0: |
32ca7a51 |
127 | Tool::Value(mySurf, myU, myV, myPnt); |
7fd59977 |
128 | break; |
129 | case 1: |
32ca7a51 |
130 | Tool::D1(mySurf, myU, myV, myPnt, myD1u, myD1v); |
7fd59977 |
131 | break; |
132 | case 2: |
32ca7a51 |
133 | Tool::D2(mySurf, myU, myV, myPnt, myD1u, myD1v, myD2u, myD2v, myDuv); |
7fd59977 |
134 | break; |
32ca7a51 |
135 | } |
136 | |
137 | myUTangentStatus = LProp_Undecided; |
138 | myVTangentStatus = LProp_Undecided; |
139 | myNormalStatus = LProp_Undecided; |
140 | myCurvatureStatus = LProp_Undecided; |
7fd59977 |
141 | } |
142 | |
143 | const gp_Pnt& LProp_SLProps::Value() const |
144 | { |
32ca7a51 |
145 | return myPnt; |
7fd59977 |
146 | } |
147 | |
148 | const gp_Vec& LProp_SLProps::D1U() |
149 | { |
32ca7a51 |
150 | if (myDerOrder < 1) |
151 | { |
152 | myDerOrder =1; |
153 | Tool::D1(mySurf,myU,myV,myPnt,myD1u,myD1v); |
7fd59977 |
154 | } |
32ca7a51 |
155 | |
156 | return myD1u; |
7fd59977 |
157 | } |
158 | |
159 | const gp_Vec& LProp_SLProps::D1V() |
160 | { |
32ca7a51 |
161 | if (myDerOrder < 1) |
162 | { |
163 | myDerOrder =1; |
164 | Tool::D1(mySurf,myU,myV,myPnt,myD1u,myD1v); |
7fd59977 |
165 | } |
32ca7a51 |
166 | |
167 | return myD1v; |
7fd59977 |
168 | } |
169 | |
170 | const gp_Vec& LProp_SLProps::D2U() |
171 | { |
32ca7a51 |
172 | if (myDerOrder < 2) |
173 | { |
174 | myDerOrder =2; |
175 | Tool::D2(mySurf,myU,myV,myPnt,myD1u,myD1v,myD2u,myD2v,myDuv); |
7fd59977 |
176 | } |
32ca7a51 |
177 | |
178 | return myD2u; |
7fd59977 |
179 | } |
180 | |
181 | const gp_Vec& LProp_SLProps::D2V() |
182 | { |
32ca7a51 |
183 | if (myDerOrder < 2) |
184 | { |
185 | myDerOrder =2; |
186 | Tool::D2(mySurf,myU,myV,myPnt,myD1u,myD1v,myD2u,myD2v,myDuv); |
7fd59977 |
187 | } |
32ca7a51 |
188 | |
189 | return myD2v; |
7fd59977 |
190 | } |
32ca7a51 |
191 | |
7fd59977 |
192 | const gp_Vec& LProp_SLProps::DUV() |
193 | { |
32ca7a51 |
194 | if (myDerOrder < 2) |
195 | { |
196 | myDerOrder =2; |
197 | Tool::D2(mySurf,myU,myV,myPnt,myD1u,myD1v,myD2u,myD2v,myDuv); |
7fd59977 |
198 | } |
32ca7a51 |
199 | |
200 | return myDuv; |
7fd59977 |
201 | } |
202 | |
203 | Standard_Boolean LProp_SLProps::IsTangentUDefined () |
204 | { |
32ca7a51 |
205 | if (myUTangentStatus == LProp_Undefined) |
7fd59977 |
206 | return Standard_False; |
32ca7a51 |
207 | else if (myUTangentStatus >= LProp_Defined) |
7fd59977 |
208 | return Standard_True; |
32ca7a51 |
209 | |
7fd59977 |
210 | // uTangentStatus == Lprop_Undecided |
211 | // we have to calculate the first non null U derivative |
32ca7a51 |
212 | return IsTangentDefined(*this, myCN, myLinTol, 0, |
213 | mySignificantFirstDerivativeOrderU, myUTangentStatus); |
7fd59977 |
214 | } |
215 | |
32ca7a51 |
216 | void LProp_SLProps::TangentU (gp_Dir& D) |
217 | { |
218 | if(!IsTangentUDefined()) |
219 | LProp_NotDefined::Raise(); |
220 | |
221 | if(mySignificantFirstDerivativeOrderU == 1) |
222 | D = gp_Dir(myD1u); |
223 | else |
224 | { |
225 | const Standard_Real DivisionFactor = 1.e-3; |
226 | Standard_Real anUsupremum, anUinfium; |
227 | Standard_Real anVsupremum, anVinfium; |
228 | Tool::Bounds(mySurf,anUinfium,anVinfium,anUsupremum,anVsupremum); |
229 | Standard_Real du; |
230 | if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) |
231 | du = 0.0; |
232 | else |
233 | du = anUsupremum-anUinfium; |
234 | |
235 | const Standard_Real aDeltaU = Max(du*DivisionFactor,MinStep); |
236 | |
237 | gp_Vec V = myD2u; |
238 | |
239 | Standard_Real u; |
240 | |
241 | if(myU-anUinfium < aDeltaU) |
242 | u = myU+aDeltaU; |
243 | else |
244 | u = myU-aDeltaU; |
245 | |
246 | gp_Pnt P1, P2; |
247 | Tool::Value(mySurf, Min(myU, u),myV,P1); |
248 | Tool::Value(mySurf, Max(myU, u),myV,P2); |
249 | |
250 | gp_Vec V1(P1,P2); |
251 | Standard_Real aDirFactor = V.Dot(V1); |
252 | |
253 | if(aDirFactor < 0.0) |
254 | V = -V; |
255 | |
256 | D = gp_Dir(V); |
7fd59977 |
257 | } |
258 | } |
259 | |
260 | Standard_Boolean LProp_SLProps::IsTangentVDefined () |
261 | { |
32ca7a51 |
262 | if (myVTangentStatus == LProp_Undefined) |
7fd59977 |
263 | return Standard_False; |
32ca7a51 |
264 | else if (myVTangentStatus >= LProp_Defined) |
7fd59977 |
265 | return Standard_True; |
32ca7a51 |
266 | |
7fd59977 |
267 | // vTangentStatus == Lprop_Undecided |
268 | // we have to calculate the first non null V derivative |
32ca7a51 |
269 | return IsTangentDefined(*this, myCN, myLinTol, 1, |
270 | mySignificantFirstDerivativeOrderV, myVTangentStatus); |
7fd59977 |
271 | } |
272 | |
32ca7a51 |
273 | void LProp_SLProps::TangentV (gp_Dir& D) |
274 | { |
275 | if(!IsTangentVDefined()) |
276 | LProp_NotDefined::Raise(); |
277 | |
278 | if(mySignificantFirstDerivativeOrderV == 1) |
279 | D = gp_Dir(myD1v); |
280 | else |
281 | { |
282 | const Standard_Real DivisionFactor = 1.e-3; |
283 | Standard_Real anUsupremum, anUinfium; |
284 | Standard_Real anVsupremum, anVinfium; |
285 | Tool::Bounds(mySurf,anUinfium,anVinfium,anUsupremum,anVsupremum); |
286 | |
287 | Standard_Real dv; |
288 | if((anVsupremum >= RealLast()) || (anVinfium <= RealFirst())) |
289 | dv = 0.0; |
290 | else |
291 | dv = anVsupremum-anVinfium; |
292 | |
293 | const Standard_Real aDeltaV = Max(dv*DivisionFactor,MinStep); |
294 | |
295 | gp_Vec V = myD2v; |
296 | |
297 | Standard_Real v; |
298 | |
299 | if(myV-anVinfium < aDeltaV) |
300 | v = myV+aDeltaV; |
301 | else |
302 | v = myV-aDeltaV; |
303 | |
304 | gp_Pnt P1, P2; |
305 | Tool::Value(mySurf, myU, Min(myV, v),P1); |
306 | Tool::Value(mySurf, myU, Max(myV, v),P2); |
307 | |
308 | gp_Vec V1(P1,P2); |
309 | Standard_Real aDirFactor = V.Dot(V1); |
310 | |
311 | if(aDirFactor < 0.0) |
312 | V = -V; |
313 | |
314 | D = gp_Dir(V); |
7fd59977 |
315 | } |
316 | } |
317 | |
318 | Standard_Boolean LProp_SLProps::IsNormalDefined() |
319 | { |
32ca7a51 |
320 | if (myNormalStatus == LProp_Undefined) |
7fd59977 |
321 | return Standard_False; |
32ca7a51 |
322 | else if (myNormalStatus >= LProp_Defined) |
7fd59977 |
323 | return Standard_True; |
32ca7a51 |
324 | |
7fd59977 |
325 | // status = UnDecided |
7fd59977 |
326 | // first try the standard computation of the normal. |
327 | CSLib_DerivativeStatus Status; |
32ca7a51 |
328 | CSLib::Normal(myD1u, myD1v, myLinTol, Status, myNormal); |
329 | if (Status == CSLib_Done ) |
330 | { |
331 | myNormalStatus = LProp_Computed; |
7fd59977 |
332 | return Standard_True; |
333 | } |
7fd59977 |
334 | |
32ca7a51 |
335 | // else solve the degenerated case only if continuity >= 2 |
7fd59977 |
336 | |
32ca7a51 |
337 | myNormalStatus = LProp_Undefined; |
7fd59977 |
338 | return Standard_False; |
339 | } |
340 | |
32ca7a51 |
341 | const gp_Dir& LProp_SLProps::Normal () |
342 | { |
343 | if(!IsNormalDefined()) |
344 | { |
345 | LProp_NotDefined::Raise(); |
346 | } |
7fd59977 |
347 | |
32ca7a51 |
348 | return myNormal; |
7fd59977 |
349 | } |
350 | |
32ca7a51 |
351 | |
7fd59977 |
352 | Standard_Boolean LProp_SLProps::IsCurvatureDefined () |
353 | { |
32ca7a51 |
354 | if (myCurvatureStatus == LProp_Undefined) |
7fd59977 |
355 | return Standard_False; |
32ca7a51 |
356 | else if (myCurvatureStatus >= LProp_Defined) |
7fd59977 |
357 | return Standard_True; |
32ca7a51 |
358 | |
359 | if(myCN < 2) |
360 | { |
361 | myCurvatureStatus = LProp_Undefined; |
7fd59977 |
362 | return Standard_False; |
363 | } |
32ca7a51 |
364 | |
7fd59977 |
365 | // status = UnDecided |
32ca7a51 |
366 | if (!IsNormalDefined()) |
367 | { |
368 | myCurvatureStatus = LProp_Undefined; |
7fd59977 |
369 | return Standard_False; |
370 | } |
32ca7a51 |
371 | |
7fd59977 |
372 | // pour eviter un plantage dans le cas du caro pointu |
373 | // en fait on doit pouvoir calculer les courbure |
374 | // avoir |
32ca7a51 |
375 | if(!IsTangentUDefined() || !IsTangentVDefined()) |
376 | { |
377 | myCurvatureStatus = LProp_Undefined; |
7fd59977 |
378 | return Standard_False; |
379 | } |
380 | |
381 | // here we compute the curvature features of the surface |
382 | |
32ca7a51 |
383 | gp_Vec Norm (myNormal); |
7fd59977 |
384 | |
32ca7a51 |
385 | Standard_Real E = myD1u.SquareMagnitude(); |
386 | Standard_Real F = myD1u.Dot(myD1v); |
387 | Standard_Real G = myD1v.SquareMagnitude(); |
7fd59977 |
388 | |
32ca7a51 |
389 | if(myDerOrder < 2) |
390 | this->D2U(); |
7fd59977 |
391 | |
32ca7a51 |
392 | Standard_Real L = Norm.Dot(myD2u); |
393 | Standard_Real M = Norm.Dot(myDuv); |
394 | Standard_Real N = Norm.Dot(myD2v); |
7fd59977 |
395 | |
396 | Standard_Real A = E * M - F * L; |
397 | Standard_Real B = E * N - G * L; |
398 | Standard_Real C = F * N - G * M; |
399 | |
400 | Standard_Real MaxABC = Max(Max(Abs(A),Abs(B)),Abs(C)); |
32ca7a51 |
401 | if (MaxABC < RealEpsilon()) // ombilic |
402 | { |
403 | myMinCurv = N / G; |
404 | myMaxCurv = myMinCurv; |
405 | myDirMinCurv = gp_Dir (myD1u); |
406 | myDirMaxCurv = gp_Dir (myD1u.Crossed(Norm)); |
407 | myMeanCurv = myMinCurv; // (Cmin + Cmax) / 2. |
408 | myGausCurv = myMinCurv * myMinCurv; // (Cmin * Cmax) |
409 | myCurvatureStatus = LProp_Computed; |
7fd59977 |
410 | return Standard_True; |
411 | } |
412 | |
413 | A = A / MaxABC; |
414 | B = B / MaxABC; |
415 | C = C / MaxABC; |
416 | Standard_Real Curv1, Curv2, Root1, Root2; |
417 | gp_Vec VectCurv1, VectCurv2; |
32ca7a51 |
418 | |
419 | if (Abs(A) > RealEpsilon()) |
420 | { |
7fd59977 |
421 | math_DirectPolynomialRoots Root (A, B, C); |
32ca7a51 |
422 | if(Root.NbSolutions() != 2) |
423 | { |
424 | myCurvatureStatus = LProp_Undefined; |
7fd59977 |
425 | return Standard_False; |
426 | } |
32ca7a51 |
427 | else |
428 | { |
429 | Root1 = Root.Value(1); |
7fd59977 |
430 | Root2 = Root.Value(2); |
431 | Curv1 = ((L * Root1 + 2. * M) * Root1 + N) / |
432 | ((E * Root1 + 2. * F) * Root1 + G); |
433 | Curv2 = ((L * Root2 + 2. * M) * Root2 + N) / |
434 | ((E * Root2 + 2. * F) * Root2 + G); |
32ca7a51 |
435 | VectCurv1 = Root1 * myD1u + myD1v; |
436 | VectCurv2 = Root2 * myD1u + myD1v; |
7fd59977 |
437 | } |
438 | } |
32ca7a51 |
439 | else if (Abs(C) > RealEpsilon()) |
440 | { |
7fd59977 |
441 | math_DirectPolynomialRoots Root(C, B, A); |
32ca7a51 |
442 | |
443 | if((Root.NbSolutions() != 2)) |
444 | { |
445 | myCurvatureStatus = LProp_Undefined; |
7fd59977 |
446 | return Standard_False; |
447 | } |
32ca7a51 |
448 | else |
449 | { |
7fd59977 |
450 | Root1 = Root.Value(1); |
451 | Root2 = Root.Value(2); |
452 | Curv1 = ((N * Root1 + 2. * M) * Root1 + L) / |
453 | ((G * Root1 + 2. * F) * Root1 + E); |
454 | Curv2 = ((N * Root2 + 2. * M) * Root2 + L) / |
455 | ((G * Root2 + 2. * F) * Root2 + E); |
32ca7a51 |
456 | VectCurv1 = myD1u + Root1 * myD1v; |
457 | VectCurv2 = myD1u + Root2 * myD1v; |
7fd59977 |
458 | } |
459 | } |
32ca7a51 |
460 | else |
461 | { |
7fd59977 |
462 | Curv1 = L / E; |
463 | Curv2 = N / G; |
32ca7a51 |
464 | VectCurv1 = myD1u; |
465 | VectCurv2 = myD1v; |
7fd59977 |
466 | } |
32ca7a51 |
467 | |
468 | if (Curv1 < Curv2) |
469 | { |
470 | myMinCurv = Curv1; |
471 | myMaxCurv = Curv2; |
472 | myDirMinCurv = gp_Dir (VectCurv1); |
473 | myDirMaxCurv = gp_Dir (VectCurv2); |
7fd59977 |
474 | } |
32ca7a51 |
475 | else |
476 | { |
477 | myMinCurv = Curv2; |
478 | myMaxCurv = Curv1; |
479 | myDirMinCurv = gp_Dir (VectCurv2); |
480 | myDirMaxCurv = gp_Dir (VectCurv1); |
7fd59977 |
481 | } |
32ca7a51 |
482 | |
483 | myMeanCurv = ((N * E) - (2. * M * F) + (L * G)) // voir Farin p.282 |
7fd59977 |
484 | / (2. * ((E * G) - (F * F))); |
32ca7a51 |
485 | myGausCurv = ((L * N) - (M * M)) |
7fd59977 |
486 | / ((E * G) - (F * F)); |
32ca7a51 |
487 | myCurvatureStatus = LProp_Computed; |
7fd59977 |
488 | return Standard_True; |
32ca7a51 |
489 | } |
7fd59977 |
490 | |
491 | Standard_Boolean LProp_SLProps::IsUmbilic () |
492 | { |
32ca7a51 |
493 | if(!IsCurvatureDefined()) |
494 | LProp_NotDefined::Raise(); |
495 | |
496 | return Abs(myMaxCurv - myMinCurv) < Abs(Epsilon(myMaxCurv)); |
7fd59977 |
497 | } |
498 | |
499 | Standard_Real LProp_SLProps::MaxCurvature () |
500 | { |
32ca7a51 |
501 | if(!IsCurvatureDefined()) |
502 | LProp_NotDefined::Raise(); |
503 | |
504 | return myMaxCurv; |
7fd59977 |
505 | } |
506 | |
507 | Standard_Real LProp_SLProps::MinCurvature () |
508 | { |
32ca7a51 |
509 | if(!IsCurvatureDefined()) |
510 | LProp_NotDefined::Raise(); |
511 | |
512 | return myMinCurv; |
7fd59977 |
513 | } |
514 | |
515 | void LProp_SLProps::CurvatureDirections(gp_Dir& Max, gp_Dir& Min) |
516 | { |
32ca7a51 |
517 | if(!IsCurvatureDefined()) |
518 | LProp_NotDefined::Raise(); |
7fd59977 |
519 | |
32ca7a51 |
520 | Max = myDirMaxCurv; |
521 | Min = myDirMinCurv; |
7fd59977 |
522 | } |
523 | |
32ca7a51 |
524 | Standard_Real LProp_SLProps::MeanCurvature () |
525 | { |
526 | if(!IsCurvatureDefined()) |
527 | LProp_NotDefined::Raise(); |
7fd59977 |
528 | |
32ca7a51 |
529 | return myMeanCurv; |
7fd59977 |
530 | } |
531 | |
32ca7a51 |
532 | Standard_Real LProp_SLProps::GaussianCurvature () |
533 | { |
534 | if(!IsCurvatureDefined()) |
535 | LProp_NotDefined::Raise(); |
7fd59977 |
536 | |
32ca7a51 |
537 | return myGausCurv; |
7fd59977 |
538 | } |