// 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;


class npoint3D
{
public:
  virtual ~npoint3D() {}
  virtual double x() const=0;
  virtual double y() const=0;
  virtual double z() const=0;
  virtual double w() const=0;
  virtual double operator()(const int i) const=0; // always able to return 4 items (3D x y z + weight)
};

class npoint : public npoint3D
{
  double coord[4];
public :
  virtual ~npoint() {}
  npoint()
  {
    coord[0]=coord[1]=coord[2]=0.0;
    coord[3]=1.0;
  }
  npoint(const npoint3 pt, const double w=1) ;
  npoint(const npoint2 pt, const double w=1) ;
  npoint(const double pt, const double w) ;
  npoint(double x,double y,double z)
  {
    coord[0]=x;
    coord[1]=y;
    coord[2]=z;
    coord[3]=1.;
  }
  npoint(double x,double y,double z,double w)
  {
    coord[0]=x;
    coord[1]=y;
    coord[2]=z;
    coord[3]=w;
  }
  double & operator[](const int i)
  {
    return coord[i];  // renvoie la ieme coordonnee (modifiable)
  }
  double operator[](const int i) const
  {
    return coord[i];  // renvoie la ieme coordonnee (constante)
  }
  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];
  }

  double & wx()
  {
    return coord[0];  // renvoie w*x
  }
  double wx() const
  {
    return coord[0];
  }
  virtual double x() const
  {
    if (coord[3]!=0.0)
      return coord[0]/coord[3];
    else
      return coord[0];
  }
  double & wy()
  {
    return coord[1];  // "" w*y
  }
  double wy() const
  {
    return coord[1];
  }
  virtual double y() const
  {
    if (coord[3]!=0.0)
      return coord[1]/coord[3];
    else
      return coord[1];
  }
  double & wz()
  {
    return coord[2];  // "" w*z
  }
  double wz() const
  {
    return coord[2];
  }
  virtual double z() const
  {
    if (coord[3]!=0.0)
      return coord[2]/coord[3];
    else
      return coord[2];
  }
  double & w()
  {
    return coord[3];  // "" w
  }
  virtual double w() const
  {
    return coord[3];
  }
  double * array()
  {
    return coord;
  };                                     // "" adresse de la premiere coordonee
  const double * array() const
  {
    return coord;  // "" adresse de la premiere coordonee
  }

  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;
  }
  npoint operator+ (npoint other) const                            // "" addition ""
  {
    for (int i=0;i<4;++i) other.coord[i]+=coord[i];
    return other;
  }
  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;
  }
  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;
  }
  npoint& operator+= (const npoint other)                                // ajout "en place"
  {
    for (int i=0;i<4;++i) coord[i]+=other.coord[i];
    return *this;
  }
  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
  double norm() const
  {
    return sqrt(coord[0]*coord[0]+coord[1]*coord[1]+coord[2]*coord[2]+coord[3]*coord[3]);
  }
  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]);
  }
  
  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;
  }
  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;
class npoint3 : public npoint3D
{
  double coord[3];
public :
  virtual ~npoint3() {}
  npoint3()
  {
    coord[0]=coord[1]=coord[2]=0.;
  }
  npoint3(const double *xyz)
  {
    coord[0]=xyz[0];
    coord[1]=xyz[1];  // initialisation
    coord[2]=xyz[2];
  }
  npoint3(double x,double y,double z)
  {
    coord[0]=x;
    coord[1]=y;  // initialisation
    coord[2]=z;
  }
  npoint3(npoint p)
  {
    p.perspective_divide();
    coord[0]=p[0];
    coord[1]=p[1];
    coord[2]=p[2];
  }
  npoint3(const npoint2 p);
  npoint3(const double pt)
  {
    coord[0]=pt;
    coord[1]=0.;
    coord[2]=0.;
  }

  double & operator[](const int i)
  {
    return coord[i];  // renvoie la ieme coordonnee (modifiable)
  }
  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;
  }

  double & x()
  {
    return coord[0];  // renvoie x
  }
  virtual double x() const
  {
    return coord[0];
  }
  double & y()
  {
    return coord[1];  // "" y
  }
  virtual double y() const
  {
    return coord[1];
  }
  double & z()
  {
    return coord[2];  // "" z
  }
  virtual double z() const
  {
    return coord[2];
  }
  virtual double w() const
  {
    return 1.;
  }

  double * array()
  {
    return coord;  // "" adresse de la premiere coordonee
  }
  const double * array() const
  {
    return coord;  // "" adresse de la premiere coordonee
  }
  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;
  }
  double dotprod(const npoint3  other) const
  {
    return coord[0]*other.coord[0]+coord[1]*other.coord[1]+coord[2]*other.coord[2];
  }
  double norm2() const
  {
    return dotprod(*this);
  }

  double norm() const
  {
    return sqrt(dotprod(*this));
  }
  double normalize()
  {
    double n=norm();
    coord[0]/=n;
    coord[1]/=n;
    coord[2]/=n;
    return n;
  }
  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;
  }
  npoint3 operator+ (npoint3 other) const                            // "" addition ""
  {
    for (int i=0;i<3;++i) other[i]+=coord[i];
    return other;
  }
  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;
  }
  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;
  }
  npoint3& operator+= (const npoint3 other)                                // ajout "en place"
  {
    for (int i=0;i<3;++i) coord[i]+=other.coord[i];
    return *this;
  }
  npoint3& operator-= (const npoint3 other)                                // soustraction "en place"
  {
    for (int i=0;i<3;++i) coord[i]-=other.coord[i];
    return *this;
  }
  npoint3& operator*= (double fact)                                // soustraction "en place"
  {
    for (int i=0;i<3;++i) coord[i]*=fact;
    return *this;
  }
  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
  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);

class npoint2 : public npoint3D
{
  double coord[2];
public :
  virtual ~npoint2() {}
  npoint2()
  {
    coord[0]=coord[1]=0.;
  }
  npoint2(double pt)
  {
    coord[0]=pt;coord[1]=0.;  // initialisation
  }
  npoint2(const npoint3 pt)
  {
    coord[0]=pt[0];coord[1]=pt[1];  // discard z
  }
  npoint2(npoint pt)
  {
    pt.perspective_divide();
    coord[0]=pt[0];coord[1]=pt[1];  // discard z
  }

  npoint2(double u,double v)
  {
    coord[0]=u;coord[1]=v;  // initialisation
  }
  double & operator[](const int i)
  {
    return coord[i];  // renvoie la ieme coordonnee (modifiable)
  }
  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;
  }

  double & u()
  {
    return coord[0];
  }
  double u() const
  {
    return coord[0];
  }
  double & v()
  {
    return coord[1];
  }
  double v() const
  {
    return coord[1];
  }
  virtual double x() const
  {
    return coord[0];
  }
  virtual double y() const
  {
    return coord[1];
  }
  virtual double z() const
  {
    return 0.;
  }
  virtual double w() const
  {
    return 1.;
  }

  double * array()
  {
    return coord;
  };

  npoint2 operator- (npoint2 other) const
  {
    for (int i=0;i<2;++i) other.coord[i]=coord[i]-other.coord[i];
    return other;
  }
  npoint2 operator+ (npoint2 other) const
  {
    for (int i=0;i<2;++i) other.coord[i]+=coord[i];
    return other;
  }
  npoint2 operator/ (const double fact) const
  {
    npoint2 buf=*this;
    for (int i=0;i<2;++i) buf[i]/=fact;
    return buf;
  }
  npoint2 operator*(const double fact) const
  {
    npoint2 buf=*this;
    for (int i=0;i<2;++i) buf[i]*=fact;
    return buf;
  }
  npoint2& operator+= (const npoint2 other)
  {
    for (int i=0;i<2;++i) coord[i]+=other.coord[i];
    return *this;
  }
  npoint2& operator-= (const npoint2 other)
  {
    for (int i=0;i<2;++i) coord[i]-=other.coord[i];
    return *this;
  }
  npoint2& operator*= (double fact)
  {
    for (int i=0;i<2;++i) coord[i]*=fact;
    return *this;
  }
  npoint2& operator/= (double fact)
  {
    for (int i=0;i<2;++i) coord[i]/=fact;
    return *this;
  }
  double dotprod(const npoint2  other) const
  {
    return coord[0]*other.coord[0]+coord[1]*other.coord[1];
  }
  double norm2() const
  {
    return dotprod(*this);
  }
  double norm() const
  {
    return sqrt(dotprod(*this));
  }
  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