0023706: Cannot project point on curve

[occt.git] / src / LProp / LProp_CLProps.gxx

// Copyright (c) 1995-1999 Matra Datavision
// Copyright (c) 1999-2012 OPEN CASCADE SAS
//
// The content of this file is subject to the Open CASCADE Technology Public
// License Version 6.5 (the "License"). You may not use the content of this file
// except in compliance with the License. Please obtain a copy of the License
// at http://www.opencascade.org and read it completely before using this file.
//
// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
//
// The Original Code and all software distributed under the License is
// distributed on an "AS IS" basis, without warranty of any kind, and the
// Initial Developer hereby disclaims all such warranties, including without
// limitation, any warranties of merchantability, fitness for a particular
// purpose or non-infringement. Please see the License for the specific terms
// and conditions governing the rights and limitations under the License.

#include <LProp_Status.hxx>
#include <LProp_NotDefined.hxx>
#include <Standard_OutOfRange.hxx>

static const Standard_Real MinStep   = 1.0e-7;



LProp_CLProps::LProp_CLProps (const Curve& C,
                              const Standard_Real U,
                              const Standard_Integer N, 
                              const Standard_Real Resolution) 
      : myCurve(C), myDerOrder(N), myCN(4), 
        myLinTol(Resolution), myTangentStatus (LProp_Undecided)
{
  Standard_OutOfRange_Raise_if (N < 0 || N > 3,
                          "LProp_CLProps::LProp_CLProps()");

  SetParameter(U);
}

LProp_CLProps::LProp_CLProps (const Curve& C, const Standard_Integer N, 
                                    const Standard_Real Resolution) 
        : myCurve(C), myU(RealLast()), myDerOrder(N), myCN(4), 
            myLinTol(Resolution), myTangentStatus (LProp_Undecided)
{
  Standard_OutOfRange_Raise_if (N < 0 || N > 3, 
                          "LProp_CLProps::LProp_CLProps()");  
}

LProp_CLProps::LProp_CLProps (const Standard_Integer N, 
                              const Standard_Real Resolution) 
        : myU(RealLast()), myDerOrder(N), myCN(0), myLinTol(Resolution),
                            myTangentStatus (LProp_Undecided)
{
  Standard_OutOfRange_Raise_if (N < 0 || N > 3, "");  
}

void LProp_CLProps::SetParameter(const Standard_Real U)
{
  myU = U;
  switch (myDerOrder)
  {
  case 0:
    Tool::Value(myCurve, myU, myPnt);
    break;
  case 1:
    Tool::D1(myCurve, myU, myPnt, myDerivArr[0]);
    break;
  case 2:
    Tool::D2(myCurve, myU, myPnt, myDerivArr[0], myDerivArr[1]);
    break;
  case 3:
    Tool::D3(myCurve, myU, myPnt, myDerivArr[0], myDerivArr[1], myDerivArr[2]);
    break;
  }

  myTangentStatus = LProp_Undecided;
}

void LProp_CLProps::SetCurve(const Curve& C)
{
  myCurve = C ; 
  myCN = 4; // Tool::Continuity(C); RLE
}

const Pnt& LProp_CLProps::Value () const
{
  return myPnt;
}

const Vec& LProp_CLProps::D1 ()
{
  if (myDerOrder < 1)
  {
    myDerOrder = 1;
    Tool::D1(myCurve, myU, myPnt, myDerivArr[0]);
  }

  return myDerivArr[0];
}

const Vec& LProp_CLProps::D2 ()
{
  if (myDerOrder < 2)
  {
    myDerOrder = 2;
    Tool::D2(myCurve, myU, myPnt, myDerivArr[0], myDerivArr[1]);
  }

  return myDerivArr[1];
}

const Vec& LProp_CLProps::D3 ()
{
  if (myDerOrder < 3)
  {
    myDerOrder = 3;
    Tool::D3(myCurve, myU, myPnt, myDerivArr[0], myDerivArr[1], myDerivArr[2]);
  }

  return myDerivArr[2];
}

Standard_Boolean LProp_CLProps::IsTangentDefined ()
{
  if (myTangentStatus == LProp_Undefined)
    return Standard_False;
  else if (myTangentStatus >= LProp_Defined)
    return Standard_True;

  // tangentStatus == Lprop_Undecided 
  // we have to calculate the first non null derivative
  const Standard_Real Tol = myLinTol * myLinTol;
  
  Vec V;
  
  Standard_Integer Order = 0;
  while (Order++ < 4)
  {
    if(myCN >= Order)
    {
      switch(Order)
      {
      case 1:
        V = D1();
        break;
      case 2:
        V = D2();
        break;
      case 3:
        V = D3();
        break;
      }//switch(Order)

      if(V.SquareMagnitude() > Tol)
      {
        mySignificantFirstDerivativeOrder = Order;
        myTangentStatus = LProp_Defined;
        return Standard_True;
      }//if(V.SquareMagnitude() > Tol)
    }//if(cn >= Order)
    else
    {
      myTangentStatus = LProp_Undefined;
      return Standard_False;
    }// else of "if(cn >= Order)" condition
  }//while (Order < 4)

  return Standard_False;
}

void  LProp_CLProps::Tangent (Dir& D)
{
  if(!IsTangentDefined())
    LProp_NotDefined::Raise();
  
  if(mySignificantFirstDerivativeOrder == 1)
    D = Dir(myDerivArr[0]);
  else if (mySignificantFirstDerivativeOrder > 1)
  {
    const Standard_Real DivisionFactor = 1.e-3;
    const Standard_Real anUsupremum = Tool::LastParameter(myCurve), 
                        anUinfium = Tool::FirstParameter(myCurve);
                        
    Standard_Real du;
    if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst()))
      du = 0.0;
    else
      du = anUsupremum-anUinfium;
    
    const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);

    Vec V = myDerivArr[mySignificantFirstDerivativeOrder - 1];
    
    Standard_Real u;
          
    if(myU-anUinfium < aDelta)
      u = myU+aDelta;
    else
      u = myU-aDelta;
    
    Pnt P1, P2;
    Tool::Value(myCurve, Min(myU, u),P1);
    Tool::Value(myCurve, Max(myU, u),P2);
    
    Vec V1(P1,P2);
    Standard_Real aDirFactor = V.Dot(V1);
    
    if(aDirFactor < 0.0)
      V = -V;
      
    D = Dir(V);
  }//else if (mySignificantFirstDerivativeOrder > 1)
}

Standard_Real LProp_CLProps::Curvature ()
{
  Standard_Boolean isDefined = IsTangentDefined();
  LProp_NotDefined_Raise_if(!isDefined,
          "LProp_CLProps::CurvatureNotDefined()");

  // if the first derivative is null the curvature is infinite.
  if(mySignificantFirstDerivativeOrder > 1)
    return RealLast();

  Standard_Real Tol = myLinTol * myLinTol;
  Standard_Real DD1 = myDerivArr[0].SquareMagnitude();
  Standard_Real DD2 = myDerivArr[1].SquareMagnitude();
  
  // if the second derivative is null the curvature is null.
  if (DD2 <= Tol)
  {
    myCurvature = 0.0;
  }
  else
  {
    Standard_Real N = myDerivArr[0].CrossSquareMagnitude(myDerivArr[1]);
    // if d[0] and d[1] are colinear the curvature is null.
    Standard_Real t = N/(DD1*DD2);
    if (t<=Tol)
    {
      myCurvature = 0.0;
    }
    else
    {
      myCurvature = sqrt(N) / (DD1*sqrt(DD1));
    }
  }

  return myCurvature;
}

void  LProp_CLProps::Normal (Dir& D)
{
  Standard_Real c = Curvature();
  if(c==RealLast() || Abs(c) <= myLinTol)
  {
    LProp_NotDefined::Raise("LProp_CLProps::Normal(...):"
        "Curvature is null or infinity"); 
  }

  // we used here the following vector relation 
  // a ^ (b ^ c) = b(ac) - c(ab)
  // Norm = d[0] ^ (d[1] ^ d[0])
  
  Vec Norm = myDerivArr[1] * (myDerivArr[0] * myDerivArr[0]) - myDerivArr[0] * (myDerivArr[0] * myDerivArr[1]);
  D = Dir(Norm);
}

void  LProp_CLProps::CentreOfCurvature (Pnt& P)
{
  if(Abs(Curvature()) <= myLinTol)
  {
    LProp_NotDefined::Raise(); 
  }

  // we used here the following vector relation 
  // a ^ (b ^ c) = b(ac) - c(ab)
  // Norm = d[0] ^ (d[1] ^ d[0])

  Vec Norm = myDerivArr[1] * (myDerivArr[0] * myDerivArr[0]) - myDerivArr[0] * (myDerivArr[0] * myDerivArr[1]);
  Norm.Normalize();
  Norm.Divide(myCurvature);
  P= myPnt.Translated(Norm);
}