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