1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
4 // This file is part of Open CASCADE Technology software library.
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
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.
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
15 #include <LProp_Status.hxx>
16 #include <LProp_NotDefined.hxx>
17 #include <Standard_OutOfRange.hxx>
18 #include <Standard_DomainError.hxx>
20 #include <CSLib_DerivativeStatus.hxx>
21 #include <CSLib_NormalStatus.hxx>
22 #include <TColgp_Array2OfVec.hxx>
23 #include <math_DirectPolynomialRoots.hxx>
25 static const Standard_Real MinStep = 1.0e-7;
27 static Standard_Boolean IsTangentDefined (LProp_SLProps& SProp,
28 const Standard_Integer cn,
29 const Standard_Real linTol,
30 const Standard_Integer Derivative,
31 Standard_Integer& Order,
32 LProp_Status& theStatus)
34 Standard_Real Tol = linTol * linTol;
55 if(V[Derivative].SquareMagnitude() > Tol)
57 theStatus = LProp_Defined;
63 theStatus = LProp_Undefined;
64 return Standard_False;
68 return Standard_False;
71 LProp_SLProps::LProp_SLProps (const Surface& S,
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)),
79 Standard_OutOfRange_Raise_if(N < 0 || N > 2,
80 "LProp_SLProps::LProp_SLProps()");
85 LProp_SLProps::LProp_SLProps (const Surface& S,
86 const Standard_Integer N,
87 const Standard_Real Resolution)
88 : mySurf(S), myU(RealLast()), myV(RealLast()), myDerOrder(N),
89 myCN(4), // (Tool::Continuity(S))
91 myUTangentStatus (LProp_Undecided),
92 myVTangentStatus (LProp_Undecided),
93 myNormalStatus (LProp_Undecided),
94 myCurvatureStatus(LProp_Undecided)
96 Standard_OutOfRange_Raise_if(N < 0 || N > 2,
97 "LProp_SLProps::LProp_SLProps()");
100 LProp_SLProps::LProp_SLProps (const Standard_Integer N,
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)
109 Standard_OutOfRange_Raise_if(N < 0 || N > 2,
110 "LProp_SLProps::LProp_SLProps() bad level");
113 void LProp_SLProps::SetSurface (const Surface& S )
116 myCN = 4; // =Tool::Continuity(S);
119 void LProp_SLProps::SetParameters (const Standard_Real U,
120 const Standard_Real V)
127 Tool::Value(mySurf, myU, myV, myPnt);
130 Tool::D1(mySurf, myU, myV, myPnt, myD1u, myD1v);
133 Tool::D2(mySurf, myU, myV, myPnt, myD1u, myD1v, myD2u, myD2v, myDuv);
137 myUTangentStatus = LProp_Undecided;
138 myVTangentStatus = LProp_Undecided;
139 myNormalStatus = LProp_Undecided;
140 myCurvatureStatus = LProp_Undecided;
143 const gp_Pnt& LProp_SLProps::Value() const
148 const gp_Vec& LProp_SLProps::D1U()
153 Tool::D1(mySurf,myU,myV,myPnt,myD1u,myD1v);
159 const gp_Vec& LProp_SLProps::D1V()
164 Tool::D1(mySurf,myU,myV,myPnt,myD1u,myD1v);
170 const gp_Vec& LProp_SLProps::D2U()
175 Tool::D2(mySurf,myU,myV,myPnt,myD1u,myD1v,myD2u,myD2v,myDuv);
181 const gp_Vec& LProp_SLProps::D2V()
186 Tool::D2(mySurf,myU,myV,myPnt,myD1u,myD1v,myD2u,myD2v,myDuv);
192 const gp_Vec& LProp_SLProps::DUV()
197 Tool::D2(mySurf,myU,myV,myPnt,myD1u,myD1v,myD2u,myD2v,myDuv);
203 Standard_Boolean LProp_SLProps::IsTangentUDefined ()
205 if (myUTangentStatus == LProp_Undefined)
206 return Standard_False;
207 else if (myUTangentStatus >= LProp_Defined)
208 return Standard_True;
210 // uTangentStatus == Lprop_Undecided
211 // we have to calculate the first non null U derivative
212 return IsTangentDefined(*this, myCN, myLinTol, 0,
213 mySignificantFirstDerivativeOrderU, myUTangentStatus);
216 void LProp_SLProps::TangentU (gp_Dir& D)
218 if(!IsTangentUDefined())
219 throw LProp_NotDefined();
221 if(mySignificantFirstDerivativeOrderU == 1)
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);
230 if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst()))
233 du = anUsupremum-anUinfium;
235 const Standard_Real aDeltaU = Max(du*DivisionFactor,MinStep);
241 if(myU-anUinfium < aDeltaU)
247 Tool::Value(mySurf, Min(myU, u),myV,P1);
248 Tool::Value(mySurf, Max(myU, u),myV,P2);
251 Standard_Real aDirFactor = V.Dot(V1);
260 Standard_Boolean LProp_SLProps::IsTangentVDefined ()
262 if (myVTangentStatus == LProp_Undefined)
263 return Standard_False;
264 else if (myVTangentStatus >= LProp_Defined)
265 return Standard_True;
267 // vTangentStatus == Lprop_Undecided
268 // we have to calculate the first non null V derivative
269 return IsTangentDefined(*this, myCN, myLinTol, 1,
270 mySignificantFirstDerivativeOrderV, myVTangentStatus);
273 void LProp_SLProps::TangentV (gp_Dir& D)
275 if(!IsTangentVDefined())
276 throw LProp_NotDefined();
278 if(mySignificantFirstDerivativeOrderV == 1)
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);
288 if((anVsupremum >= RealLast()) || (anVinfium <= RealFirst()))
291 dv = anVsupremum-anVinfium;
293 const Standard_Real aDeltaV = Max(dv*DivisionFactor,MinStep);
299 if(myV-anVinfium < aDeltaV)
305 Tool::Value(mySurf, myU, Min(myV, v),P1);
306 Tool::Value(mySurf, myU, Max(myV, v),P2);
309 Standard_Real aDirFactor = V.Dot(V1);
318 Standard_Boolean LProp_SLProps::IsNormalDefined()
320 if (myNormalStatus == LProp_Undefined)
321 return Standard_False;
322 else if (myNormalStatus >= LProp_Defined)
323 return Standard_True;
325 // status = UnDecided
326 // first try the standard computation of the normal.
327 CSLib_DerivativeStatus aStatus = CSLib_Done;
328 CSLib::Normal(myD1u, myD1v, myLinTol, aStatus, myNormal);
329 if (aStatus == CSLib_Done)
331 myNormalStatus = LProp_Computed;
332 return Standard_True;
335 // else solve the degenerated case only if continuity >= 2
337 myNormalStatus = LProp_Undefined;
338 return Standard_False;
341 const gp_Dir& LProp_SLProps::Normal ()
343 if(!IsNormalDefined())
345 throw LProp_NotDefined();
352 Standard_Boolean LProp_SLProps::IsCurvatureDefined ()
354 if (myCurvatureStatus == LProp_Undefined)
355 return Standard_False;
356 else if (myCurvatureStatus >= LProp_Defined)
357 return Standard_True;
361 myCurvatureStatus = LProp_Undefined;
362 return Standard_False;
365 // status = UnDecided
366 if (!IsNormalDefined())
368 myCurvatureStatus = LProp_Undefined;
369 return Standard_False;
372 // pour eviter un plantage dans le cas du caro pointu
373 // en fait on doit pouvoir calculer les courbure
375 if(!IsTangentUDefined() || !IsTangentVDefined())
377 myCurvatureStatus = LProp_Undefined;
378 return Standard_False;
381 // here we compute the curvature features of the surface
383 gp_Vec Norm (myNormal);
385 Standard_Real E = myD1u.SquareMagnitude();
386 Standard_Real F = myD1u.Dot(myD1v);
387 Standard_Real G = myD1v.SquareMagnitude();
392 Standard_Real L = Norm.Dot(myD2u);
393 Standard_Real M = Norm.Dot(myDuv);
394 Standard_Real N = Norm.Dot(myD2v);
396 Standard_Real A = E * M - F * L;
397 Standard_Real B = E * N - G * L;
398 Standard_Real C = F * N - G * M;
400 Standard_Real MaxABC = Max(Max(Abs(A),Abs(B)),Abs(C));
401 if (MaxABC < RealEpsilon()) // ombilic
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;
410 return Standard_True;
416 Standard_Real Curv1, Curv2, Root1, Root2;
417 gp_Vec VectCurv1, VectCurv2;
419 if (Abs(A) > RealEpsilon())
421 math_DirectPolynomialRoots Root (A, B, C);
422 if(Root.NbSolutions() != 2)
424 myCurvatureStatus = LProp_Undefined;
425 return Standard_False;
429 Root1 = Root.Value(1);
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);
435 VectCurv1 = Root1 * myD1u + myD1v;
436 VectCurv2 = Root2 * myD1u + myD1v;
439 else if (Abs(C) > RealEpsilon())
441 math_DirectPolynomialRoots Root(C, B, A);
443 if((Root.NbSolutions() != 2))
445 myCurvatureStatus = LProp_Undefined;
446 return Standard_False;
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);
456 VectCurv1 = myD1u + Root1 * myD1v;
457 VectCurv2 = myD1u + Root2 * myD1v;
472 myDirMinCurv = gp_Dir (VectCurv1);
473 myDirMaxCurv = gp_Dir (VectCurv2);
479 myDirMinCurv = gp_Dir (VectCurv2);
480 myDirMaxCurv = gp_Dir (VectCurv1);
483 myMeanCurv = ((N * E) - (2. * M * F) + (L * G)) // voir Farin p.282
484 / (2. * ((E * G) - (F * F)));
485 myGausCurv = ((L * N) - (M * M))
486 / ((E * G) - (F * F));
487 myCurvatureStatus = LProp_Computed;
488 return Standard_True;
491 Standard_Boolean LProp_SLProps::IsUmbilic ()
493 if(!IsCurvatureDefined())
494 throw LProp_NotDefined();
496 return Abs(myMaxCurv - myMinCurv) < Abs(Epsilon(myMaxCurv));
499 Standard_Real LProp_SLProps::MaxCurvature ()
501 if(!IsCurvatureDefined())
502 throw LProp_NotDefined();
507 Standard_Real LProp_SLProps::MinCurvature ()
509 if(!IsCurvatureDefined())
510 throw LProp_NotDefined();
515 void LProp_SLProps::CurvatureDirections(gp_Dir& Max, gp_Dir& Min)
517 if(!IsCurvatureDefined())
518 throw LProp_NotDefined();
524 Standard_Real LProp_SLProps::MeanCurvature ()
526 if(!IsCurvatureDefined())
527 throw LProp_NotDefined();
532 Standard_Real LProp_SLProps::GaussianCurvature ()
534 if(!IsCurvatureDefined())
535 throw LProp_NotDefined();