Warning, /include/opencascade/LProp_SLProps.gxx is written in an unsupported language. File is not indexed.
0001 // Copyright (c) 1995-1999 Matra Datavision
0002 // Copyright (c) 1999-2014 OPEN CASCADE SAS
0003 //
0004 // This file is part of Open CASCADE Technology software library.
0005 //
0006 // This library is free software; you can redistribute it and/or modify it under
0007 // the terms of the GNU Lesser General Public License version 2.1 as published
0008 // by the Free Software Foundation, with special exception defined in the file
0009 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
0010 // distribution for complete text of the license and disclaimer of any warranty.
0011 //
0012 // Alternatively, this file may be used under the terms of Open CASCADE
0013 // commercial license or contractual agreement.
0014
0015 #include <LProp_Status.hxx>
0016 #include <LProp_NotDefined.hxx>
0017 #include <Standard_OutOfRange.hxx>
0018 #include <Standard_DomainError.hxx>
0019 #include <CSLib.hxx>
0020 #include <CSLib_DerivativeStatus.hxx>
0021 #include <CSLib_NormalStatus.hxx>
0022 #include <TColgp_Array2OfVec.hxx>
0023 #include <math_DirectPolynomialRoots.hxx>
0024
0025 static const Standard_Real MinStep = 1.0e-7;
0026
0027 static Standard_Boolean IsTangentDefined (LProp_SLProps& SProp,
0028 const Standard_Integer cn,
0029 const Standard_Real linTol,
0030 const Standard_Integer Derivative,
0031 Standard_Integer& Order,
0032 LProp_Status& theStatus)
0033 {
0034 Standard_Real Tol = linTol * linTol;
0035 gp_Vec V[2];
0036 Order = 0;
0037
0038 while (Order < 3)
0039 {
0040 Order++;
0041 if(cn >= Order)
0042 {
0043 switch(Order)
0044 {
0045 case 1:
0046 V[0] = SProp.D1U();
0047 V[1] = SProp.D1V();
0048 break;
0049 case 2:
0050 V[0] = SProp.D2U();
0051 V[1] = SProp.D2V();
0052 break;
0053 }//switch(Order)
0054
0055 if(V[Derivative].SquareMagnitude() > Tol)
0056 {
0057 theStatus = LProp_Defined;
0058 return Standard_True;
0059 }
0060 }//if(cn >= Order)
0061 else
0062 {
0063 theStatus = LProp_Undefined;
0064 return Standard_False;
0065 }
0066 }
0067
0068 return Standard_False;
0069 }
0070
0071 LProp_SLProps::LProp_SLProps (const Surface& S,
0072 const Standard_Real U,
0073 const Standard_Real V,
0074 const Standard_Integer N,
0075 const Standard_Real Resolution)
0076 : mySurf(S),myDerOrder(N), myCN(4), // (Tool::Continuity(S)),
0077 myLinTol(Resolution)
0078 {
0079 Standard_OutOfRange_Raise_if(N < 0 || N > 2,
0080 "LProp_SLProps::LProp_SLProps()");
0081
0082 SetParameters(U, V);
0083 }
0084
0085 LProp_SLProps::LProp_SLProps (const Surface& S,
0086 const Standard_Integer N,
0087 const Standard_Real Resolution)
0088 : mySurf(S), myU(RealLast()), myV(RealLast()), myDerOrder(N),
0089 myCN(4), // (Tool::Continuity(S))
0090 myLinTol(Resolution),
0091 myUTangentStatus (LProp_Undecided),
0092 myVTangentStatus (LProp_Undecided),
0093 myNormalStatus (LProp_Undecided),
0094 myCurvatureStatus(LProp_Undecided)
0095 {
0096 Standard_OutOfRange_Raise_if(N < 0 || N > 2,
0097 "LProp_SLProps::LProp_SLProps()");
0098 }
0099
0100 LProp_SLProps::LProp_SLProps (const Standard_Integer N,
0101 const Standard_Real Resolution)
0102 : myU(RealLast()), myV(RealLast()), myDerOrder(N), myCN(0),
0103 myLinTol(Resolution),
0104 myUTangentStatus (LProp_Undecided),
0105 myVTangentStatus (LProp_Undecided),
0106 myNormalStatus (LProp_Undecided),
0107 myCurvatureStatus(LProp_Undecided)
0108 {
0109 Standard_OutOfRange_Raise_if(N < 0 || N > 2,
0110 "LProp_SLProps::LProp_SLProps() bad level");
0111 }
0112
0113 void LProp_SLProps::SetSurface (const Surface& S )
0114 {
0115 mySurf = S;
0116 myCN = 4; // =Tool::Continuity(S);
0117 }
0118
0119 void LProp_SLProps::SetParameters (const Standard_Real U,
0120 const Standard_Real V)
0121 {
0122 myU = U;
0123 myV = V;
0124 switch (myDerOrder)
0125 {
0126 case 0:
0127 Tool::Value(mySurf, myU, myV, myPnt);
0128 break;
0129 case 1:
0130 Tool::D1(mySurf, myU, myV, myPnt, myD1u, myD1v);
0131 break;
0132 case 2:
0133 Tool::D2(mySurf, myU, myV, myPnt, myD1u, myD1v, myD2u, myD2v, myDuv);
0134 break;
0135 }
0136
0137 myUTangentStatus = LProp_Undecided;
0138 myVTangentStatus = LProp_Undecided;
0139 myNormalStatus = LProp_Undecided;
0140 myCurvatureStatus = LProp_Undecided;
0141 }
0142
0143 const gp_Pnt& LProp_SLProps::Value() const
0144 {
0145 return myPnt;
0146 }
0147
0148 const gp_Vec& LProp_SLProps::D1U()
0149 {
0150 if (myDerOrder < 1)
0151 {
0152 myDerOrder =1;
0153 Tool::D1(mySurf,myU,myV,myPnt,myD1u,myD1v);
0154 }
0155
0156 return myD1u;
0157 }
0158
0159 const gp_Vec& LProp_SLProps::D1V()
0160 {
0161 if (myDerOrder < 1)
0162 {
0163 myDerOrder =1;
0164 Tool::D1(mySurf,myU,myV,myPnt,myD1u,myD1v);
0165 }
0166
0167 return myD1v;
0168 }
0169
0170 const gp_Vec& LProp_SLProps::D2U()
0171 {
0172 if (myDerOrder < 2)
0173 {
0174 myDerOrder =2;
0175 Tool::D2(mySurf,myU,myV,myPnt,myD1u,myD1v,myD2u,myD2v,myDuv);
0176 }
0177
0178 return myD2u;
0179 }
0180
0181 const gp_Vec& LProp_SLProps::D2V()
0182 {
0183 if (myDerOrder < 2)
0184 {
0185 myDerOrder =2;
0186 Tool::D2(mySurf,myU,myV,myPnt,myD1u,myD1v,myD2u,myD2v,myDuv);
0187 }
0188
0189 return myD2v;
0190 }
0191
0192 const gp_Vec& LProp_SLProps::DUV()
0193 {
0194 if (myDerOrder < 2)
0195 {
0196 myDerOrder =2;
0197 Tool::D2(mySurf,myU,myV,myPnt,myD1u,myD1v,myD2u,myD2v,myDuv);
0198 }
0199
0200 return myDuv;
0201 }
0202
0203 Standard_Boolean LProp_SLProps::IsTangentUDefined ()
0204 {
0205 if (myUTangentStatus == LProp_Undefined)
0206 return Standard_False;
0207 else if (myUTangentStatus >= LProp_Defined)
0208 return Standard_True;
0209
0210 // uTangentStatus == Lprop_Undecided
0211 // we have to calculate the first non null U derivative
0212 return IsTangentDefined(*this, myCN, myLinTol, 0,
0213 mySignificantFirstDerivativeOrderU, myUTangentStatus);
0214 }
0215
0216 void LProp_SLProps::TangentU (gp_Dir& D)
0217 {
0218 if(!IsTangentUDefined())
0219 throw LProp_NotDefined();
0220
0221 if(mySignificantFirstDerivativeOrderU == 1)
0222 D = gp_Dir(myD1u);
0223 else
0224 {
0225 const Standard_Real DivisionFactor = 1.e-3;
0226 Standard_Real anUsupremum, anUinfium;
0227 Standard_Real anVsupremum, anVinfium;
0228 Tool::Bounds(mySurf,anUinfium,anVinfium,anUsupremum,anVsupremum);
0229 Standard_Real du;
0230 if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst()))
0231 du = 0.0;
0232 else
0233 du = anUsupremum-anUinfium;
0234
0235 const Standard_Real aDeltaU = Max(du*DivisionFactor,MinStep);
0236
0237 gp_Vec V = myD2u;
0238
0239 Standard_Real u;
0240
0241 if(myU-anUinfium < aDeltaU)
0242 u = myU+aDeltaU;
0243 else
0244 u = myU-aDeltaU;
0245
0246 gp_Pnt P1, P2;
0247 Tool::Value(mySurf, Min(myU, u),myV,P1);
0248 Tool::Value(mySurf, Max(myU, u),myV,P2);
0249
0250 gp_Vec V1(P1,P2);
0251 Standard_Real aDirFactor = V.Dot(V1);
0252
0253 if(aDirFactor < 0.0)
0254 V = -V;
0255
0256 D = gp_Dir(V);
0257 }
0258 }
0259
0260 Standard_Boolean LProp_SLProps::IsTangentVDefined ()
0261 {
0262 if (myVTangentStatus == LProp_Undefined)
0263 return Standard_False;
0264 else if (myVTangentStatus >= LProp_Defined)
0265 return Standard_True;
0266
0267 // vTangentStatus == Lprop_Undecided
0268 // we have to calculate the first non null V derivative
0269 return IsTangentDefined(*this, myCN, myLinTol, 1,
0270 mySignificantFirstDerivativeOrderV, myVTangentStatus);
0271 }
0272
0273 void LProp_SLProps::TangentV (gp_Dir& D)
0274 {
0275 if(!IsTangentVDefined())
0276 throw LProp_NotDefined();
0277
0278 if(mySignificantFirstDerivativeOrderV == 1)
0279 D = gp_Dir(myD1v);
0280 else
0281 {
0282 const Standard_Real DivisionFactor = 1.e-3;
0283 Standard_Real anUsupremum, anUinfium;
0284 Standard_Real anVsupremum, anVinfium;
0285 Tool::Bounds(mySurf,anUinfium,anVinfium,anUsupremum,anVsupremum);
0286
0287 Standard_Real dv;
0288 if((anVsupremum >= RealLast()) || (anVinfium <= RealFirst()))
0289 dv = 0.0;
0290 else
0291 dv = anVsupremum-anVinfium;
0292
0293 const Standard_Real aDeltaV = Max(dv*DivisionFactor,MinStep);
0294
0295 gp_Vec V = myD2v;
0296
0297 Standard_Real v;
0298
0299 if(myV-anVinfium < aDeltaV)
0300 v = myV+aDeltaV;
0301 else
0302 v = myV-aDeltaV;
0303
0304 gp_Pnt P1, P2;
0305 Tool::Value(mySurf, myU, Min(myV, v),P1);
0306 Tool::Value(mySurf, myU, Max(myV, v),P2);
0307
0308 gp_Vec V1(P1,P2);
0309 Standard_Real aDirFactor = V.Dot(V1);
0310
0311 if(aDirFactor < 0.0)
0312 V = -V;
0313
0314 D = gp_Dir(V);
0315 }
0316 }
0317
0318 Standard_Boolean LProp_SLProps::IsNormalDefined()
0319 {
0320 if (myNormalStatus == LProp_Undefined)
0321 return Standard_False;
0322 else if (myNormalStatus >= LProp_Defined)
0323 return Standard_True;
0324
0325 // status = UnDecided
0326 // first try the standard computation of the normal.
0327 CSLib_DerivativeStatus aStatus = CSLib_Done;
0328 CSLib::Normal(myD1u, myD1v, myLinTol, aStatus, myNormal);
0329 if (aStatus == CSLib_Done)
0330 {
0331 myNormalStatus = LProp_Computed;
0332 return Standard_True;
0333 }
0334
0335 // else solve the degenerated case only if continuity >= 2
0336
0337 myNormalStatus = LProp_Undefined;
0338 return Standard_False;
0339 }
0340
0341 const gp_Dir& LProp_SLProps::Normal ()
0342 {
0343 if(!IsNormalDefined())
0344 {
0345 throw LProp_NotDefined();
0346 }
0347
0348 return myNormal;
0349 }
0350
0351
0352 Standard_Boolean LProp_SLProps::IsCurvatureDefined ()
0353 {
0354 if (myCurvatureStatus == LProp_Undefined)
0355 return Standard_False;
0356 else if (myCurvatureStatus >= LProp_Defined)
0357 return Standard_True;
0358
0359 if(myCN < 2)
0360 {
0361 myCurvatureStatus = LProp_Undefined;
0362 return Standard_False;
0363 }
0364
0365 // status = UnDecided
0366 if (!IsNormalDefined())
0367 {
0368 myCurvatureStatus = LProp_Undefined;
0369 return Standard_False;
0370 }
0371
0372 // pour eviter un plantage dans le cas du caro pointu
0373 // en fait on doit pouvoir calculer les courbure
0374 // avoir
0375 if(!IsTangentUDefined() || !IsTangentVDefined())
0376 {
0377 myCurvatureStatus = LProp_Undefined;
0378 return Standard_False;
0379 }
0380
0381 // here we compute the curvature features of the surface
0382
0383 gp_Vec Norm (myNormal);
0384
0385 Standard_Real E = myD1u.SquareMagnitude();
0386 Standard_Real F = myD1u.Dot(myD1v);
0387 Standard_Real G = myD1v.SquareMagnitude();
0388
0389 if(myDerOrder < 2)
0390 this->D2U();
0391
0392 Standard_Real L = Norm.Dot(myD2u);
0393 Standard_Real M = Norm.Dot(myDuv);
0394 Standard_Real N = Norm.Dot(myD2v);
0395
0396 Standard_Real A = E * M - F * L;
0397 Standard_Real B = E * N - G * L;
0398 Standard_Real C = F * N - G * M;
0399
0400 Standard_Real MaxABC = Max(Max(Abs(A),Abs(B)),Abs(C));
0401 if (MaxABC < RealEpsilon()) // ombilic
0402 {
0403 myMinCurv = N / G;
0404 myMaxCurv = myMinCurv;
0405 myDirMinCurv = gp_Dir (myD1u);
0406 myDirMaxCurv = gp_Dir (myD1u.Crossed(Norm));
0407 myMeanCurv = myMinCurv; // (Cmin + Cmax) / 2.
0408 myGausCurv = myMinCurv * myMinCurv; // (Cmin * Cmax)
0409 myCurvatureStatus = LProp_Computed;
0410 return Standard_True;
0411 }
0412
0413 A = A / MaxABC;
0414 B = B / MaxABC;
0415 C = C / MaxABC;
0416 Standard_Real Curv1, Curv2, Root1, Root2;
0417 gp_Vec VectCurv1, VectCurv2;
0418
0419 if (Abs(A) > RealEpsilon())
0420 {
0421 math_DirectPolynomialRoots Root (A, B, C);
0422 if(Root.NbSolutions() != 2)
0423 {
0424 myCurvatureStatus = LProp_Undefined;
0425 return Standard_False;
0426 }
0427 else
0428 {
0429 Root1 = Root.Value(1);
0430 Root2 = Root.Value(2);
0431 Curv1 = ((L * Root1 + 2. * M) * Root1 + N) /
0432 ((E * Root1 + 2. * F) * Root1 + G);
0433 Curv2 = ((L * Root2 + 2. * M) * Root2 + N) /
0434 ((E * Root2 + 2. * F) * Root2 + G);
0435 VectCurv1 = Root1 * myD1u + myD1v;
0436 VectCurv2 = Root2 * myD1u + myD1v;
0437 }
0438 }
0439 else if (Abs(C) > RealEpsilon())
0440 {
0441 math_DirectPolynomialRoots Root(C, B, A);
0442
0443 if((Root.NbSolutions() != 2))
0444 {
0445 myCurvatureStatus = LProp_Undefined;
0446 return Standard_False;
0447 }
0448 else
0449 {
0450 Root1 = Root.Value(1);
0451 Root2 = Root.Value(2);
0452 Curv1 = ((N * Root1 + 2. * M) * Root1 + L) /
0453 ((G * Root1 + 2. * F) * Root1 + E);
0454 Curv2 = ((N * Root2 + 2. * M) * Root2 + L) /
0455 ((G * Root2 + 2. * F) * Root2 + E);
0456 VectCurv1 = myD1u + Root1 * myD1v;
0457 VectCurv2 = myD1u + Root2 * myD1v;
0458 }
0459 }
0460 else
0461 {
0462 Curv1 = L / E;
0463 Curv2 = N / G;
0464 VectCurv1 = myD1u;
0465 VectCurv2 = myD1v;
0466 }
0467
0468 if (Curv1 < Curv2)
0469 {
0470 myMinCurv = Curv1;
0471 myMaxCurv = Curv2;
0472 myDirMinCurv = gp_Dir (VectCurv1);
0473 myDirMaxCurv = gp_Dir (VectCurv2);
0474 }
0475 else
0476 {
0477 myMinCurv = Curv2;
0478 myMaxCurv = Curv1;
0479 myDirMinCurv = gp_Dir (VectCurv2);
0480 myDirMaxCurv = gp_Dir (VectCurv1);
0481 }
0482
0483 myMeanCurv = ((N * E) - (2. * M * F) + (L * G)) // voir Farin p.282
0484 / (2. * ((E * G) - (F * F)));
0485 myGausCurv = ((L * N) - (M * M))
0486 / ((E * G) - (F * F));
0487 myCurvatureStatus = LProp_Computed;
0488 return Standard_True;
0489 }
0490
0491 Standard_Boolean LProp_SLProps::IsUmbilic ()
0492 {
0493 if(!IsCurvatureDefined())
0494 throw LProp_NotDefined();
0495
0496 return Abs(myMaxCurv - myMinCurv) < Abs(Epsilon(myMaxCurv));
0497 }
0498
0499 Standard_Real LProp_SLProps::MaxCurvature ()
0500 {
0501 if(!IsCurvatureDefined())
0502 throw LProp_NotDefined();
0503
0504 return myMaxCurv;
0505 }
0506
0507 Standard_Real LProp_SLProps::MinCurvature ()
0508 {
0509 if(!IsCurvatureDefined())
0510 throw LProp_NotDefined();
0511
0512 return myMinCurv;
0513 }
0514
0515 void LProp_SLProps::CurvatureDirections(gp_Dir& Max, gp_Dir& Min)
0516 {
0517 if(!IsCurvatureDefined())
0518 throw LProp_NotDefined();
0519
0520 Max = myDirMaxCurv;
0521 Min = myDirMinCurv;
0522 }
0523
0524 Standard_Real LProp_SLProps::MeanCurvature ()
0525 {
0526 if(!IsCurvatureDefined())
0527 throw LProp_NotDefined();
0528
0529 return myMeanCurv;
0530 }
0531
0532 Standard_Real LProp_SLProps::GaussianCurvature ()
0533 {
0534 if(!IsCurvatureDefined())
0535 throw LProp_NotDefined();
0536
0537 return myGausCurv;
0538 }
