#ifndef POINT4D_HH 
#define POINT4D_HH
////////////////////////////////////////////////////////////////////////
//! file="amma/Geometry/AnGeo4/Point4d.hh"
//! lib=Mag4
//! userlevel=Normal
//! author="Radek Marik"
//! date="26.10.1993"
//! docentry="Geometry.4-D"
//! rcsid="$Id: Point4d.hh,v 1.8 2000/10/04 12:48:59 ees1cg Exp $"

#include "amma/StdMath.hh" //::Abs(),::Sqrt()
#include "amma/Boolean.hh"
#include "amma/Index.hh"
#include "amma/StdType.hh"

class Vector4dC;
class Matrix4d4C;
class Line4dPVC;
class PlanePVV4dC;
class istream;
class ostream;

// -------------------------------------------------------------------
// ******** Point4dC *************************************************
// -------------------------------------------------------------------


//: Point in 4D space
// The class Point4dC represents a geometrical entity point in 4 dimensional
// Euclidian space. A point is represented by a quadruple of coordinates.
// The index of a coordinate can be 0 for x, 1 for y, 2 for z, 3 for t.

class Point4dC
{
  public:

    // Type definitions
    // ----------------
    typedef StdTypeC::RealT CoordT;    // An coordinate.
    typedef StdTypeC::RealT DistanceT; // Distance of points.
    typedef StdTypeC::RealT AngleT;    // An angle in rad.
    
    typedef StdTypeC::RealT CosAngleT;
    // The cosinus of an angle of three points.

  public:

    // Constructors, assigment, copy, and destructor
    // ---------------------------------------------

    Point4dC(istream & inS);
    // Reads and creates the point from the input stream.

    inline Point4dC();
    // Creates the point (0.0, 0.0, 0.0, 0.0).

    inline Point4dC(CoordT xCoor, CoordT yCoor, CoordT zCoor, CoordT tCoor);
    // Creates the point (xCoor, yCoor, zCoor, tCoor).
 
    inline Point4dC(const Point4dC & point); 
    // Copy constructor.

    inline const Point4dC & operator=(const Point4dC & p);
    // Assigment.

    // Access to the object
    // --------------------

    inline IndexT N() const;
    // Returns the number of coordinates.

    inline const CoordT & operator[](IndexT i) const;
    // Access to the i-th constant coordinate.
 
    inline CoordT & operator[](IndexT i);
    // Access to the i-th coordinate.

    inline CoordT & X();
    // Access to the first coordinate.

    inline CoordT & Y();
    // Access to the second coordinate.
 
    inline CoordT & Z();
    // Access to the third coordinate.
 
    inline CoordT & T();
    // Access to the fourth coordinate.
 
    inline const CoordT & X() const; 
    // Access to the first constant coordinate.

    inline const CoordT & Y() const;
    // Access to the second constant coordinate.
 
    inline const CoordT & Z() const;
    // Access to the third constant coordinate.

    inline const CoordT & T() const;
    // Access to the fourth constant coordinate.

    inline const Point4dC & Point4d() const;
    // Access to the constant point.

    inline Point4dC & Point4d(); 
    // Access to the point.

    // Logical operators with points
    // -----------------------------

    inline BooleanT operator==(const Point4dC & point) const;
    // Returns TRUE if 2 points are the same point. Two points
    // are claimed to be same iff their relative distance error
    // is small (under 'relZero' threshold). If both points
    // are close to the origin of the coordinate system, their
    // absolute distance error is compared to 'zeroDistance'.  

    inline BooleanT operator!=(const Point4dC & point) const;
    // Returns TRUE if 2 points are different. 2 points are different
    // iff they are not equal in the sense of the operator '=='.

    inline BooleanT operator>(const Point4dC & point) const;
    // Returns TRUE if all coordinates are bigger then the coordinates
    // of the point.

    inline BooleanT OneBigger(const Point4dC & point) const;
    // Returns TRUE if one coordinate is bigger then the related one
    // of the point.

    // Setting of the point
    // --------------------

    inline Point4dC & SetZero(); 
    // Sets all coordinates to be 0.0.

    inline void Set(const Point4dC & point);
    // Sets coordinates according to the values of 'point'.

    inline void Set(const CoordT xCoor,
                      const CoordT yCoor,
                      const CoordT zCoor,
                      const CoordT tCoor);
    // Sets the point to the new coordinates.

    inline Point4dC & SetMinimum(const Point4dC & point);
    // Sets the point coordinates to minimal value from both points.

    inline Point4dC & SetMaximum(const Point4dC & point);
    // Sets the point coordinates to maximal value from both points.

    // Information about extrems
    // -------------------------

    inline IndexT MinIndex() const;
    // Returns the index of the minimum value.

    inline IndexT MaxIndex() const;
    // Returns the index of the maximum value.

    inline IndexT AbsMinIndex() const;
    // Returns the index of the minimum coordinate in absolute value.

    inline IndexT AbsMaxIndex() const;
    // Returns the index of the maximum coordinate in absolute value.

    inline CoordT Min() const;
    // Returns the minimum coordinate.

    inline CoordT Max() const;
    // Returns the maximum coordinate.
 
    // Arithmetical operations
    // -----------------------

    inline const Point4dC & operator+=(const Point4dC & point);
    // Adds the point.

    inline const Point4dC & operator-=(const Point4dC & point);
    // Subtracts the point.

    inline const Point4dC & operator*=(const Point4dC & point);
    // Multiplies items by point items.

    inline const Point4dC & operator/=(const Point4dC & point);
    // Divides items by point items.

    inline const Point4dC & operator*=(const RealT alpha);
    // Multiplies items by 'alpha'.

    inline const Point4dC & operator/=(const RealT alpha);
    // Divides items by 'alpha'.

    inline Point4dC operator+(const Point4dC & v) const;
    // Adds 2 points.

    Vector4dC operator-(const Point4dC & v) const;
    // Returns the difference between 2 points.

    inline Point4dC operator*(const Point4dC & v) const;
    // Returns the points which coordinates are results of multiplication
    // of this point and 'v' points item by item.

    inline Point4dC operator/(const Point4dC & v) const;
    // Returns the points which coordinates are results of division
    // of this point and 'v' points item by item.

    inline Point4dC operator*(const RealT alpha) const;
    // Returns the points which coordinates are these coordinates
    // multiplied by 'alpha'.

    inline Point4dC operator/(const RealT alpha) const;
    // Returns the points which coordinates are these coordinates
    // divided by 'alpha'.

    Point4dC operator+(const Vector4dC & v) const;
    // Returns the point = 'this point' + v;

    inline RealT operator&(const Point4dC & p) const;
    // Returns the sum X()*p.X() + Y()*p.Y() + Z()*p.Z(). 
    // In 3 dimensional Euclidian vector space the result is 
    // equivalent to dot product of 2 vectors.

    // Euclidian geometry
    // ------------------

    inline Point4dC Translation(const Point4dC & newOrigin) const;
    // Returns the point with coordinates related to the new origin
    // 'newOrigin'.

    inline RealT SqrEuclidDistance(const Point4dC & point) const;
    // Returns the square of Euclid distance between this point
    // and 'point'.

    inline RealT EuclidDistance(const Point4dC & point) const;
    // Returns the Euclid distance between this point and point.

    inline RealT EDistanceFromOrigin() const;
    // Returns the Euclid distance of this point from the coordinate origin.

    inline CosAngleT CosAngle(const Point4dC & m, const Point4dC & b) const;
    // Returns the cosinus of the angle (this point, m, b).

    inline AngleT Angle(const Point4dC & m, const Point4dC & b) const;
    // Returns the angle (this point, m, b). The angle range in rad is 0 to PI
    // and the angle is not oriented.

    AngleT AngleOrigin(const Point4dC & p) const;
    // Returns the angle in rad between points 'this point', (0,0,0), and
    // the point 'p'. The angle range is 0 to PI
    // and the angle is not oriented.

    inline Point4dC MidPoint(const Point4dC & p) const;
    // Returns the point which is in the middle of this point 
    // and the point 'p'.

    inline Point4dC Center(const Point4dC & point) const;
    // Returns the point which is in the middle of both points.

    // Useful functions with points
    // ----------------------------

    inline RealT AbsDistance(const Point4dC & point) const;
    // Returns the block distance between this point and 'point'.

    RealT MahalanobisDistance(const Point4dC   & mean, 
                              const Matrix4d4C & invCov);
    // Returns the Mahalanobis distance from the 'mean'. The 'invCov'
    // is the inversion of the covariance matrix.

    inline CoordT Sum() const;
    // Returns the sum of coordinates.

    inline CoordT SumAbs() const;
    // Returns the sum of absolute value of coordinates.

    inline Point4dC & Sqrt();
    // Changes all coordinates into their square root values.

    inline Point4dC & Abs();
    // Changes all coordinates into their absolute values.

    inline Point4dC & Normalize();
    // Scales the coordinates of the point to have the distance
    // from the origin equal to 1. The point is supposed to be different
    // from the origin.

    inline Point4dC & NormalizeSum();
    // Scales the coordinates of the point to have unit sum.
    // The point is supposed to be different from the origin.

    // Relationship with a line
    // ------------------------

    Point4dC OrthoProjectionInto(const Line4dPVC & l) const;
    // Returns the orthogonal projection of the point into the line 'l'.

    // Relationship with a plane
    // -------------------------

    Point4dC OrthoProjectionInto(const PlanePVV4dC & plane) const;
    // Returnss the orthogonal projection of the point into the plane.

    // Computational precision
    // -----------------------

    inline void FixAlmostZero();
    // If the absolute value  of a coordinate is smaller than zoroDistance
    // the function sets the coordinate to be zero.

    static RealT SetZeroDistance(RealT zero);
    // Sets the precision of test of zero distance of points.
    // The function returns the previous value.

    static inline BooleanT IsAlmostZero(RealT a);
    // Returns TRUE if 'a' is smaller or equal to zeroDistance.

    static RealT SetRelZero(RealT zero);
    // Sets the precision of the relative error test.
    // The function returns the previous value.

    static inline BooleanT IsRelZero(RealT a);
    // Returns TRUE if 'a' is smaller than or equal to 'relZero'.

  protected:

    // Special functions
    // -----------------

    static inline RealT RAbs(RealT t);
    // Returns the absolute value of the parameter t.

  private:
    enum {dim = 4};                // number of coordinates     
    static const RealT defRelZero; // default value of 'relZero'
    static RealT relZero;          // small relative error
    static const RealT defZeroDistance; // default value of 'zeroDistance'
    static RealT zeroDistance;     // distance of 2 points smaller than
                                   // 'zeroDistance' is treated as zero.
    static RealT sqrZeroDistance;  // square of 'zeroDistance'.

  private:
    CoordT  a[dim];

  friend inline Point4dC operator*(RealT alpha, const Point4dC & point);
  friend inline Point4dC operator/(RealT alpha, const Point4dC & point);
  friend ostream & operator<<(ostream & outS, const Point4dC & point);
  friend istream & operator>>(istream & inS, Point4dC & point);
};

inline 
Point4dC operator*(RealT alpha, const Point4dC & point);
// Returns the points which coordinates are 'point' coordinates multiplied
// by 'alpha'.

inline 
Point4dC operator/(RealT alpha, const Point4dC & point);
// Returns the points which coordinates are equal 'alpha' divided by 
// 'point' coordinates.


ostream & operator<<(ostream & outS, const Point4dC & point);
// Saves the 'point' into the output stream.

istream & operator>>(istream & inS, Point4dC & point);
// Set the coordinates of the 'point' according to the information
// stored in the input stream.

#include <iostream.h>

// -----------------------------------------------------------------
// ******** Point4dC ***********************************************
// -----------------------------------------------------------------

inline 
const Point4dC & 
Point4dC::operator=(const Point4dC & p)
{
  a[0] = p.a[0];
  a[1] = p.a[1];
  a[2] = p.a[2];
  a[3] = p.a[3];
  return *this;
}

inline 
RealT 
Point4dC::RAbs(RealT t)
{ return (t >= 0) ? t : -t; }

inline 
BooleanT 
Point4dC::IsAlmostZero(RealT a)
{ return RAbs(a) <= zeroDistance; }

inline 
BooleanT 
Point4dC::IsRelZero(RealT a)
{ return RAbs(a) <= relZero; }

inline 
IndexT 
Point4dC::N() const
{ return dim; }

inline 
Point4dC::CoordT &
Point4dC::X() 
{ return a[0]; }

inline 
Point4dC::CoordT & 
Point4dC::Y() 
{ return a[1];}

inline 
Point4dC::CoordT & 
Point4dC::Z()
{ return a[2]; }

inline 
Point4dC::CoordT & 
Point4dC::T()
{ return a[3]; }

inline 
const Point4dC::CoordT & 
Point4dC::X() const
{ return a[0]; }

inline 
const Point4dC::CoordT & 
Point4dC::Y() const
{ return a[1]; }

inline 
const Point4dC::CoordT & 
Point4dC::Z() const
{ return a[2]; }

inline 
const Point4dC::CoordT & 
Point4dC::T() const
{ return a[3]; }

inline 
Point4dC::Point4dC()
{
  X() = 0;
  Y() = 0;
  Z() = 0;
  T() = 0;
}

inline 
Point4dC::Point4dC(Point4dC::CoordT xCoor,
                   Point4dC::CoordT yCoor,
                   Point4dC::CoordT zCoor,
                   Point4dC::CoordT tCoor)
{
  X() = xCoor;
  Y() = yCoor;
  Z() = zCoor;
  T() = tCoor;
}

inline 
Point4dC::Point4dC(const Point4dC & point)
{
  X() = point.X();
  Y() = point.Y();
  Z() = point.Z();
  T() = point.T();
}

inline 
Point4dC::CoordT 
Point4dC::SumAbs() const
{ return RAbs(X()) + RAbs(Y()) + RAbs(Z()) + RAbs(T()); }

inline 
Point4dC 
Point4dC::Translation(const Point4dC & newOrigin) const
{
  return Point4dC( X()-newOrigin.X(), Y()-newOrigin.Y(),
                   Z()-newOrigin.Z(), T()-newOrigin.T() );
}

inline
BooleanT 
Point4dC::operator==(const Point4dC & p) const
{
  const CoordT diffSum = Translation(p).SumAbs(); 
  const CoordT sum = SumAbs() + p.SumAbs();
  return IsAlmostZero(sum) 
         ? IsAlmostZero(diffSum)
         : IsRelZero(diffSum/sum);
}

inline
BooleanT 
Point4dC::operator!=(const Point4dC & p) const
{ return ! operator==(p); }

inline 
BooleanT 
Point4dC::operator>(const Point4dC & point) const
{
  return    (X() > point.X())
         && (Y() > point.Y())
         && (Z() > point.Z())
         && (T() > point.T());
}

inline 
BooleanT 
Point4dC::OneBigger(const Point4dC & point) const
{
  return    (X() > point.X())
         || (Y() > point.Y())
         || (Z() > point.Z())
         || (T() > point.T());
}

inline 
const Point4dC::CoordT & 
Point4dC::operator[](IndexT i) const
{ return a[i.V()]; }

inline 
Point4dC::CoordT & 
Point4dC::operator[](IndexT i)
{ return a[i.V()]; }

inline 
const Point4dC &
Point4dC::Point4d() const
{ return *this; }

inline 
Point4dC &
Point4dC::Point4d()
{ return *this; }

inline 
IndexT 
Point4dC::MinIndex() const
{
  UIntT ind = 0;
  for (UIntT i = 1; i < dim; i++)
    if (a[ind] > a[i])
      ind = i;
  return ind;
}

inline 
IndexT 
Point4dC::MaxIndex() const
{
  UIntT ind = 0;
  for (UIntT i = 1; i < dim; i++)
    if (a[ind] < a[i])
      ind = i;
  return ind;
}

inline 
IndexT 
Point4dC::AbsMinIndex() const
{
  UIntT ind = 0;
  for (UIntT i = 1; i < dim; i++)
    if (RAbs(a[ind]) > RAbs(a[i]))
      ind = i;
  return ind;
}

inline 
IndexT 
Point4dC::AbsMaxIndex() const
{
  UIntT ind = 0;
  for (UIntT i = 1; i < dim; i++)
    if (RAbs(a[ind]) < RAbs(a[i]))
      ind = i;
  return ind;
}


inline 
Point4dC::CoordT 
Point4dC::Min() const
{
  CoordT ext = a[0];
  for (UIntT i = 1; i < dim; i++)
    if (a[i] < ext)
      ext = a[i];
  return ext;
}

inline 
Point4dC::CoordT 
Point4dC::Max() const
{
  CoordT ext = a[0];
  for (UIntT i = 1; i < dim; i++)
    if (a[i] > ext)
      ext = a[i];
  return ext;
}
 
inline 
Point4dC & 
Point4dC::SetZero()
{
  X() = 0;
  Y() = 0;
  Z() = 0;
  T() = 0;
  return *this;
}

inline 
void 
Point4dC::Set(const Point4dC::CoordT xCoor, const Point4dC::CoordT yCoor,
              const Point4dC::CoordT zCoor, const Point4dC::CoordT tCoor)
{
  X() = xCoor;
  Y() = yCoor;
  Z() = zCoor;
  T() = tCoor;
}

inline
void
Point4dC::Set(const Point4dC & point)
{
  X() = point.X();
  Y() = point.Y();
  Z() = point.Z();
  T() = point.T();
}

inline 
Point4dC & 
Point4dC::SetMinimum(const Point4dC & point)
{
  if (X() > point.X()) X() = point.X();
  if (Y() > point.Y()) Y() = point.Y();
  if (Z() > point.Z()) Z() = point.Z();
  if (T() > point.T()) T() = point.T();
  return *this;
}

inline 
Point4dC & 
Point4dC::SetMaximum(const Point4dC & point)
{
  if (X() < point.X()) X() = point.X();
  if (Y() < point.Y()) Y() = point.Y();
  if (Z() < point.Z()) Z() = point.Z();
  if (T() < point.T()) T() = point.T();
  return *this;
}


inline 
const Point4dC & 
Point4dC::operator+=(const Point4dC & point)
{
  X() += point.X();
  Y() += point.Y();
  Z() += point.Z();
  T() += point.T();
  return *this;
}

inline 
const Point4dC & 
Point4dC::operator-=(const Point4dC & point)
{
  X() -= point.X();
  Y() -= point.Y();
  Z() -= point.Z();
  T() -= point.T();
  return *this;
}

inline 
const Point4dC & 
Point4dC::operator*=(const Point4dC & point)
{
  X() *= point.X();
  Y() *= point.Y();
  Z() *= point.Z();
  T() *= point.T();
  return *this;
}

inline 
const Point4dC & 
Point4dC::operator/=(const Point4dC & point)
{
  X() /= point.X();
  Y() /= point.Y();
  Z() /= point.Z();
  T() /= point.T();
  return *this;
}

inline 
const Point4dC & 
Point4dC::operator*=(const RealT alpha)
{
  X() *= alpha;
  Y() *= alpha;
  Z() *= alpha;
  T() *= alpha;
  return *this;
}

inline 
const Point4dC & 
Point4dC::operator/=(const RealT alpha)
{
  X() /= alpha;
  Y() /= alpha;
  Z() /= alpha;
  T() /= alpha;
  return *this;
}

inline 
Point4dC 
Point4dC::operator+(const Point4dC & p) const
{ return Point4dC( X()+p.X(), Y()+p.Y(), Z()+p.Z(), T()+p.T() ); }

inline 
Point4dC 
Point4dC::operator*(const Point4dC & v) const
{ return Point4dC( X()*v.X(), Y()*v.Y(), Z()*v.Z(), T()*v.T() ); }

inline 
Point4dC 
Point4dC::operator/(const Point4dC & v) const
{ return Point4dC( X()/v.X(), Y()/v.Y(), Z()/v.Z(), T()/v.T() ); }

inline 
Point4dC 
Point4dC::operator*(const RealT alpha) const
{ return Point4dC( X() * alpha, Y() * alpha, Z() * alpha, T() * alpha); }

inline 
Point4dC 
Point4dC::operator/(const RealT alpha) const
{ return Point4dC( X() / alpha, Y() / alpha, Z() / alpha, T() / alpha); }

inline 
RealT 
Point4dC::operator&(const Point4dC & p) const
{ return X()*p.X() + Y()*p.Y() + Z()*p.Z() + T()*p.T(); }

inline 
RealT 
Point4dC::SqrEuclidDistance(const Point4dC & point) const
{
  Point4dC p(Translation(point));
  return p & p;
}

inline 
RealT 
Point4dC::EuclidDistance(const Point4dC & point) const
{ return ::Sqrt(SqrEuclidDistance(point)); }

inline 
RealT 
Point4dC::EDistanceFromOrigin() const
{ return ::Sqrt(*this & *this); }

inline 
Point4dC::CosAngleT 
Point4dC::CosAngle(const Point4dC & m, const Point4dC & b) const
{
  Point4dC relA(Translation(m));
  Point4dC relB(b.Translation(m));
  return relA & relB /relA.EDistanceFromOrigin()/relB.EDistanceFromOrigin();
}

inline 
Point4dC::AngleT 
Point4dC::Angle(const Point4dC & m, const Point4dC & b) const
{ return ::ACos(CosAngle(m,b)); }

inline 
Point4dC 
Point4dC::MidPoint(const Point4dC & point) const
{ return (*this + point)/2.0; }

inline 
Point4dC 
Point4dC::Center(const Point4dC & point) const
{ return (*this + point)/2.0; }

inline 
RealT  
Point4dC::AbsDistance(const Point4dC & point) const
{
  return   ::Abs(X()-point.X())
         + ::Abs(Y()-point.Y())
         + ::Abs(Z()-point.Z())
         + ::Abs(T()-point.T());
}

inline 
RealT 
Point4dC::Sum() const
//===================
{
  return X()+Y()+Z()+T();
}

inline 
Point4dC & 
Point4dC::Sqrt()
//==============
{
  X() = ::Sqrt(X());
  Y() = ::Sqrt(Y());
  Z() = ::Sqrt(Z());
  T() = ::Sqrt(T());
  return *this;
}

inline 
Point4dC & 
Point4dC::Abs()
{
  X() = ::Abs(X());
  Y() = ::Abs(Y());
  Z() = ::Abs(Z());
  T() = ::Abs(T());
  return *this;
}

inline 
Point4dC & 
Point4dC::Normalize()
{
  *this /= EDistanceFromOrigin();
  return *this;
}

inline 
Point4dC & 
Point4dC::NormalizeSum()
{
  *this /= Sum();
  return *this;
}

inline 
void 
Point4dC::FixAlmostZero()
{
  if (::Abs(X()) <= zeroDistance) X() = 0;
  if (::Abs(Y()) <= zeroDistance) Y() = 0;
  if (::Abs(Z()) <= zeroDistance) Z() = 0;
  if (::Abs(T()) <= zeroDistance) T() = 0;
}

inline 
Point4dC 
operator*(RealT alpha, const Point4dC & point)
{
  return Point4dC(alpha * point.X(), alpha * point.Y(),
                  alpha * point.Z(), alpha * point.T());
}

inline 
Point4dC 
operator/(RealT alpha, const Point4dC & point)
{
  return Point4dC(alpha / point.X(), alpha / point.Y(),
                  alpha / point.Z(), alpha / point.T());
}

#endif

// IAPS - Image analysis program system. 
// End of include file Point4d.hh