src/LProp/LProp_CLProps.gxx

   1 // Copyright (c) 1995-1999 Matra Datavision
   2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   3 //
   4 // This file is part of Open CASCADE Technology software library.
   5 //
   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.
  11 //
  12 // Alternatively, this file may be used under the terms of Open CASCADE
  13 // commercial license or contractual agreement.
  14
  15 #include <LProp_Status.hxx>
  16 #include <LProp_NotDefined.hxx>
  17 #include <Standard_OutOfRange.hxx>
  18
  19 static const Standard_Real MinStep   = 1.0e-7;
  20
  21
  22
  23 LProp_CLProps::LProp_CLProps (const Curve& C,
  24                               const Standard_Real U,
  25                               const Standard_Integer N,
  26                               const Standard_Real Resolution)
  27       : myCurve(C), myDerOrder(N), myCN(4),
  28         myLinTol(Resolution), myTangentStatus (LProp_Undecided)
  29 {
  30   Standard_OutOfRange_Raise_if (N < 0 || N > 3,
  31                           "LProp_CLProps::LProp_CLProps()");
  32
  33   SetParameter(U);
  34 }
  35
  36 LProp_CLProps::LProp_CLProps (const Curve& C, const Standard_Integer N,
  37                                     const Standard_Real Resolution)
  38         : myCurve(C), myU(RealLast()), myDerOrder(N), myCN(4),
  39             myLinTol(Resolution), myTangentStatus (LProp_Undecided)
  40 {
  41   Standard_OutOfRange_Raise_if (N < 0 || N > 3,
  42                           "LProp_CLProps::LProp_CLProps()");
  43 }
  44
  45 LProp_CLProps::LProp_CLProps (const Standard_Integer N,
  46                               const Standard_Real Resolution)
  47         : myU(RealLast()), myDerOrder(N), myCN(0), myLinTol(Resolution),
  48                             myTangentStatus (LProp_Undecided)
  49 {
  50   Standard_OutOfRange_Raise_if (N < 0 || N > 3, "");
  51 }
  52
  53 void LProp_CLProps::SetParameter(const Standard_Real U)
  54 {
  55   myU = U;
  56   switch (myDerOrder)
  57   {
  58   case 0:
  59     Tool::Value(myCurve, myU, myPnt);
  60     break;
  61   case 1:
  62     Tool::D1(myCurve, myU, myPnt, myDerivArr[0]);
  63     break;
  64   case 2:
  65     Tool::D2(myCurve, myU, myPnt, myDerivArr[0], myDerivArr[1]);
  66     break;
  67   case 3:
  68     Tool::D3(myCurve, myU, myPnt, myDerivArr[0], myDerivArr[1], myDerivArr[2]);
  69     break;
  70   }
  71
  72   myTangentStatus = LProp_Undecided;
  73 }
  74
  75 void LProp_CLProps::SetCurve(const Curve& C)
  76 {
  77   myCurve = C ;
  78   myCN = 4; // Tool::Continuity(C); RLE
  79 }
  80
  81 const Pnt& LProp_CLProps::Value () const
  82 {
  83   return myPnt;
  84 }
  85
  86 const Vec& LProp_CLProps::D1 ()
  87 {
  88   if (myDerOrder < 1)
  89   {
  90     myDerOrder = 1;
  91     Tool::D1(myCurve, myU, myPnt, myDerivArr[0]);
  92   }
  93
  94   return myDerivArr[0];
  95 }
  96
  97 const Vec& LProp_CLProps::D2 ()
  98 {
  99   if (myDerOrder < 2)
 100   {
 101     myDerOrder = 2;
 102     Tool::D2(myCurve, myU, myPnt, myDerivArr[0], myDerivArr[1]);
 103   }
 104
 105   return myDerivArr[1];
 106 }
 107
 108 const Vec& LProp_CLProps::D3 ()
 109 {
 110   if (myDerOrder < 3)
 111   {
 112     myDerOrder = 3;
 113     Tool::D3(myCurve, myU, myPnt, myDerivArr[0], myDerivArr[1], myDerivArr[2]);
 114   }
 115
 116   return myDerivArr[2];
 117 }
 118
 119 Standard_Boolean LProp_CLProps::IsTangentDefined ()
 120 {
 121   if (myTangentStatus == LProp_Undefined)
 122     return Standard_False;
 123   else if (myTangentStatus >= LProp_Defined)
 124     return Standard_True;
 125
 126   // tangentStatus == Lprop_Undecided
 127   // we have to calculate the first non null derivative
 128   const Standard_Real Tol = myLinTol * myLinTol;
 129
 130   Vec V;
 131
 132   Standard_Integer Order = 0;
 133   while (Order++ < 4)
 134   {
 135     if(myCN >= Order)
 136     {
 137       switch(Order)
 138       {
 139       case 1:
 140         V = D1();
 141         break;
 142       case 2:
 143         V = D2();
 144         break;
 145       case 3:
 146         V = D3();
 147         break;
 148       }//switch(Order)
 149
 150       if(V.SquareMagnitude() > Tol)
 151       {
 152         mySignificantFirstDerivativeOrder = Order;
 153         myTangentStatus = LProp_Defined;
 154         return Standard_True;
 155       }//if(V.SquareMagnitude() > Tol)
 156     }//if(cn >= Order)
 157     else
 158     {
 159       myTangentStatus = LProp_Undefined;
 160       return Standard_False;
 161     }// else of "if(cn >= Order)" condition
 162   }//while (Order < 4)
 163
 164   return Standard_False;
 165 }
 166
 167 void  LProp_CLProps::Tangent (Dir& D)
 168 {
 169   if(!IsTangentDefined())
 170     LProp_NotDefined::Raise();
 171
 172   if(mySignificantFirstDerivativeOrder == 1)
 173     D = Dir(myDerivArr[0]);
 174   else if (mySignificantFirstDerivativeOrder > 1)
 175   {
 176     const Standard_Real DivisionFactor = 1.e-3;
 177     const Standard_Real anUsupremum = Tool::LastParameter(myCurve),
 178                         anUinfium = Tool::FirstParameter(myCurve);
 179
 180     Standard_Real du;
 181     if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst()))
 182       du = 0.0;
 183     else
 184       du = anUsupremum-anUinfium;
 185
 186     const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
 187
 188     Vec V = myDerivArr[mySignificantFirstDerivativeOrder - 1];
 189
 190     Standard_Real u;
 191
 192     if(myU-anUinfium < aDelta)
 193       u = myU+aDelta;
 194     else
 195       u = myU-aDelta;
 196
 197     Pnt P1, P2;
 198     Tool::Value(myCurve, Min(myU, u),P1);
 199     Tool::Value(myCurve, Max(myU, u),P2);
 200
 201     Vec V1(P1,P2);
 202     Standard_Real aDirFactor = V.Dot(V1);
 203
 204     if(aDirFactor < 0.0)
 205       V = -V;
 206
 207     D = Dir(V);
 208   }//else if (mySignificantFirstDerivativeOrder > 1)
 209 }
 210
 211 Standard_Real LProp_CLProps::Curvature ()
 212 {
 213   Standard_Boolean isDefined = IsTangentDefined();
 214   (void)isDefined; // trick to avoid compiler warning on variable unised in Release mode; note that IsTangentDefined() must be called always
 215   LProp_NotDefined_Raise_if(!isDefined,
 216           "LProp_CLProps::CurvatureNotDefined()");
 217
 218   // if the first derivative is null the curvature is infinite.
 219   if(mySignificantFirstDerivativeOrder > 1)
 220     return RealLast();
 221
 222   Standard_Real Tol = myLinTol * myLinTol;
 223   Standard_Real DD1 = myDerivArr[0].SquareMagnitude();
 224   Standard_Real DD2 = myDerivArr[1].SquareMagnitude();
 225
 226   // if the second derivative is null the curvature is null.
 227   if (DD2 <= Tol)
 228   {
 229     myCurvature = 0.0;
 230   }
 231   else
 232   {
 233     Standard_Real N = myDerivArr[0].CrossSquareMagnitude(myDerivArr[1]);
 234     // if d[0] and d[1] are colinear the curvature is null.
 235     Standard_Real t = N/(DD1*DD2);
 236     if (t<=Tol)
 237     {
 238       myCurvature = 0.0;
 239     }
 240     else
 241     {
 242       myCurvature = sqrt(N) / (DD1*sqrt(DD1));
 243     }
 244   }
 245
 246   return myCurvature;
 247 }
 248
 249 void  LProp_CLProps::Normal (Dir& D)
 250 {
 251   Standard_Real c = Curvature();
 252   if(c==RealLast() || Abs(c) <= myLinTol)
 253   {
 254     LProp_NotDefined::Raise("LProp_CLProps::Normal(...):"
 255         "Curvature is null or infinity");
 256   }
 257
 258   // we used here the following vector relation
 259   // a ^ (b ^ c) = b(ac) - c(ab)
 260   // Norm = d[0] ^ (d[1] ^ d[0])
 261
 262   Vec Norm = myDerivArr[1] * (myDerivArr[0] * myDerivArr[0]) - myDerivArr[0] * (myDerivArr[0] * myDerivArr[1]);
 263   D = Dir(Norm);
 264 }
 265
 266 void  LProp_CLProps::CentreOfCurvature (Pnt& P)
 267 {
 268   if(Abs(Curvature()) <= myLinTol)
 269   {
 270     LProp_NotDefined::Raise();
 271   }
 272
 273   // we used here the following vector relation
 274   // a ^ (b ^ c) = b(ac) - c(ab)
 275   // Norm = d[0] ^ (d[1] ^ d[0])
 276
 277   Vec Norm = myDerivArr[1] * (myDerivArr[0] * myDerivArr[0]) - myDerivArr[0] * (myDerivArr[0] * myDerivArr[1]);
 278   Norm.Normalize();
 279   Norm.Divide(myCurvature);
 280   P= myPnt.Translated(Norm);
 281 }