ImplicitPolynomial3Shape.ih

#include <cstdlib>

// IMPLEMENTATION of inline methods.

// ----------------------- Standard services ------------------------------

//-----------------------------------------------------------------------------
template <typename TSpace>
inline
DGtal::ImplicitPolynomial3Shape<TSpace>::~ImplicitPolynomial3Shape()
{
}
//-----------------------------------------------------------------------------
template <typename TSpace>
inline
DGtal::ImplicitPolynomial3Shape<TSpace>::
ImplicitPolynomial3Shape( const Polynomial3 & poly )
{
  init( poly );
}
//-----------------------------------------------------------------------------
template <typename TSpace>
inline
DGtal::ImplicitPolynomial3Shape<TSpace> &
DGtal::ImplicitPolynomial3Shape<TSpace>::
operator=( const ImplicitPolynomial3Shape & other )
{
  if ( this != &other )
  {
    myPolynomial = other.myPolynomial;

    myFx= other.myFx;
    myFy= other.myFy;
    myFz= other.myFz;

    myFxx= other.myFxx;
    myFxy= other.myFxy;
    myFxz= other.myFxz;

    myFyx= other.myFyx;
    myFyy= other.myFyy;
    myFyz= other.myFyz;

    myFzx= other.myFzx;
    myFzy= other.myFzy;
    myFzz= other.myFzz;
  }
  return *this;
}
//-----------------------------------------------------------------------------
template <typename TSpace>
inline
void
DGtal::ImplicitPolynomial3Shape<TSpace>::
init( const Polynomial3 & poly )
{
  myPolynomial = poly;

  myFx= derivative<0>( poly );
  myFy= derivative<1>( poly );
  myFz= derivative<2>( poly );

  myFxx= derivative<0>( myFx );
  myFxy= derivative<1>( myFx );
  myFxz= derivative<2>( myFx);

  myFyx= derivative<0>( myFy );
  myFyy= derivative<1>( myFy );
  myFyz= derivative<2>( myFy );

  myFzx= derivative<0>( myFz );
  myFzy= derivative<1>( myFz );
  myFzz= derivative<2>( myFz );

  myUpPolynome = myFx*(myFx*myFxx+myFy*myFyx+myFz*myFzx)+
                                myFy*(myFx*myFxy+myFy*myFyy+myFz*myFzy)+
                                myFz*(myFx*myFxz+myFy*myFyz+myFz*myFzz)-
                                ( myFx*myFx +myFy*myFy+myFz*myFz )*(myFxx+myFyy+myFzz);

  myLowPolynome = myFx*myFx +myFy*myFy+myFz*myFz;
}
//-----------------------------------------------------------------------------
template <typename TSpace>
inline
double
DGtal::ImplicitPolynomial3Shape<TSpace>::
operator()(const RealPoint &aPoint) const
{
  return myPolynomial( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] );
}
//-----------------------------------------------------------------------------
template <typename TSpace>
inline
bool
DGtal::ImplicitPolynomial3Shape<TSpace>::
isInside(const RealPoint &aPoint) const
{
  return (this->operator()(aPoint) > (Ring)0);
}
//-----------------------------------------------------------------------------
template <typename TSpace>
inline
DGtal::Orientation
DGtal::ImplicitPolynomial3Shape<TSpace>::
orientation(const RealPoint &aPoint) const
{
  Ring v = this->operator()(aPoint);
  if ( v > (Ring)0 )
    return INSIDE;
  else if ( v < (Ring)0 )
    return OUTSIDE;
  else
    return ON;
}
//-----------------------------------------------------------------------------
template <typename TSpace>
inline
typename DGtal::ImplicitPolynomial3Shape<TSpace>::RealVector
DGtal::ImplicitPolynomial3Shape<TSpace>::
gradient( const RealPoint &aPoint ) const
{
  // ISO C++ tells that an object created at return time will not be
  // copied into the caller context, but will be already defined in
  // the correct context.
  return RealVector
      ( myFx ( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] ),
        myFy ( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] ),
        myFz ( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] ) );

}


// ------------------------------------------------------------ Added by Anis Benyoub
//-----------------------------------------------------------------------------

template <typename TSpace>
inline
double
DGtal::ImplicitPolynomial3Shape<TSpace>::
meanCurvature( const RealPoint &aPoint ) const
{
  double temp= myLowPolynome( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] );
  temp = sqrt(temp);
  double downValue = 2*(temp*temp*temp);
  double upValue = myUpPolynome( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] );


  return -(upValue/downValue);
}



//-----------------------------------------------------------------------------
template <typename TSpace>
inline
double
DGtal::ImplicitPolynomial3Shape<TSpace>::
gaussianCurvature( const RealPoint &aPoint ) const
{

  double vFx= myFx( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] );
  double vFy= myFy( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] );
  double vFz= myFz( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] );

  double vFxx= myFxx( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] );
  double vFxy= myFxy( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] );
  double vFxz= myFxz( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] );

  //double vFyx= myFyx( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] );
  double vFyy= myFyy( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] );
  double vFyz= myFyz( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] );

  
  /*double vFzx = myFzx( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] );
  double vFzy = myFzy( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] );
  */
  double vFzz = myFzz( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] );
 

  double A = vFz*(vFxx*vFz-2.0*vFx*vFxz)+vFx*vFx*vFzz;

  double B = vFz*(vFyy*vFz-2.0*vFy*vFyz)+vFy*vFy*vFzz;

  double C = vFz*(-vFx*vFyz+vFxy*vFz-vFxz*vFy)+vFx*vFy*vFzz;

  double D = vFz*(vFx*vFx+vFy*vFy+vFz*vFz);

  double G= (A*B-C*C)/(D*D);

  return G;

}


template <typename TSpace>
inline
typename DGtal::ImplicitPolynomial3Shape<TSpace>::RealPoint 
DGtal::ImplicitPolynomial3Shape<TSpace>::nearestPoint( const RealPoint &aPoint, const double accuracy, 
                                                       const int maxIter, const double gamma ) const
{
   RealPoint agradient= (*this).gradient(aPoint);
   RealPoint X=aPoint;
   int numberIter=0;
   while((fabs((*this)(X))>=accuracy) && (numberIter<maxIter))
   {
        double norm=agradient.norm();
        RealPoint normalizedGradient= RealPoint(agradient[0]/norm,agradient[1]/norm,agradient[2]/norm);
        double alpha =gamma;
        if((*this)(X)>0)
        {
                alpha=-alpha;
        }
        else
        {

        }
        
        X=X+normalizedGradient*alpha;
        agradient=  (*this).gradient(X);
        numberIter++;
   }
        return X;
}

// Interface - public :

template <typename TSpace>
inline
void
DGtal::ImplicitPolynomial3Shape<TSpace>::selfDisplay ( std::ostream & out ) const
{
  out << "[ImplicitPolynomial3Shape] P(x,y,z) = " << myPolynomial;
}

template <typename TSpace>
inline
bool
DGtal::ImplicitPolynomial3Shape<TSpace>::isValid() const
{
  return true;
}



// Implementation of inline functions                                        //

template <typename TSpace>
inline
std::ostream&
DGtal::operator<< ( std::ostream & out,
                    const ImplicitPolynomial3Shape<TSpace> & object )
{
  object.selfDisplay( out );
  return out;
}

//                                                                           //