// nUtil - An utility Library for gnurbs
// Copyright (C) 2008-2026 Eric Bechet
//
// See the LICENSE file for contributions and license information.
// Please report all bugs and problems to <bechet@cadxfem.org>.
//

#ifndef __NPOINT_H
#define __NPOINT_H

#include <cmath>
#include <iostream>

#ifndef n_pi
#define n_pi (3.1415926535897932384626433832795)
#endif

class npoint2;
class npoint3;

/// 3D point (generic class).
class npoint3D
{
public:
  /// \brief Destructor.
  virtual ~npoint3D() {}

  /// \brief Returns x coordinate (constant version).
  /// \return Coordinate.
  virtual double x() const=0;

  /// \brief Returns y coordinate (constant version).
  /// \return Coordinate.
  virtual double y() const=0;

  /// \brief Returns z coordinate (constant version).
  /// \return Coordinate.
  virtual double z() const=0;

  /// \brief Returns homogeneous coordinate (constant version).
  /// \return Coordinate.
  virtual double w() const=0;

  /// \brief Returns the ith coordinate (from 0 to 3).
  /// \return Coordinate.
  virtual double operator()(const int i) const=0; // always able to return 4 items (3D x y z + weight)
};

/// 4D point.
class npoint : public npoint3D
{
  double coord[4]; //!< Homogeneous coordinates.
public :

  /// \brief Destructor.
  virtual ~npoint() {}

  /// \brief Default constructor.
  npoint()
  {
    coord[0]=coord[1]=coord[2]=0.0;
    coord[3]=1.0;
  }

  /// \brief Constructor a npoint from a npoint3.
  /// \param[in] pt 3D point.
  /// \param[in] w weight.
  npoint(const npoint3 pt, const double w=1) ;

  /// \brief Constructor a npoint from a npoint2.
  /// \param[in] pt 2D point.
  /// \param[in] w weight.
  npoint(const npoint2 pt, const double w=1) ;

  /// \brief Constructor a npoint from a value.
  /// \param[in] pt 1D point.
  /// \param[in] w weight.
  npoint(const double pt, const double w) ;

  /// \brief Constructor a npoint from three scalars.
  /// \param[in] x x-coordinate.
  /// \param[in] y y-coordinate.
  /// \param[in] z z-coordinate.
  npoint(double x,double y,double z)
  {
    coord[0]=x;
    coord[1]=y;
    coord[2]=z;
    coord[3]=1.;
  }

  /// \brief Constructor a npoint from four scalars.
  /// \param[in] x x-coordinate.
  /// \param[in] y y-coordinate.
  /// \param[in] z z-coordinate.
  /// \param[in] w weight.
  npoint(double x,double y,double z,double w)
  {
    coord[0]=x;
    coord[1]=y;
    coord[2]=z;
    coord[3]=w;
  }

  /// \brief Returns the ith homogeneous coordinate (from 0 to 3) - mutable version.
  /// \param[in] i wanted coordinate number.
  /// \return Coordinate.
  double & operator[](const int i)
  {
    return coord[i];  // renvoie la ieme coordonnee (modifiable)
  }

  /// \brief Returns the ith homogeneous coordinate (from 0 to 3) - constant version.
  /// \param[in] i wanted coordinate number.
  /// \return Coordinate.
  double operator[](const int i) const
  {
    return coord[i];  // renvoie la ieme coordonnee (constante)
  }

  /// \brief Returns the ith coordinate (from 0 to 3).
  /// \warning If \f$ i \in [0, 2] \f$ a perspective division is performed on the returned coordinate.
  /// \param[in] i wanted coordinate number.
  /// \return Coordinate.
  virtual double operator()(const int i) const
  {
    if (i<3)
      if (coord[3]!=0.0) return coord[i]/coord[3]; else return coord[i];
    else return coord[3];
  }

  /// \brief Returns the homogeneous x coordinate (mutable version).
  /// \return Coordinate.
  double & wx()
  {
    return coord[0];  // renvoie w*x
  }

  /// \brief Returns the homogeneous x coordinate (constant version).
  /// \return Coordinate.
  double wx() const
  {
    return coord[0];
  }

  /// \brief Returns the x coordinate.
  /// \warning A perspective division is performed.
  /// \return Coordinate.
  virtual double x() const
  {
    if (coord[3]!=0.0)
      return coord[0]/coord[3];
    else
      return coord[0];
  }

  /// \brief Returns the homogeneous y coordinate (mutable version).
  /// \return Coordinate.
  double & wy()
  {
    return coord[1];  // "" w*y
  }

  /// \brief Returns the homogeneous y coordinate (constant version).
  /// \return Coordinate.
  double wy() const
  {
    return coord[1];
  }

  /// \brief Returns the y coordinate.
  /// \warning A perspective division is performed.
  /// \return Coordinate.
  virtual double y() const
  {
    if (coord[3]!=0.0)
      return coord[1]/coord[3];
    else
      return coord[1];
  }

  /// \brief Returns the homogeneous z coordinate (mutable version).
  /// \return Coordinate.
  double & wz()
  {
    return coord[2];  // "" w*z
  }

  /// \brief Returns the homogeneous z coordinate (constant version).
  /// \return Coordinate.
  double wz() const
  {
    return coord[2];
  }

  /// \brief Returns the z coordinate.
  /// \warning A perspective division is performed.
  /// \return Coordinate.
  virtual double z() const
  {
    if (coord[3]!=0.0)
      return coord[2]/coord[3];
    else
      return coord[2];
  }

  /// \brief Returns the weight (mutable version).
  /// \return Coordinate.
  double & w()
  {
    return coord[3];  // "" w
  }

  /// \brief Returns the weight (constant version).
  /// \return Coordinate.
  virtual double w() const
  {
    return coord[3];
  }

  /// \brief Returns array on data (mutable version).
  /// \return Array.
  double * array()
  {
    return coord;
  };                                     // "" adresse de la premiere coordonee

  /// \brief Returns array on data (constant version).
  /// \return Array.
  const double * array() const
  {
    return coord;  // "" adresse de la premiere coordonee
  }

  /// \brief Substract this point \f$ P \f$ with another point: \f$ P_{out} = P - P_{other} \f$.
  /// \param[in] other \f$ P_{other} \f$.
  /// \return \f$ P_{out} \f$.
  npoint operator- (npoint other) const                            // renvoie la soustraction des coords. de deux points
  {
    for (int i=0;i<4;++i) other.coord[i]=coord[i]-other.coord[i];
    return other;
  }

  /// \brief Adds this point \f$ P \f$ with another point: \f$ P_{out} = P + P_{other} \f$.
  /// \param[in] other \f$ P_{other} \f$.
  /// \return \f$ P_{out} \f$.
  npoint operator+ (npoint other) const                            // "" addition ""
  {
    for (int i=0;i<4;++i) other.coord[i]+=coord[i];
    return other;
  }

  /// \brief Divides this point \f$ P \f$ with a scalar: \f$ P_{out} = P / a \f$.
  /// \param[in] fact \f$ a \f$.
  /// \return \f$ P_{out} \f$.
  npoint operator/ (const double fact) const                              // "" division par une constante
  {
    npoint buf=*this;
    for (int i=0;i<4;++i) buf[i]/=fact;
    return buf;
  }

  /// \brief Right-multiply this point \f$ P \f$ with a scalar: \f$ P_{out} = P\ a \f$.
  /// \param[in] fact \f$ a \f$.
  /// \return \f$ P_{out} \f$.
  npoint operator*(const double fact) const                               // "" multiplication (a droite) par un scalaire
  {
    npoint buf=*this;
    for (int i=0;i<4;++i) buf[i]*=fact;
    return buf;
  }

  /// \brief Adds this point with another, this point is modified.
  /// param[in] other other point.
  /// \return Modified point.
  npoint& operator+= (const npoint other)                                // ajout "en place"
  {
    for (int i=0;i<4;++i) coord[i]+=other.coord[i];
    return *this;
  }

  /// \brief Substract this point with another, this point is modified.
  /// param[in] other other point.
  /// \return Modified point.
  npoint& operator-= (const npoint other)                                // soustraction "en place"
  {
    for (int i=0;i<4;++i) coord[i]-=other.coord[i];
    return *this;
  }

  
  friend npoint operator * (const double fact,const npoint other);      // "" multiplication (a gauche) par un scalaire

  /// Euclidean norm in 4D/homogeneous space.
  /// \return Norm.
  double norm() const
  {
    return sqrt(coord[0]*coord[0]+coord[1]*coord[1]+coord[2]*coord[2]+coord[3]*coord[3]);
  }

  /// Euclidean norm in 3D space.
  /// \return Norm.
  double norm3D() const
  {
    if (w()!=1.0)
    {
      if (w()!=0.0)
        return sqrt((coord[0]*coord[0]+coord[1]*coord[1]+coord[2]*coord[2])/(coord[3]*coord[3]));
      else
        return sqrt(coord[0]*coord[0]+coord[1]*coord[1]+coord[2]*coord[2]);
    }
    else return sqrt(coord[0]*coord[0]+coord[1]*coord[1]+coord[2]*coord[2]);
  }
  
  /// \brief Compute the perspective division of this point.
  /// \return True if point in 3D space, false if vector.
  bool perspective_divide() // returns true if point in 3D space, false if vector
  {
    double w=coord[3];
    if (w!=1.0)
    {
      if (w!=0.0)
      {
        for (int i=0;i<4;++i) coord[i]/=w;
        return true;
      }
      else return false;
    }
    else return true;
  }

  /// \brief Print this point to output stream.
  /// \param[in,out] os output stream.
  /// \param[in] prec printing precision.
  void print(std::ostream &os,int prec=5) const ;                        // sort les coordonees dans le flux os
};

npoint operator * (const double fact,const npoint other);
class npoint2;

/// 3D point.
class npoint3 : public npoint3D
{
  double coord[3]; //!< Array of coordinates.
public :

  /// \brief Destructor.
  virtual ~npoint3() {}

  /// \brief Default constructor.
  npoint3()
  {
    coord[0]=coord[1]=coord[2]=0.;
  }

  /// \brief Initializes point from 3D array.
  /// \param[in] xyz 3D array.
  npoint3(const double *xyz)
  {
    coord[0]=xyz[0];
    coord[1]=xyz[1];  // initialisation
    coord[2]=xyz[2];
  }

  /// \brief Initializes point from 3 scalars.
  /// \param x first coordinate.
  /// \param y second coordinate.
  /// \param z third coordinate.
  npoint3(double x,double y,double z)
  {
    coord[0]=x;
    coord[1]=y;  // initialisation
    coord[2]=z;
  }

  /// \brief Initializes 3D point from a 4D point.
  /// \warning A perspective division is performed.
  /// \param[in] p 4D point.
  npoint3(npoint p)
  {
    p.perspective_divide();
    coord[0]=p[0];
    coord[1]=p[1];
    coord[2]=p[2];
  }

  npoint3(const npoint2 p);

  /// \brief Initialization from a scalar.
  /// Coordinate x = pt. Other coordinates are set to zero.
  /// \param[in] pt scalar.
  npoint3(const double pt)
  {
    coord[0]=pt;
    coord[1]=0.;
    coord[2]=0.;
  }

  /// \brief Gets the ith coordinate (mutable version).
  /// \param[in] i coordinate number.
  /// \return Coordinate.
  double & operator[](const int i)
  {
    return coord[i];  // renvoie la ieme coordonnee (modifiable)
  }

  /// \brief Gets the ith coordinate (constant version).
  /// \param[in] i coordinate number.
  /// \return Coordinate.
  double operator[](const int i) const
  {
    return coord[i];  // renvoie la ieme coordonnee (constante)
  }

  virtual double operator()(const int i) const
  {
    if (i<3)
      return coord[i];
    else return 1.0;
  }

  /// \brief Returns the x coordinate (mutable version).
  /// \return Coordinate.
  double & x()
  {
    return coord[0];  // renvoie x
  }
  virtual double x() const
  {
    return coord[0];
  }

  /// \brief Returns the y coordinate (mutable version).
  /// \return Coordinate.
  double & y()
  {
    return coord[1];  // "" y
  }
  virtual double y() const
  {
    return coord[1];
  }

  /// \brief Returns the z coordinate (mutable version).
  /// \return Coordinate.
  double & z()
  {
    return coord[2];  // "" z
  }
  virtual double z() const
  {
    return coord[2];
  }
  virtual double w() const
  {
    return 1.;
  }

  /// \brief Returns the data array of coordinates (mutable version).
  /// \return Data array.
  double * array()
  {
    return coord;  // "" adresse de la premiere coordonee
  }

  /// \brief Returns the data array of coordinates (constant version).
  /// \return Data array.
  const double * array() const
  {
    return coord;  // "" adresse de la premiere coordonee
  }

  /// Computes the cross product between two points \f$ this = v_1 \wedge v_2 \f$
  /// and stores the result into this point.
  /// \param[in] V1 left point.
  /// \param[in] V2 right point.
  /// \return Cross product.
  npoint3& crossprod(const npoint3  V1,const npoint3 V2)
  {
    coord[0]=V1[1]* V2[2]-V1[2]* V2[1];
    coord[1]=V1[2]* V2[0]-V1[0]* V2[2];
    coord[2]=V1[0]* V2[1]-V1[1]* V2[0];
    return *this;
  }

  /// \brief Computes the dot product with another point.
  /// \param[in] other other point.
  /// \return Dot product.
  double dotprod(const npoint3  other) const
  {
    return coord[0]*other.coord[0]+coord[1]*other.coord[1]+coord[2]*other.coord[2];
  }

  /// \brief Computes the squared euclidean norm.
  /// \return Squared norm.
  double norm2() const
  {
    return dotprod(*this);
  }

  /// \brief Computes the euclidean norm.
  /// \return Norm.
  double norm() const
  {
    return sqrt(dotprod(*this));
  }

  /// \brief Normalizes this point (w.r.t the euclidean norm).
  /// \return Euclidean norm.
  double normalize()
  {
    double n=norm();
    coord[0]/=n;
    coord[1]/=n;
    coord[2]/=n;
    return n;
  }

  /// \brief Substracts this point with another. This point is not modified.
  /// \param[in] other other point.
  /// \return Result.
  npoint3 operator- (npoint3 other) const                            // renvoie la soustraction des coords. de deux points
  {
    for (int i=0;i<3;++i) other[i]=coord[i]-other.coord[i];
    return other;
  }

  /// \brief Adds this point with another. This point is not modified.
  /// \param[in] other other point.
  /// \return Result.
  npoint3 operator+ (npoint3 other) const                            // "" addition ""
  {
    for (int i=0;i<3;++i) other[i]+=coord[i];
    return other;
  }

  /// \brief Divides this point by a factor. This point is not modified.
  /// \param[in] fact factor.
  /// \return Result.
  npoint3 operator/ (const double fact) const                              // "" division par une constante
  {
    npoint3 buf=*this;
    for (int i=0;i<3;++i) buf[i]/=fact;
    return buf;
  }

  /// \brief Multiplies this point by a factor. This point is not modified.
  /// \param[in] fact factor.
  /// \return Result.
  npoint3 operator*(const double fact) const                               // "" multiplication (a droite) par un scalaire
  {
    npoint3 buf=*this;
    for (int i=0;i<3;++i) buf[i]*=fact;
    return buf;
  }

  /// \brief Adds this point with another. This point is modified.
  /// \param[in] other other point.
  /// \return Result.
  npoint3& operator+= (const npoint3 other)                                // ajout "en place"
  {
    for (int i=0;i<3;++i) coord[i]+=other.coord[i];
    return *this;
  }

  /// \brief Substracts this point with another. This point is modified.
  /// \param[in] other other point.
  /// \return Result.
  npoint3& operator-= (const npoint3 other)                                // soustraction "en place"
  {
    for (int i=0;i<3;++i) coord[i]-=other.coord[i];
    return *this;
  }

  /// \brief Multiplies this point by a factor. This point is modified.
  /// \param[in] fact factor.
  /// \return Result.
  npoint3& operator*= (double fact)                                // soustraction "en place"
  {
    for (int i=0;i<3;++i) coord[i]*=fact;
    return *this;
  }

  /// \brief Divides this point by a factor. This point is modified.
  /// \param[in] fact factor.
  /// \return Result.
  npoint3& operator/= (double fact)                                // soustraction "en place"
  {
    for (int i=0;i<3;++i) coord[i]/=fact;
    return *this;
  }
  friend npoint3 operator * (const double fact,const npoint3 other);      // "" multiplication (a gauche) par un scalaire

  /// \brief Print this point to output stream.
  /// \param[in,out] os output stream.
  /// \param[in] prec printing precision.
  void print(std::ostream &os,int prec=5) const;                         // sort les coordonees dans le flux os
};


npoint3 operator * (const double fact,const npoint3 other);

double randr(void);

double operator * (const npoint3 p1, const npoint3 p2);

npoint3 crossprod(const npoint3  V1,const npoint3 V2);

/// 2D point.
class npoint2 : public npoint3D
{
  double coord[2]; //!< Array of coordinates.
public :

  /// \brief Destructor.
  virtual ~npoint2() {}

  /// \brief default constructor.
  npoint2()
  {
    coord[0]=coord[1]=0.;
  }

  /// Initializes the first coordinate by a scalar.
  /// The second coordinate is set to zero.
  /// \param[in] pt scalar.
  npoint2(double pt)
  {
    coord[0]=pt;coord[1]=0.;  // initialisation
  }

  /// \brief Initializes this point from a 3D point.
  /// The z coordinate is discarded.
  /// \param[in] pt 3D point.
  npoint2(const npoint3 pt)
  {
    coord[0]=pt[0];coord[1]=pt[1];  // discard z
  }

  /// \brief Initializes this point from a 3D point.
  /// The z coordinate is discarded.
  /// \warning A perspective division is performed.
  /// \param[in] pt 4D point.
  npoint2(npoint pt)
  {
    pt.perspective_divide();
    coord[0]=pt[0];coord[1]=pt[1];  // discard z
  }

  /// \brief Initialization from two scalars.
  /// \param[in] u x coordinate.
  /// \param[in] v y coordinate.
  npoint2(double u,double v)
  {
    coord[0]=u;coord[1]=v;  // initialisation
  }

  /// \brief Gets the ith coordinate (mutable version).
  /// \param[in] i coordinate number.
  /// \return Coordinate.
  double & operator[](const int i)
  {
    return coord[i];  // renvoie la ieme coordonnee (modifiable)
  }

  /// \brief Gets the ith coordinate (constant version).
  /// \param[in] i coordinate number.
  /// \return Coordinate.
  double operator[](const int i) const
  {
    return coord[i];  // renvoie la ieme coordonnee (constante)
  }

  virtual double operator()(const int i) const
  {
    if (i<2)
      return coord[i];
    else return 1.0;
  }

  /// \brief Returns the u (x) coordinate (mutable version).
  /// \return Coordinate.
  double & u()
  {
    return coord[0];
  }

  /// \brief Returns the u (x) coordinate (constant version).
  /// \return Coordinate.
  double u() const
  {
    return coord[0];
  }

  /// \brief Returns the v (y) coordinate (mutable version).
  /// \return Coordinate.
  double & v()
  {
    return coord[1];
  }

  /// \brief Returns the v (y) coordinate (constant version).
  /// \return Coordinate.
  double v() const
  {
    return coord[1];
  }

  /// \brief Returns the x coordinate (constant version).
  /// \return Coordinate.
  virtual double x() const
  {
    return coord[0];
  }

  /// \brief Returns the y coordinate (constant version).
  /// \return Coordinate.
  virtual double y() const
  {
    return coord[1];
  }

  /// \brief Returns the z coordinate (constant version).
  /// \return Coordinate.
  virtual double z() const
  {
    return 0.;
  }

  /// \brief Returns the weigh (constant version).
  /// For 3D points, this is always 1.0.
  /// \return Coordinate.
  virtual double w() const
  {
    return 1.;
  }

  /// \brief Returns the data array of coordinates (mutable version).
  /// \return Data array.
  double * array()
  {
    return coord;
  };

  /// \brief Substracts this point with another. This point is not modified.
  /// \param[in] other other point.
  /// \return Result.
  npoint2 operator- (npoint2 other) const
  {
    for (int i=0;i<2;++i) other.coord[i]=coord[i]-other.coord[i];
    return other;
  }

  /// \brief Adds this point with another. This point is not modified.
  /// \param[in] other other point.
  /// \return Result.
  npoint2 operator+ (npoint2 other) const
  {
    for (int i=0;i<2;++i) other.coord[i]+=coord[i];
    return other;
  }

  /// \brief Divides this point by a factor. This point is not modified.
  /// \param[in] fact factor.
  /// \return Result.
  npoint2 operator/ (const double fact) const
  {
    npoint2 buf=*this;
    for (int i=0;i<2;++i) buf[i]/=fact;
    return buf;
  }

  /// \brief Multiplies this point by a factor. This point is not modified.
  /// \param[in] fact factor.
  /// \return Result.
  npoint2 operator*(const double fact) const
  {
    npoint2 buf=*this;
    for (int i=0;i<2;++i) buf[i]*=fact;
    return buf;
  }

  /// \brief Adds this point with another. This point is modified.
  /// \param[in] other other point.
  /// \return Result.
  npoint2& operator+= (const npoint2 other)
  {
    for (int i=0;i<2;++i) coord[i]+=other.coord[i];
    return *this;
  }

  /// \brief Substracts this point with another. This point is modified.
  /// \param[in] other other point.
  /// \return Result.
  npoint2& operator-= (const npoint2 other)
  {
    for (int i=0;i<2;++i) coord[i]-=other.coord[i];
    return *this;
  }

  /// \brief Multiplies this point by a factor. This point is modified.
  /// \param[in] fact factor.
  /// \return Result.
  npoint2& operator*= (double fact)
  {
    for (int i=0;i<2;++i) coord[i]*=fact;
    return *this;
  }

  /// \brief Divides this point by a factor. This point is modified.
  /// \param[in] fact factor.
  /// \return Result.
  npoint2& operator/= (double fact)
  {
    for (int i=0;i<2;++i) coord[i]/=fact;
    return *this;
  }

  /// \brief Computes the dot product with another point.
  /// \param[in] other other point.
  /// \return Dot product.
  double dotprod(const npoint2  other) const
  {
    return coord[0]*other.coord[0]+coord[1]*other.coord[1];
  }

  /// \brief Computes the squared euclidean norm.
  /// \return Squared norm.
  double norm2() const
  {
    return dotprod(*this);
  }

  /// \brief Computes the euclidean norm.
  /// \return Norm.
  double norm() const
  {
    return sqrt(dotprod(*this));
  }

  /// \brief Normalizes this point (w.r.t the euclidean norm).
  /// \return Euclidean norm.
  double normalize()
  {
    double n=norm();
    coord[0]/=n;
    coord[1]/=n;
    return n;
  }
};

npoint2 operator * (const double fact,const npoint2 other);
double operator * (const npoint2 p1, const npoint2 p2);
double crossprod(const npoint2  V1,const npoint2 V2);


#endif // __NPOINT_H